Hydrological Variability in the El Cielo Biosphere Reserve, Mexico: A Watershed-Scale Analysis Using Tree-Ring Records

Abstract

The El Cielo Biosphere Reserve (CBR) stands as a vital forested region in eastern Mexico due to its high biodiversity in flora and fauna and provision of environmental services. This study established a network of 10 ring-width chronologies of different species within the CBR and adjacent watersheds. The objective was to analyze their climatic response and reconstruct the seasonal streamflow contribution of each sub-basin to the main stream, utilizing data from a gauge network of eight hydrological stations located at strategic locations of the CBR. With chronologies ranging from 116 to 564 years, most exhibited association with the accumulated streamflow between January and June. Based on the adjusted R², Akaike Information Criteria, and Variance Inflation Factor, the stepwise regression procedure was selected among different statistical methods for developing the reconstruction model. In spite of differences in the seasonal reconstructed periods, all the species showed potential to develop hydrological reconstructions as indicated by their common response to streamflow variability, as occurred in the wet years of 1976, 1993, 2000, and 2008, and dry years of 1980, 1982, 1996, and 2011. It was found that the response of the chronologies to gauge records increased as a function of the chronologies’ interseries correlation, average mean sensitivity, and distance of the tree-ring series to the gauge station. Streamflow reconstructions at the sub-basin level allowed a better understanding of the hydroclimatic variability characterizing the CBR, but also suggested the need to increase the network of chronologies for some particular sub-basins lacking tree-ring series to improve the reconstructed models.

1. Introduction

Hydrological resources are some of the most essential environmental services provided by forest ecosystems [1]; however, in the last few decades, an increase in deforestation, overgrazing, logging, and other land-use changes, along with climate warming, has disrupted the hydrological cycle in most of the ecosystems worldwide [2,3]. These changes affect the amount and distribution of rainfall, increasing evapotranspiration rates and decreasing water quality for different human uses [4,5].

The El Cielo Biosphere Reserve (CBR) is northeastern Mexico’s most important protected area. It was created in 1985 by the State government of Tamaulipas, Mexico, and recognized as a Biosphere Reserve by the United Nations in 1987 [6]. The CBR covers an area of 144,530 ha, and it is recognized for its high biodiversity in flora and fauna richness favored by the interaction of Neotropical and Nearctic ecosystems [7].

The importance of this protected area involves different ecological and socioeconomic issues, given that rural communities settled on the reserve made their way of living based on ecotourism services. However, ecotourism could have negative impacts on the reserve’s biodiversity if the design of interfaces that promote knowledge exchange among different sectors, systems, and stakeholders is not considered [8].

The CBR has been the focus of a fair amount of ecological research studies, mainly biodiversity [9,10]; however, no studies are available to analyze the impact of climate variability on the water resources yielded in this watershed that could explain much of the biodiversity of the area and the use of water for productive activities.

At the same time, a long-term understanding of climate variability in the CBR is constrained by the availability of climatic and gauge records that are scarce for northeastern Mexico, short in length (usually <50 years), poor quality, and representativeness. For example, only one gauge station, located within the CBR limits, has records from the middle of the 20th century. An alternative approach to obtain climate and hydrological information back in time is the use of proxy records, where tree rings constitute one of the most suitable sources of information given their annual resolution and the potential to develop centuries-long time series that are helpful in reconstructing climate and streamflow volumes and determining their interannual and long-term variability, as well as potential trends [11].

Streamflow reconstructions from tree rings are essential to analyze periods larger than the available gauge records, to evaluate the amount of flow produced at a watershed level, and to determine the frequency and trends of extreme hydroclimatic events, the influence of major circulatory modes on climate variability, and the impact they have on the proper planning and management of available water resources, mainly in dry regions [12,13].

Long-term information on water availability helps in designing strategies for proper ecosystem management, preserving biodiversity, and developing conservation strategies for the forested areas with the highest water yield, mainly where downstream settlements depend on the water produced at the watershed headwaters.

The objective of this study was to reconstruct the streamflow variability of the El Cielo Biosphere Reserve at the basin level using a set of tree-ring-width chronologies. For this, we did not restrict our focus to the boundaries of the CBR, but instead included information from all sub-basins that, at some point, influenced the hydrological behavior of the CBR. Finally, the variables we examined were distance to the gauge station, tree species, and statistic quality of the tree-ring series, to determine their relationship with the strength of the correlation between streamflow and ring-width indices. The relationship between tree rings and streamflow is essential to analyze their interannual variability and the potential trends as affected by major atmospheric circulatory modes and climate change.

2. Materials and Methods

2.1. Site Description

The El Cielo Biosphere Reserve (CBR) is part of the Sierra Madre Oriental (Eastern Sierra Madre), located in the southeastern portion of the State of Tamaulipas, in the municipalities of Gomez Farias, Llera, Jaumave, and Ocampo in northeastern Mexico. Geographically, it encompasses the extreme coordinates from 22.925° to 23.431° N and −99.097° to −99.441° W, at an elevation that ranges from 200 to 2320 m a.s.l. (Figure 1).

Geologically, limestone rocks of sedimentary origin dominate the region and have formed a karstic topography with rocky outcrops, sinkholes, and caves. Dominant soils are lithosols, morphologically characterized by their poor development, low fertility, and limited water-holding capacity [14]. Most rainfall infiltrated appears at lower elevations as springs, representing the main water source for the Rio Sabinas and Rio Frio.

The CBR is part of a larger basin known as “The Bajo Rio Panuco” (Lower Panuco River), Hydrological Region 26 (RH-26), that covers an area of 15,256.5 km², where the main stream Rio Panuco (Panuco River) drains a mean annual volume of circa 200,000 × 106 m³ into the Gulf of Mexico. The CBR encompasses four sub-basins: Guayalejo, Sabinas, Comandante, and Drenaje Subterraneo, covering an area ranging from 160 to 5519 km². Notably, the CBR is characterized by its significant hydrological potential, contributing 72% of the water flow in the basin. This substantial water input originates from the Sabinas, Frio, and Comandante rivers and a portion of the Guayalejo River, the primary collector for the entire basin.

The climate of the reserve has different variations according to the elevation, with a dominance of warm–subhumid conditions at elevations from 300 to 800 m, where the lowest temperature for the coldest month is >18 °C with an annual precipitation range from 1500 to 2000 mm. Elevations between 700 and 2000 m dominate a semi-warm humid climate with a temperature of the coldest month between 15 and 18 °C and annual precipitation between 1000 and 1500 mm [15]. Most of the rainfall (55 to 65% of the annual total) occurs in the summer months (June to September), with the influence of the Bermuda–Azores (a large subtropical semi-permanent center of high atmospheric pressure favoring tropical cyclogenesis) [16]. In contrast, the rest of the precipitation occurs in the winter months, influenced by the entrance of cold fronts [17].

The variability in environmental conditions, topography, and elevation characterizing the reserve favors the presence of over ten vegetation types [18]; among them are tropical deciduous and semi-deciduous forest, cloud forest, mixed-conifer forest, oak–pine forest, desert shrublands, grasslands, and aquatic and subaquatic vegetation, where riparian vegetation is included [7].

2.2. Tree-Ring Chronologies

We employed ten dendrochronological sites distributed within and outside the CBR limits, considering the sub-basin distribution (Figure 1). Accordingly, we used tree-ring measurements of Montezuma bald cypress (Taxodium mucronatum Ten., HUA and RSA sites), vejar fir (Abies vejarii Martínez, LEN site), nut pine (Pinus cembroides Zucc., MIP and JAU sites), Montezuma pine (Pinus montezumae Lamb., RDC site), nelson pine (Pinus nelsonii Shaw, MIQ site), weeping pinyon pine (Pinus pinceana Gordon, MAT site), Mexican mountain pine (Pinus hartwegii Lindl., MAR site), and Douglas fir (Pseudotsuga menziesii Mirb. Franco, PEN site) (Figure 1). Most of our sites were developed by the authors, except for one being an updated version of the ring-width database developed by Stahle and collaborators for the Rio Sabinas [19] https://www.ncei.noaa.gov/access/paleo-search/study/4935 (accessed on 16 March 2024). To update this site, we sampled 66 bald cypress trees located along the Rio Sabinas stream banks (23.148° N, −99.15° W; 317 m a.s.l.). At all sites, we collected two to three increment cores per tree with a Pressler increment borer, in the opposite direction to the stem and avoiding sampling rotten wood.

We processed the cores from the rest of the sites following standard dendrochronological techniques [20]. The dating of the raw measurement file was conducted with the COFECHA program [21], and for detrending purposes, we used dplR software version 3.4.3 [22,23]. A ring-width index (RWI) was calculated by dividing the raw measurements with a fitted value according to a smoothing spline function with a frequency response of 0.50 at a wavelength of 0.67 series length. To build a site chronology, we averaged the RWI using Tukey’s biweight robust mean [11]. Finally, to analyze the quality of the chronologies, we used the standard version of the chronology and considered different parameters such as series intercorrelation, average mean sensitivity, mean interseries correlation (Rbar) [24], and Expressed Population Signal (EPS) [25].

2.3. Streamflow Data and Imputation Methods

Gauge data were acquired through the National Water Commission (CONAGUA) hydrometric database, which provided access to historical records from gauge stations located in the sub-basins of interest [26]. We compiled data from eight gauge stations in the studied area (

Table 1. Gauge stations located within or near the El Cielo Biosphere Reserve (CBR). The stations are arranged from low to high average streamflow.

No.	Gauge Identification	Sub-Basin	Information PERIOD	Latitude (N)	Longitude (W)
1	Zapata	Rio Guayalejo	1987–2014	23.2714°	99.5417°
2	Servilleta	Rio Comandante	1960–2014	22.8536°	99.1211°
3	Ahualulco	Rio Comandante	1945–2014	22.9869°	99.1489°
4	La Encantada	Rio Guayalejo	1960–2014	23.3855°	99.0791°
5	Sabinas	Rio Sabinas	1960–2014	23.0331°	99.0942°
6	Frio	Rio Tamesi	1960–2013	22.8431°	99.0261°
7	San Gabriel	Rio Guayalejo	1984–2014	23.0853°	98.7878°
8	Tamesi	Rio Comandante	1973–2014	22.8942°	99.4319°

).

The La Encantada, Sabinas, and La Servilleta gauge stations showed no missed data during the recorded periods. The Rio Frio and San Gabriel stations had less than 1% and 4% data gaps, respectively. The Tamesi station exhibited an 11% data gap for 2006–2010, while the Zapata station had a 15% absence in daily records for 1995–1998 and 2010–2013. To address this issue, we employed the imputTS package from R to impute missing data. This package offers several imputation methods suitable for univariate time series. The “na.interpolation” method with the “stine” option was selected due to its effectiveness in adjusting the curvature of the interpolation based on the slope of adjacent points, thereby providing more accurate estimations than simpler methods. Theoretical details and the implementation of this method are explained in the imputeTS manual [27].

2.4. Association between Ring-Width Chronologies and Streamflow Data and Reconstruction Models

The ring-width chronologies of the CBR were compared against accumulated monthly streamflow records from September of the previous year to October of the current year of growth. This comparison was conducted using a bootstrap Pearson correlation with the treeclim package in R [28], where we correlated all tree-ring sites with gauge records to determine the optimal period of association between variables.

We accumulated streamflow records of consecutive months with a significant association, and then compared them to the ring-width chronology to develop a seasonal reconstruction model for each gauge station. After identifying the optimal period for each station, we tested single dendrochronological sites, the first component of all dendrochronological sites (PC1), and stepwise regressions using site-level information for runoff reconstruction. We selected the best model according to the adjusted R², Akaike Information Criteria (AIC), and a Variance Inflation Factor (VIF) lower than ten, when appropriate. Residuals for all models were inspected for normality, trends, and autocorrelation. Afterwards, the gauged streamflow records were divided into two subperiods, one for calibration and one for verification, and if reconstruction statistics were adequate (reduction of error—RE; sign test—ST; correlation coefficient—r), the whole period of recorded data was used to build the final model.

For the final model, we followed the approach described by [13], using a set of 1000 bootstrap regressions. The bootstrap approach allows more accurate information to be obtained among predictors and is particularly useful when small databases are used and do not come from a normal distribution. In this case, we used bootstrapping with resampling, which results in multiple replicates of the original data [29]. This process was performed for each gauge station.

2.5. Influence of Explanatory Variables in the Strength Signal of Tree Rings and Streamflow Data

To assess the impact of various explanatory variables on the strength of the correlation between tree rings and streamflow data, we performed an analysis of covariance (ANCOVA) in R, version 4.1.0 [22]. This analysis considered categorical variables (tree species and land use) as well as continuous variables (distance between tree-ring sites and gauge stations, elevation, stream order, intercorrelation, and sensitivity).

To explore the common variability and detect periods of low and high streamflow at each gauge station, yearly records for both observed (annual cumulative runoff) and reconstructed periods were transformed into z-scores. This was conducted by subtracting the mean value from each record and dividing it by the standard deviation, allowing for comparison among gauge stations even if the streamflow differed. Furthermore, we used a principal component analysis to determine a common behavior among gauge reconstruction runoffs, using the common period from 1902 to 1995.

3. Results

3.1. Length and Quality of the Ring-Width Chronologies

The dendrochronological series spanned from 1450 to 2021, covering a 572-year period. The PEN site had the longest record (564 years, Pseudotsuga menziesii), while JAU had the shortest (116 years, Pinus cembroides) (Figure 2). All sites exhibited high interseries correlation values, ranging from 0.42 to 0.75, exceeding the minimum 0.328 threshold to be considered properly dated [21]. Additionally, all sites had high rbar and EPS values, indicating suitable characteristics of the tree-ring data. Spatially, only two sites were located at lower elevations, near the bottom of the sub-basins, while the rest were at higher elevations (

Table 2. Dendrochronological sites within the influence area of the El Cielo Biosphere Reserve (CBR). Sites are listed from the highest to the lowest intercorrelation.

Site	Species	Elevation m a.s.l.	Period	n (cores)	Series Intercorrelation	Average Mean Sensitivity	rbar	EPS
MIQ	Pinus nelsonii	1950	1846–2019	87	0.75	0.50	0.59	0.99
MIP	Pinus cembroides	2080	1882–2019	114	0.72	0.42	0.53	0.99
LEN	Abies vejarii	2700	1894–2021	295	0.69	0.34	0.51	0.99
MAR	Pinus hartwegii	3300	1838–2019	90	0.68	0.41	0.50	0.99
PEN	Pseudotsuga menziesii	2700	1450–2013	108	0.66	0.31	0.48	0.99
RDC	Pinus montezumae	1800	1772–1995	27	0.61	0.43	0.42	0.95
JAU	Pinus cembroides	1760	1902–2017	69	0.59	0.33	0.37	0.98
MAT	Pinus pinceana	1070	1810–2019	52	0.57	0.56	0.35	0.96
RSA	Taxodium mucronatum	317	1474–2017	118	0.44	0.44	0.19	0.97
HUA	Taxodium mucronatum	400	1810–2019	68	0.42	0.39	0.19	0.94

).

3.2. Daily and Monthly Runoff for the Gauge Stations

The gauge stations Gabriel and Tamesi recorded the highest streamflow (m³ day⁻¹), while Zapata and Servilleta reported the lowest. Despite this, the region experienced common extreme events, such as a notorious peak observed in several stations in 1976 (Figure 3). Most sites begin their highest runoff season in June, at the onset of the rainy season, concluding by October, except for Zapata, which extends until November. We detected critical years with low annual runoff in 1973, 1993, 2000, and 2014. In contrast, years with high annual runoff were identified in 1974, 1991, and 1992.

3.3. Response of the Chronology to Monthly Streamflow Records

A comparison between the ring-width indices and monthly streamflow from the previous September to the current October indicates significant correlations (p < 0.05) across all gauge stations (Figure 4). However, some dendrochronological sites were less correlated with the gauge data. For instance, the Zapata station exhibited the highest correlation values (up to 0.62 in April, MIQ site), followed by Gabriel (r = 0.55 in May, HUA site) and Sabinas (r = 0.54 in April, PEN site). Generally, the April to June period had better correlations for most of the tree-ring sites. However, specific sites were also correlated with previous conditions (previous November to current May, LEN, MAR, and PEN with the Ahualulco station). Interestingly, the sign of the association changes from a positive correlation to a negative one starting in August, although it is not statistically significant.

3.4. Development of the Streamflow Models for Reconstruction Purposes

Based on the information from Figure 4, we selected the best correlated period at each gauge station for reconstruction purposes. Consequently, we chose to accumulate runoff from January to June for the Zapata station, February to April for Servilleta, January to April for Ahualulco, January to May for Encantada, January to June for Sabinas, February to June for Frio, April to May for Gabriel, and April to June for Tamesi.

After testing all models, they demonstrated statistical superiority over each station’s pooled version, principal component model, or single version (

Table 3. Statistical comparison among different models explaining the runoff within or near the El Cielo Biosphere Reserve (CBR). Note that some statistics are given in intervals for single tree species because they represent values across sites. The symbol (⸸) indicates the best model for each gauge station based on R² and AIC.

Gauge Station	Period	Model	Variables	p-Value	R2	R2adjusted	AIC
Zapata	January to June	Single tree species	HUA to MIQ	[0.02, 0.96]	[0.02, 0.18]	[0.01, 0.14]	[940.2, 945.6]
		Principal component model	PC1	0.54	0.02	0.01	945.6
		Pooled model	HUA to MIQ	0.45	0.94	0.56	327.7
		⸸Stepwise model	HUA, LEN, JAU, RSA	0.02	0.93	0.85	109.2
Servilleta	February to April	Single tree species	HUA to MIQ	[0.001, 0.44]	[0.02, 0.23]	[0.01, 0.22]	[1352.5, 1362.2]
		Principal component model	PC1	0.007	0.14	0.13	1355.6
		Pooled model	HUA to MIQ	0.08	0.48	0.25	852.5
		⸸Stepwise model	JAU, MAR, RSA	0.002	0.43	0.36	831.4
Ahualulco	January to April	Single tree species	HUA to MIQ	[0.01, 0.96]	[0.02, 0.13]	[0.01, 0.12]	[2468.3, 2478.3]
		Principal component model	PC1	0.05	0.05	0.04	2474.4
		Pooled model	HUA to MIQ	0.08	0.30	0.15	1792.5
		⸸Stepwise model	HUA, MAR, PEN, RSA	<0.001	0.28	0.22	1778.9
Encantada	January to May	Single tree species	HUA to MIQ	[0.001, 0.20]	[0.02, 0.23]	[0.02, 0.22]	[2435, 2450.4]
		Principal component model	PC1	0.002	0.14	0.12	2442.3
		Pooled model	HUA to MIQ	<0.001	0.51	0.39	1731.7
		⸸Stepwise model	JAU, MAR, MAT, PEN, MIQ, RSA	<0.001	0.51	0.43	1722.4
Sabinas	January to June	Single tree species	HUA to MIQ	[0.001, 0.96]	[0.01, 0.21]	[0.01, 0.20]	[2160.5, 2164.1]
		Principal component model	PC1	0.008	0.12	0.11	2160
		Pooled model	HUA to MIQ	0.003	0.60	0.44	1407.3
		⸸Stepwise model	HUA, MAR, MAT, MIP, MIQ, RSA, RDC	<0.001	0.56	0.45	1398.2
Frio	February to June	Single tree species	HUA to MIQ	[0.001,0.66]	[0.03, 0.13]	[0.02, 0.11]	[2207.6, 2212.4]
		Principal component model	PC1	0.036	0.08	0.06	2209.7
		Pooled model	HUA to MIQ	<0.01	0.55	0.39	1450.2
		⸸Stepwise model	HUA, MAT, MIP, MIQ, RSA, RDC	<0.001	0.53	0.44	1439.9
Gabriel	April to May	Single tree species	HUA to MIQ	[0.001,0.96]	[0.02,0.28]	[0.01,0.27]	[1223.4, 1229.4]
		Principal component model	PC1	0.03	0.15	0.12	1225.4
		Pooled model	HUA to MIQ	0.08	0.47	0.24	1239.9
		⸸Stepwise model	JAU, MAR, RSA	0.003	0.40	0.33	1220.1
Tamesi	April to June	Single tree species	HUA to MIQ	[0.02, 0.58]	[0.01, 0.17]	[0.01, 0.14]	[1353.7, 1358.6]
		Principal component model	PC1	0.03	0.14	0.11	1354.8
		Pooled model	HUA to MIQ	0.33	0.47	0.11	1003.6
		⸸Stepwise model	JAU, LEN, MAR, MIP	0.05	0.39	0.25	978.2

). Additionally, these models exhibited homoscedasticity, normally distributed residuals, and VIF < 10 (Table S1 in Supplementary Materials). As mentioned earlier, some gauge stations were more suitable for reconstruction purposes, as indicated by high R² and R²_adjusted values (up to 0.93 and 0.85, respectively, for the Zapata station) and the lowest AIC. In contrast, Ahualulco exhibited the lowest performance, with R² and R²_adjusted values of 0.28 and 0.22, respectively (Table 3).

All stepwise models exhibited satisfactory reconstruction skills, as indicated by positive values for RE, ST, and R² during both the calibration and verification phases. For instance, the Pearson correlation coefficient ranged from 0.46 (Ahualulco) to 0.86 (Sabinas). Additionally, the RE ranged from 0.211 to 0.742, indicating that the performance exceeds that of simply using the mean of the records (Table S2 in Supplementary Materials). Thus, we used the whole period at each gauge station to conduct a set of 1000 replications. The dependent variable was the runoff in the selected period, ranging from January to June for the more extended season (Zapata and Sabinas) to the shorter period from April to May (Gabriel). The independent variables included the sites and, consequently, the species-level chronologies. The statistics reveal that tree-ring information explained the variance in the runoff in CBR, ranging from 21% to 64% (see

Table 4. Bootstrap performance of selected periods at each gauge station. The statistics represent the median for 1000 replicates. See Table 2 for a complete reference of site names.

St	Pe	Final Model	R2	R2adjusted	AIC
Za 1	J-J	$Y_i={-2.97\times 10 }^8 {+1.46\times 10 }^7\left(HUA\right){+3.97\times 10 }^6\left(JAU\right) {-6.48\times 10 }^6\left(LEN\right){-4.01\times 10 }^6\left(RSA\right)$	0.39	0.28	934.4
Se 2	F-A	$Y_i={-3.84\times 10 }^5 {-2.78\times 10 }^5\left(JAU\right){+3.32\times 10 }^5\left(MAR\right) {+5.85\times 10 }^5\left(RSA\right)$	0.36	0.32	1341.6
Ah 3	J-A	$Y_i={2.88\times 10 }^7 {-1.64\times 10 }^7\left(HUA\right){-8.70\times 10 }^6\left(MAR\right) {{+3.51\times 10 }^7\left(PEN\right)+4.50\times 10 }^6\left(RSA\right)$	0.21	0.16	2465.1
En 4	J-M	$Y_i={7.10\times 10 }^6 {+1.85\times 10 }^7\left(JAU\right){-2.93\times 10 }^7\left(MAR\right) {+3.44\times 10 }^7\left(MAT\right) {-3.13\times 10 }^7\left(MIQ\right){{+4.94\times 10 }^7\left(PEN\right)+3.30\times 10 }^7\left(RSA\right)$	0.43	0.37	2421.2
Sa 5	J-J	$Y_i={-1.44\times 10 }^8 {+1.74\times 10 }^8\left(HUA\right){-7.18\times 10 }^7\left(MAR\right) {+5.52\times 10 }^7\left(MAT\right) {+9.08\times 10 }^7\left(MIP\right){-5.85\times 10 }^7\left(MIQ\right)$ ${-1.05\times 10 }^8\left(PEN\right)+{1.72\times 10 }^8\left(RSA\right)$	0.64	0.55	1386.9
Fr 6	F-J	$Y_i={-1.68\times 10 }^8 {+3.69\times 10 }^8\left(HUA\right){+1.05\times 10 }^8\left(MAT\right) {+1.51\times 10 }^8\left(MIP\right) {-1.48\times 10 }^8\left(MIQ\right){{-2.60\times 10 }^8\left(RDC\right)+1.54\times 10 }^8\left(RSA\right)$	0.60	0.52	1430.1
Ga 7	A-M	$Y_i={-1.25\times 10 }^8 {-1.29\times 10 }^8\left(JAU\right){+1.14\times 10 }^8\left(MAR\right) {+2.50\times 10 }^8\left(RSA\right)$	0.41	0.35	1213.9
Ta 8	A-J	$Y_i={2.45\times 10 }^8 {-1.06\times 10 }^9\left(JAU\right){-6.83\times 10 }^8\left(LEN\right) {+1.06\times 10 }^9\left(MAR\right) {+7.84\times 10 }^8\left(MIP\right)$	0.34	0.25	1350.3

St: station; Pe: period; ¹: Zapata (January–June); ²: Servilleta (February–April); ³: Ahualulco (January–April); ⁴: Encantada (January–May); ⁵: Sabinas (January–June); ⁶: Frio (February–June); ⁷: Gabriel (April–May); ⁸: Tamesi (April–June).

).

Our results demonstrate that, irrespective of high or low runoff volumes, all gauge stations were correlated to tree-ring information, enabling the reconstruction of runoff events beyond the available records. In some cases, specific sites exhibited the capability to detect both extremely high runoff and years with low runoff, as exemplified by the Servilleta and Gabriel gauge stations (Figure 5b,g). Although some tree-ring chronologies spanned centuries (Table 2), merging them with shorter yet highly sensitive sites resulted in the oldest reconstructed runoff dating back to the mid-19th century, while the shortest spanned from the beginning of the 20th century onwards (see Figure 5). Furthermore, in specific instances, a lack of tree-ring information was observed despite the availability of gauge data, highlighting the importance of updated information, particularly in Sabinas and Frio (see Figure 5e,f).

3.5. The Common Variability within the El Cielo Biosphere Reserve

Overall, all gauge stations exhibited a common variability in their behavior when considering the annual cumulative runoff. Consequently, years with high runoff volumes, such as 1976, 1993, 2000, and 2008, were observed. Similarly, years with low runoff were concurrent in 1980, 1982, 1996, and 2011 (see Figure 6a). However, it is important to note that in recent years, lower synchrony has been observed, particularly with some stations like Zapata.

The common variability in their reconstructed versions showed mixed results, with few common periods concurrent among all gauge stations. For instance, the period 1907–1912 was observed in almost all stations except for Zapata, Frio, and Sabinas. Similar results were observed around 1920, 1940, and 1975. This apparent mismatch may originate from the fact that the reconstructions represent different periods of the year (e.g., months). For example, while Sabinas represents the accumulated runoff for six months, the Gabriel station represents only two months (Figure 6b). This difference was mirrored in two different groups according to a principal component analysis, with the first group involving Ahualulco, Tamesi, Gabriel, and Servilleta, and the second group comprising Encantada, Sabinas, Frio, and Zapata (See Figure 6c). However, the higher synchrony at the yearly level implies a common influence of climatic drivers at the watershed scale, such as those related to ocean–atmospheric modes.

3.6. The Effect of Species Composition on the Strength of Correlations

The variables influencing the strength of the correlation between tree-ring series and gauge records involved the tree species, series intercorrelation, average mean sensitivity, and the distance to the gauge stations (

Table 5. Analysis of covariance (ANCOVA) to detect variables influencing the correlation between tree-ring series and gauge data.

Variable	F-Value	p-Value
Tree species	8.33	<0.001
Elevation	0.75	0.387
Series intercorrelation	8.95	0.002
Average mean sensitivity	25.78	<0.001
Stream order	0.01	0.905
Distance to gauge station	52.56	<0.001

). Notably, Taxodium mucronatum, Pseudotsuga menziesii, and Pinus hartwegii exhibited higher correlations, whereas Pinus montezumae and Pinus nelsonii showed the lowest correlations. Regarding dendrochronological statistics, higher intercorrelation values—indicating similar climatic responses among trees within a site—and average mean sensitivity were associated with statistically significant correlations with runoff data. Similarly, a shorter distance between tree-ring data and gauge stations resulted in a higher correlation, highlighting the importance of having tree-ring information close to each one of the gauge stations. Notably, the amount of runoff, which can be indirectly related to the stream order, does not influence the strength of the correlation. Therefore, any station, regardless of its position in the watershed (highlands or outlets), has the potential to be related to the tree-ring records.

4. Discussion

4.1. The Importance of Dendrochronological Networks for Hydrometric Reconstructions

This study has explored runoff variability in the CBR and adjacent sub-basins using tree-ring chronologies (Table 2). Our findings reveal significant patterns of hydrological variability, reflecting the complex interaction between climate and water resources in this region. Despite the spatial location of the gauge stations within the sub-basins, all were capable of being related to tree-ring records, albeit reflecting different seasons throughout the year. While some stations showed correlations for six months (e.g., the Zapata and Sabina stations), others only reflected a limited period over the year, such as Gabriel station, but with a strong correlation with the tree-ring data (Figure 4).

Research in dendrohydrology has been distinguished by implementing increasingly extensive networks of moisture-sensitive trees [30,31]. These networks have been important in developing proxy records that reflect historical regional hydroclimate and interannual variability and trends. A notable example is a series of studies that have reconstructed the paleohydrology of the Colorado River [32,33,34,35]. These studies have enriched our understanding due to the addition of new tree-ring chronologies from various sites in the river basin. Moreover, applying innovative statistical methods to tree-ring data and hydroclimatic information has enabled more profound and detailed insights. Here, we demonstrate the importance of certain tree-ring series that are highly sensitive to runoff records, such as RSA (i.e., Taxodium mucronatum, Table 4), which must be regularly updated to enhance our capability to understand the hydrological variability of critical regions such as the CBR. Furthermore, our results also highlighted the significance of the distance of tree-ring data from the gauge location. Thus, selecting sensitive tree species closer to the gauge locations may ensure a higher chance of reconstructing the historical runoff of a sub-basin; nevertheless, this approach will be based on a statistical relationship between secondary growth and runoff excluding physiological processes, whose understanding may contribute to explaining the influence of internal factors that trigger the size of the annual ring [36].

In the context of Mexico, the studies conducted by Osorio et al. [37] have been fundamental in understanding climate dynamics over several centuries. These researchers identified the recurrent presence of droughts at the beginning and at the end of each century, from the 16th to the 21st, based on a historical reconstruction of precipitation in the CBR from tree rings of Taxodium mucronatum. Furthermore, they highlighted the occurrence of cyclical droughts approximately every 50 to 60 years. Other findings with moisture-sensitive species distributed in this region like Pseudotsuga menziesii and Pinus cembroides indicate a winter–spring (January–April or January–May) precipitation response [38,39], similar to the one detected with Taxodium mucronatum in the CBR, which implies that they have a greater response to the cool season precipitation extending from the previous October to the current May or June. When comparing these findings with our results, we observed a synchrony between the periods of low precipitation (Table S3 in Supplementary Materials) and the reduction in runoff at the same time intervals, i.e., 1870s, 1900s, 1930s, 1950s, 1980s, 2010s (Figure 5). This pattern reinforces the idea that variations in precipitation, particularly during winter and spring, directly and significantly impact runoff levels in the CBR, underscoring the importance of these studies for water resource management and environmental planning in the region.

With climate change, precipitation patterns are expected to become more erratic, potentially leading to greater variability in runoff [40]. This could mean more prolonged periods of drought interspersed with extreme precipitation events. The alteration in runoff patterns would have significant implications for water resource management [41]. A more erratic water supply could affect water availability for agriculture, human consumption, and biodiversity conservation, especially in regions that depend on natural sources like the CBR.

Changes in runoff will also influence aquatic and terrestrial ecosystems. Species that depend on a specific water regime could face elevated risks, potentially leading to changes in species abundance, distribution, and possibly local extinctions [40].

4.2. The Runoff Variability beyond the Hydrometric Records in the El Cielo Biosphere Reserve

The correlation between tree rings and gauge records reveals that other environmental factors may become limiting in years with higher moisture, as noted by [11]. This dynamic could explain why certain high-flow events are not fully recorded in the tree-ring series, particularly occurring during the rainy season, which coincides with the highest runoff season (Figure S1 in Supplementary Materials). It is probable that these may have occurred during periods when trees are physiologically inactive or in soils with a low water-holding capacity or low infiltration rate. When developing an annual streamflow reconstruction, tree growth is limited during the winter season, making it difficult to capture the entire streamflow behavior in the reconstruction [42,43]. However, our results showed that specific sites had the ability to detect extreme positive flow peaks, as well as reflect years of low runoff. We observed this behavior at stations like Servilleta and Gabriel, where the recording capability of tree-ring chronologies is notable and coincides with the highest association between runoff and precipitation (Figure 5b,g, Table S3 in Supplementary Materials). Nevertheless, it is clear that even with 1000 replications through the bootstrap procedure, there were years with overestimated/underestimated runoff levels, which can limit management strategies in the region. This could be addressed by establishing new sampling sites close to the gauge stations and using moisture-sensitive species such as Taxodium mucronatum or Pseudotsuga menziesii. Additionally, exploring different detrending methods, such as signal-free or Bayesian detrending, could maximize the climatic signal contained in tree rings [44].

4.3. Outlook in the El Cielo Biosphere Reserve

Our findings significantly impact the management of water resources in the El Cielo Biosphere Reserve. Understanding long-term hydrological variability can assist in water conservation planning and adaptation to future climate changes, especially considering the climate change scenarios in precipitation. Pichardo et al. [45] forecasted a decrease in annual precipitation in the Guayalejo River sub-basin, which could result in modifications to the water balance.

A potential limitation of our study could be the geographic representativeness of the tree-ring samples. Future research could expand the sampling network to include more sub-basins, tree species, and sampling sites near gauge stations, thus improving the accuracy of our reconstructions. Additionally, integrating more climatic and environmental variables into the flow reconstruction models, as Goeking and Tarboton [5] suggested, would allow for a better understanding of the factors influencing hydrological variability, including common climatic conditions (mainly precipitation and evapotranspiration), which control tree growth and streamflow reconstruction [46,47].

Consequently, dendrochronology is an instrumental methodology for reconstructing historical climatic and hydrological patterns, providing indispensable insights for current decision-making processes and future strategic planning. Its implementation within the CBR, as well as in analogous regions, could play a fundamental role in formulating water management strategies that are both more effective and sustainable, adapted to the unique challenges and environmental conditions inherent to each locality.

5. Conclusions

Most of the published papers on dendrohydrological reconstruction are based on a representative network of tree-ring chronologies for a whole basin, but in this study, we emphasize the importance of considering partial reconstructions at the sub-basin level, as independent contributors to the total streamflow of the basin. This finding is important when the tree-ring series developed from different species correspond to a common period of streamflow record, which in this case corresponded to April, May, and June.

The use of different approaches (stepwise, bootstrap) allowed us to find a proper reconstruction model explaining over 50% of the seasonal streamflow variability for some single sub-basins.

In order to improve the chronologies’ response to streamflow records, most of them should be developed closer to the gauge stations, and these data should be complemented with hydrological and ecological information of the sub-basin for a better interpretation of the results.