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: