Methodology
Study
site
The research was carried out in five headwater micro catchments located
between 1241 m a.s.l. and 1713 m a.s.l. in the TMCF zone in central
Veracruz, Mexico (Figure 1). A detailed description of the
characteristics of each micro catchment is provided in Table 1. The
micro catchments are drained by first or second-order perennial streams
and are located within the catchments of the Pixquiac and Gavilanes
rivers (106 and 42 km2, respectively), which are, in
turn, sub-catchments of the Antigua River basin (1565
km2). The micro-catchments where chosen based on the
dominance of five LULC categories of interest: primary forest (PF; 100%
of total cover; Table 1), intermediate age secondary forest (IF; 77%),
young secondary forest (YF; 68%), shade coffee plantations (SC; 94%),
and high-intensity pasture (IP; 63%). Together these land uses comprise
89% of LULC in the Pixquiac and Gavilanes sub-watersheds (Von
Thaden-Ugalde, unpublished data). The most important land use change in
the last 40 years in the Pixquiac and Gavilanes catchments is TMCF
regeneration. For example, in the Pixquiac catchment the forest cover
increased from 6561 ha in 1975 to 7685 ha by 2004 (Paré & Gerez, 2012).
Recent land use maps have confirmed that the forest cover has remained
constant during the last decade and accounts for around 79 % of the
study area (CONANP et al., 2015).
Figure 1. Location of the study area in central Veracruz, Mexico, and
maps of the study catchments showing the instrumentation and
experimental sites. Sources: (INEGI; 2013). See the text for catchment
abbreviations.
The PF micro-catchment is dominated by old-growth TMCF (>
50 years old) with low disturbance; less than 10% of the area is
covered by a 15-year-old secondary forest. The IF site is mainly covered
by 40-year-old TMCFs (>77 % of the area). The remaining
proportion of the IF micro-catchment is pastureland and annual crops
(mainly maize) located in the upper parts. The YF micro-catchment is
covered (> 68% of the area) by 20-year-old TMCF, with the
remaining area of pasture and crops (Von-Thaden U., unpublished data).
dominated by Quercus spp, Oreomunnea mexicana, Turpinia insignis,
Liquidambar styraxiflua, Carpinus tropicalis, Clethra macrophylla(Williams-Linera & Vizcaíno-Bravo, 2016, Garcı́a-Franco et al., 2008),
with greater heterogeneity in the IF and YF micro-catchments. The mean
tree height in the three catchments is between 20-25 m, with a few
larger trees reaching more than 30 m in height.
The original vegetation in the IP micro-catchment was TMCF, which was
cleared more than 40 years ago (Paré & Gerez, 2012; and local
inhabitants, personal communication). Since then, the pasture has been
heavily grazed by sheep, goats and cows (>63% of the land
area), with grass height generally less than 5 cm. In addition,
approximately 10% of IP is covered by cropland (maize and, more
recently, potatoes) and associated shrub-dominated fallows. Also,
pastures within IP feature compacted cow trails running roughly along
the contours, which likely results in reduced infiltration in those
areas and promotes Hortonian overland flow during intense rainfall
events (Bruijnzeel et al., 2006; Tobón et al., 2010). The SC
micro-catchment has been covered mainly by shaded coffee for more than
80 years (Marín-Castro et al., 2016), with 94 % of the land area
dominated by this land cover. This production system retains some trees
to provide shade to the coffee and reduce water stress. However, coffee
cultivation in the area includes management practices such as removal of
the herbaceous groundcover, pruning, fertilization, and agrochemical
applications, which may affect runoff generation processes.
soils in the micro-catchments are classified as Umbric Andosols derived
from volcanic ash, with clay and silty clay as dominant textures
(Campos-Cascaredo, 2010; Paré & Gerez, 2012). Topsoils in all
micro-catchments are characterized by low bulk densities (<
0.7 g cm-3) due to the abundance of poorly crystalline
materials and organic matter. In general, soils in TMCFs exhibit lower
bulk densities in comparison to pasture and coffee, revealing less soil
compaction (Looker N., unpublished data; Table 1). In addition, soil
profiles are generally deep (A + B horizons >1 m and C + Cr
horizons > 10 m on ridges and backslopes) and moderately
well developed (Karlsen, 2010) in all micro-catchments, favoring good
water storage. Moreover, the soils in the region are generally
underlined by andesitic saprolite, with high permeability ranging from
0.05 to 0.08 mm h−1 (Gabrielli & McDonnell, 2011;
Karlsen, 2010; Muñoz-Villers & McDonnell, 2012). Although we did not
measure bedrock hydraulic properties here, we observed the presence of
saturated saprolite on various road cuts in our study sites.
The general climate is classified as humid temperate with abundant rains
during the summer (Koppen classification modified by García, 2004). The
average temperature and precipitation between 1000 and 2000 m asl in
this area ranges between 13 to 21 °C and 1500 to 3000 mm (Shinbrot et
al., in prep.; SMN, 2015; Williams-Linera & Vizcaíno-Bravo, 2016).
Approximately 80% of the rainfall falls as convective storms during the
wet season (May–October), when the region is under the influence of the
easterly trade wind flow. During the dry season (November–April)
rainfall is generally associated with cold fronts and characterized by
light rains and drizzle (Holwerda et al., 2010).
Table 1. Topographic, land cover and soil physical characteristics of
the five study catchments. Where available, the standard deviation (SD)
is provided. Mean Kfs is the geometric mean
because this variable is log normally distributed, whereas for the other
variables the mean is an arithmetic mean.
rainfall, P (top y axis; grey bars) and streamflow,Q (bottom y axis; blue lines), as measured at the five
study catchments from Feb 26, 2015 to Aug 23, 2019. See the text for
catchment abbreviations.
Hydrometeorological
measurements
Rainfall was measured with climate stations located in clearings
(> 30 m from forest edge) in each of the five study
micro-catchments (Figure 1). The stations were equipped with tipping
buckets rain gauges, with dataloggers (Campbell Scientific and Davis),
Campbells had 0.1 mm resolution and Davis stations had 0.2 mm, recording
readings every 10 and 15 minutes (Figure 2).
Streamflow was measured using V-notch weirs at the catchment outlets
(90° angle for the YF, IF, PF, IP and 120° for the SC catchment). The
120° angle was preferred in SC because this catchment has a larger area
and the design peak discharge was larger. Water levels were registered
every 1.5 min using Solinst water level sensors paired with barometric
pressure recorders. Water levels were converted to streamflow (L
s-1) using experimental stage-discharge relationships
for these weirs calibrated with field-derived rating curves generated
via volumetric and salt dilution measurements of discharge (ASTM
D5242-92, 2013; Moore, 2005).
Rainfall measurements in the five catchments started in February 2015,
whereas discharge data collection began in July 2016. The monitoring of
discharge continued only to July 2018 in the SC and YF. For PF, IF and
IP, rainfall and streamflow measurements continued until August
2019 (Figure 2). Damage to the
equipment by animals and problems with access to the sites result in
periodic gaps in the streamflow and rainfall data; on average, these
gaps account for around 25 % of the period of record, for both
variables.
A wide range of hydrological conditions were observed in each
micro-catchment allowing for standardized comparisons of metrics and
statistics across the land covers (Nainar et al., 2018; Ochoa-Tocachi et
al., 2016). The streamflow data was normalized by the contributing
catchment area (Figure 2). In mountainous regions, Zhang et al., (2011)
and Kienzle (2010) demonstrated that the true surface area should always
be considered when more than 30% of the studied area has slopes steeper
than 18.2 degrees. In our study area, the PF micro-catchment has a mean
slope around 20 degrees. Therefore, we standardized our information
using true surface areas; even though the other catchments have more
gentle slopes. Surface areas were estimated using a triangulated
irregular network (TIN) model, created by transforming a 15-m DEM
(INEGI) into a continuous 3-dimensional surface, according to the method
proposed by Jenness (2004) and Zhang et al., (2011). This procedure was
conducted using 3D Analyst Tools in ArcMap version 10.5.1.
Data process and
analysis
Streamflow and rainfall data were resampled (mean streamflow and sum of
rainfall) to daily timesteps and a set of metrics were calculated for
each catchment to compare their hydrologic regimes. Table 2 presents
definitions for each of these metrics. To estimate rainfall to runoff
ratios, only days containing paired values of rainfall and discharge
were selected, while all the available data was used to explore
variability of streamflow and rainfall. The following parameters were
computed: average daily rainfall (MP), average daily discharge (MQ), and
rainfall to runoff ratio (RR) (Muñoz-Villers & McDonnell, 2013; Sawicz
et al., 2011). All the available data of rainfall was used to estimate
the annual ratio of days with zero precipitation (DAYP0), and daily
rainfall variability (PVAR). Likewise, the complete dataset of discharge
was used to compute the daily flow variability (QVAR), flow duration
curve (FDC), slope of the flow duration curve (SFDC), mean annual high
flow (MAHF), and mean annual low flow (MALF)
(Kelleher et al., 2015; Nainar et al., 2018; Ochoa-Tocachi et al., 2016;
Ogden et al., 2013; Price et al., 2011; Tetzlaff & Soulsby, 2008).
Additionally, the baseflow recession constant \((k)\) was obtained from
the master recession curve (MRC), which was constructed using the
matching strip method (Toebes & Strang, 1964), during the period where
the minimum rainfall inputs were identified in the dry season (Chapman,
1999); a large k value indicates slow drainage and greater
storage capacity (Murphy & Stallard, 2012).
Table 2. Definition of the hydrological indices analyzed in the study
conducted in central Veracruz, Mexico.
To assess stream response to precipitation and the influence of
antecedent wetness conditions on runoff response in the studied
catchments, 266 storm events were examined during the period from July
1st, 2016 to August 10th, 2018. For
this analysis, ten-minute rainfall and streamflow data were used to
graphically separate streamflow into quickflow (direct flow in response
to a rainfall event) and baseflow (the delayed flow from storage)
following the approach of Hewlett and Hibbert (1967). The hydrograph
separation was performed using the slope constant method; further
details can be found in Mosley (1979). Storms were defined as periods
with more than 0.2 mm of rainfall, separated by a dry period of at least
3 h (Muñoz-Villers & McDonnell, 2013).
For each storm event, the following parameters were calculated: total
rainfall (Pev [mm]), total runoff
(Qt [mm]), quickflow (Qqf [mm]),
baseflow (Qbf [mm]), peak discharge
(Qpeak [mm h−1) and the 24 h
antecedent precipitation (AP [mm]). Furthermore, we calculated the
lag time (TL [h]), defined as the time between peak
rainfall and peak discharge, and the time to peak (Tp
[h]), defined as the time between the onset of storm discharge and
peak discharge (Mosley, 1979).
Statistical
methods
Statistical tests were applied to detect differences in the hydrologic
responses across the investigated land cover catchments. Since the
distribution of events did not follow the normality assumption
(Shapiro-Wilk test), we chose the non-parametric Kruskal-Wallis H test
followed by the post-hoc Dunn’s test (α = 0.05) (Nainar et al., 2018).
The latter test may be understood as a test for median difference, which
makes multiple pairwise comparisons based on approximations to the
actual rank statistics. These statistics were carried out in the R
software v.3.5.3 (Figure 4).
The storm parameters (Qt, Qqf,
Qbf, Qpeak) were positively correlated
with the total event rainfall (Pev). Consequently, to
control for the effects of this covariate, two approaches were applied.
First, the parameters were normalized by rainfall event
(Qt/Pev,
Qqf/Pev,
Qbf/Pev) and compared with the
non-parametric methods listed above. Second, a multiple regression
analysis was conducted on each hydrologic variable versus
Pev across the five land covers, followed by an analysis
of covariance (ANCOVA) of the resulting linear models in R software
v.3.5.3. (Baumer et al., 2017; Nainar et al., 2018; Staelens et al.,
2008; Verma et al., 2018).
In the ANCOVA we modeled the land cover-specific relationship between
Pev and each of the storm parameters using a multiple
linear regression model (Equation 1) with a distinct slope and intercept
for each land cover: