Historic Variability of the Water Inﬂow to the Lazaro Cardenas Dam and Water Allocation in the Irrigation District 017, Comarca Lagunera, Mexico

: An assemblage of tree-ring chronologies for the Nazas (NZW) in the Western Sierra Madre (WSM), Mexico was developed to determine water inﬂow at the Lazaro Cardenas Dam (LCD), the main source of water for surface irrigation in the Irrigation District 017 (DDR 017), Comarca Lagunera. A Principal Component Analysis of the ring-width chronologies was conducted to determine a common climate signal, and a stepwise model based on selected chronologies of the PC1 (CBA, COC) and PC2 (ARN) were used to develop a water inﬂow reconstruction to the Lazaro Cardenas Dam (LCD) extending from 1753 to 2003 (251 years), resulting in the following signiﬁcant ﬁndings: the warm phase of the El Niño Southern Oscillation (ENSO) in the winter-spring season had a signiﬁcant inﬂuence (SOI; Dec–Feb = − 0.24, p < 0.01), but the North American Monsoon System (NAMS) was the most important in determining the water yield in the summer season (r = 0.48, p < 0.01). Water gauge inﬂow records (77 years) at the LCD used to determine the annual allocation of water for agriculture in the irrigation district 017 was an average of 1676 × 10 6 m 3 , where the maximum annual water outﬂow allowed of 1100 × 10 6 m 3 for safety reasons, the dam infrastructure was released in 74% of the years and increasing to 78% when considering the reconstructed inﬂow. Prolonged drought episodes lasting more than 10 consecutive years were detected in the reconstructed inﬂow, information that could be used by decision makers to establish proper irrigation management strategies to ameliorate the economic and social impact when these extreme hydroclimatic events may occur.


Introduction
The Western Sierra Madre (WSM) is the most important source of water in northern Mexico, where watersheds that drain towards the Pacific Coastal Plane and the highlands of the Chihuahuan Desert yield the water needed to irrigate extensive areas located on the lowlands of the western and eastern slopes of this mountain range [1,2]. In this respect, the gauge water records in dams constitute an important source of information to determine the agricultural area and the water volumes to be allocated for irrigation in a particular year. There is significant interannual variability for most of the irrigated agricultural districts in northern Mexico. These districts depend mainly on the amount of water stored in reservoirs from the water year-streamflow (October 1st of the previous year to September 30 of the current year), which is released in the growing season for irrigation purposes. A specific 2 of 20 case is the Irrigation District 017 (DDR 017), located in the state limits of Durango and Coahuila, Mexico, which is known as "Comarca Lagunera" (Place of Lagoons). During the period of 1995-2014, this region had a mean irrigated area of 48,550 ha year− 1 , whereas the long-term average  is of 66,718 ha [3].
The annual variability in runoff volumes is influenced by ocean-atmospheric phenomena, such as El Niño Southern Oscillation (ENSO) and the North American Monsoon System (NAMS), both of which perform an important role in determining the amount of seasonal rainfall, the volume of groundwater recharged, and the amount of runoff reaching some of the main streams [4].
These phenomena (ENSO, NAMS) usually are interconnected, and their influence is enhanced when they occur simultaneously [5][6][7][8][9][10][11]. Some of the most intense droughts in Mexico, like those occurred in the 1860's and 1950's, took place when some of these circulatory patterns interacted [6]. However, further studies are needed to analyze details of this influence and other atmospheric phenomena, such as the Pacific Decadal Oscillation (PDO), Atlantic Multidecadal Oscillation (AMO), and the Artic Oscillation (AO).
Land-use changes at the upper-forested Nazas Watershed derived from intense logging activities, overgrazing, clearcutting, and human-induced fires to expand the agricultural frontier have altered the hydrologic cycle, decreased water quality, and enhanced soil erosion problems [12]. Global warming as forecasted by models developed by the Intergovernmental Panel on Climate Change [6], where northern Mexico basins are subject to more frequent and intensive droughts [7], can exacerbate these land-use changes. Dominant drier conditions may favor an increase in desertification, and, particularly, shortages in water availability for most of the human settlements established in semiarid regions [8].
The volume of water released from a reservoir depends on several factors; some are intrinsic to the dam design, but most are due to the amount of water stored during a particular water year. In this respect, in order to estimate a proper allocation of water for a sustainable use and to stabilize the irrigated area considering the interannual streamflow variability, it is important to have updated information about historical water inflow records, which are useful for a better understanding of its long-term variability, and to implement proper measurements to ameliorate the impact of extensive droughts detected in long-term streamflow reconstructions [13].
The pressure on ground water resources to fulfill the demand of a population near 1.3 million people on the Comarca Lagunera has been increased in the last decades due to a greater demand in forage production for the dairy industry, human consumption, and other industrial uses [14]. In order to mitigate this problem, the local water administration is in the process of using 200 × 10 6 m 3 of the authorized annual water outflow for irrigation to alleviate the water deficit for human consumption. Nevertheless, this action requires improvement in infrastructure and irrigation techniques in order to minimize the impact of this volume reduction and, at the same time, maintain the agricultural productivity of the region [15].
The extent of streamflow records for most of the basins in Mexico is limited in time and coverage; therefore, an alternative option to extend back in time the hydrological information is to develop dendrohydrological reconstructions from tree rings. This information could be useful to feed algorithms for modelling water allocation in a particular year, and to release an additional ecological volume to maintain the biodiversity of riparian ecosystems and to increase groundwater recharge [16].
Historical water inflow records in the "Lazaro Cardenas" Dam (LCD) have changed through time, as indicated by documented fluctuations in the agricultural area of the DDR 017. A minimum surface irrigated area of near 20,000 ha was recorded in 1952, 1953, 1957, 1958, and 1963, which is a significant decrease in the amount of land irrigated, which can be attributed to the severe droughts recorded in a broad region of northern Mexico [10,17].
The volume of the water allocated per year for most irrigated districts in northern Mexico depends on the water stored on the reservoirs in the water year, but there is a limitation on the amount and extent of the inflow gauge records. Therefore, in order to have a better understanding of the frequency and intensity of extreme hydroclimatic events and their potential impact on the availability of water resources, it is important to develop a more extended record of hydroclimatic information based on streamflow reconstructions.
The objective of this study was to develop a historic dendrohydrological reconstruction of the streamflow yield at the upper NZW, recorded as water inflow in the PLC dam, to determine the amount of water to be allocated in the growing season for agricultural purposes and compare it to the authorized water outflow by the district administrators of the DDR 017.

Study Area
The Nazas watershed (NZW) is part of the RH10 an enclosed basin of the WSM draining toward the Chihuahuan Desert. The NZW runoff is stored in the LCD, from where it is allocated for irrigation of the DDR 017. Over the last 20 years, the irrigated area in this irrigation district covers an average of 48,550 ha year −1 , where the allocated annual volume has been 775 × 10 6 m 3 . Water is yield is from the upper-forested area of the watershed, which is composed by three main sub-basins: the Sardinas, with an average annual runoff volume yield of 284 × 10 6 m 3 ; Rio Ramos, with a volume of 205 × 10 6 m 3 ; and Rio Santiago Papasquiaro and Tepehuanes, with a volume of 286 × 10 6 m 3 . The water yielded at these sub-basins represent 92% of the total inflow to the LCD [14].
The NZW has a summer precipitation regime, where 88 to 90% of the annual precipitation falls from May to September, and the remaining (12 to 10 %) occurs in the winter months. The elevation gradient (1000 to 3200 m.a.s.l) influences the amount of precipitation. At lower elevations, characterized by semiarid climate and dominant desert vegetation, annual rainfall is around 200 mm, whereas at higher elevations, with a temperate and subhumid climate and dominated by a mixed-conifer forest and oak-pine communities, annual rainfall reaches over 800 mm [1].
The interannual and multiannual precipitation variability of the NZW has a significant influence of the ENSO warm phase (Niño), particularly in the cold season [17][18][19][20]. The amount of rainfall occurring in this seasonal period significantly affects the size of the earlywood and total ring-width of conifer species present in mixed conifer forests of the WSM, favoring a greater radial growth during the warm phase, and decreasing size in the cold phase [21,22]. On the other hand, intense El Niño events favor drier conditions, decreasing precipitation and diminishing water yield in the summer period [11,23]. Largescale atmospheric phenomena influencing summer precipitation in this region, which accounts for 90% of the annual water yield, have been less studied. One specific case is the North American Monsoon System (NAMS), a circulation pattern that favors convective rains at higher-elevation portions of the western (winward) and eastern (leeward) slopes of the WSM [24]. In this system, summer rainfall is produced by an interaction between the continental warming, regional topography, and by the arrival of moisture from the Pacific and Atlantic oceans [25,26]. The variability of NAMS on inter-annual bases is more stable than ENSO variability, and this is probably the reason why there are few studies that model and forecast the behavior of this phenomenon [26][27][28].
In the Mexican Monsoon region, a portion of the rainfall occurs during a period where the ring formation is near completion. Therefore, the rainfall occurring at this seasonal period may not be fully integrated into the annual ring. Currently, this period has not been detected in the latewood portion of the annual rings of conifer species growing in the WSM [29]. In addition, the water holding capacity of the soil constitutes another constraint, since the amount of rainfall that often occurs in short periods exceeds the soil infiltration capacity of the dominant shallow soils present at this elevation. Nevertheless, summer rainfall is of great importance for rainfed agriculture, range production, groundwater recharge, and water storage purposes [2].
In the Southwestern United States, in contrast, where the influence of NAMS takes place from mid-June to late September, the latewood growth of conifer trees has been related to the amount of rainfall occurring in the summer season. Thus, research efforts have been focused on separating the climatic influence, or autocorrelation of the earlywood growth on the latewood growth, in order to capture the influence of this phenomenon [30][31][32].

Hydrometric Records
Weather data for this study was recorded from a single gauge station known as Lazaro Cardenas or Palmito (25 • 25 42 N−105 • 00 53 W; 1521 m elevation) that records the water inflow to the Lazaro Cardenas Dam. The construction of this dam began in 1936 and ended in 1946 with the objective to store the water yield at the headwater subbasins of the Nazas watershed [33]. The available water inflow and water outflow records at the LCD extended from 1947 to 2006 ( Figure 1). The annual water inflow is in average of 1676 × 10 6 m 3 and represents 92% of the water yielded at the upper watershed. In the DDR 017, the water volume entry recorded in the water year constitutes the main source to determine the irrigated agricultural area in the following growing season. constraint, since the amount of rainfall that often occurs in short periods exceeds the soil infiltration capacity of the dominant shallow soils present at this elevation. Nevertheless, summer rainfall is of great importance for rainfed agriculture, range production, groundwater recharge, and water storage purposes [2].
In the Southwestern United States, in contrast, where the influence of NAMS takes place from mid-June to late September, the latewood growth of conifer trees has been related to the amount of rainfall occurring in the summer season. Thus, research efforts have been focused on separating the climatic influence, or autocorrelation of the earlywood growth on the latewood growth, in order to capture the influence of this phenomenon [30][31][32].

Hydrometric Records
Weather data for this study was recorded from a single gauge station known as Lazaro Cardenas or Palmito (25°25′42″ N−105°00′53″ W; 1521 m elevation) that records the water inflow to the Lazaro Cardenas Dam. The construction of this dam began in 1936 and ended in 1946 with the objective to store the water yield at the headwater subbasins of the Nazas watershed [33]. The available water inflow and water outflow records at the LCD extended from 1947 to 2006 ( Figure 1). The annual water inflow is in average of 1676 × 10 6 m 3 and represents 92% of the water yielded at the upper watershed. In the DDR 017, the water volume entry recorded in the water year constitutes the main source to determine the irrigated agricultural area in the following growing season  The water volume to be allocated to the DDR 017 of a particular year is based on the volume stored at the LCD on October 1 of the previous year. The allocation of water is based on an algorithm proposed by the Mexican Institute of Water Technology (IMTA) and is based on 80 gauge records of water inflow for a policy of non-deficit water operation. However, the availability of long-term records of inflow is important to improve the algorithm used to define the amount of water outflow to be allocated for irrigation purposes of the DDR 017. The proposal indicates that, if the stored volume by the October 1 is <1650.5 × 10 6 m 3 , then the volume to be allocated could be 400 × 10 6 m 3 . On the other hand, if the stored volume is >1650.5 × 10 6 m 3 but <2355.2 × 10 6 m 3 , then the water outflow is determined by multiplying the stored volume times 0.9933 and subtracting 1296.46 × 10 6 m 3 ; a volume near 1300 × 10 6 m 3 should remain in the reservoir for safety reasons of the dam infrastructure. According to this algorithm, the maximum volume to be allocated in any case could not exceed 1100 × 10 6 m 3 , although in extreme wet years (>2355.2 × 10 6 m 3 ), some additional volume could be released to the main stream to protect the reservoir infrastructure and to favor ground water recharge.
Pluvial events with extraordinary streamflow (a standard deviation above the mean) have been recorded in the NZW in years 1968, 1974, 1991-1992, 2008, 2010, and 2016 [17]. Some of them are seen in Figure 1 and are reflected in the inflow gauge record of the same year or in the following year.

Dendrochronological Network
Hydroclimatic reconstructions for large basins usually require the development of an extensive dendrochronological assemblage of tree species sensitive to climatic variations, and when they are integrated into a representative chronology, be able to capture most of the climate variability that characterizes it [34]. The NZW and neighbor watersheds, due to their biodiversity, conifer richness, and oak species diversity constitute a region where, in the last two decades, over 50 tree-ring chronologies have been developed for climate reconstructions, and most of them have been included in the International Tree Ring Data Bank (https://www.ncei.noaa.gov/products/paleoclimatology/tree-ring, accessed on 27 May 2021). This network of tree-ring series is well distributed in the NZW and in other neighbor watersheds, but most of them are far away (over 50 km apart) from the upper-forested area of the Nazas watershed where water is yielded; therefore, in this study, most of those series were discarded for streamflow reconstruction purposes.
The site tree-ring series selected for this study (nine sites) are from two conifer species, Pseudotsuga menziesii (8 chronologies) and Picea chihuahuana (one chronology), both are climate-sensitive species with well-defined annual rings, where earlywood and latewood growth are influenced by different climatic conditions present in the previous or during the year of formation [29,35].
The sites for development the tree-ring series are located at the upper elevations of the NZW and in neighbor watersheds draining toward the Pacific Coastal Plane and to the Chihuahuan Desert ( Figure 2). The forest area located at the upper section of the NZW yield over 95% of the runoff used at the lower section of the watershed where the irrigation area is located [14]. Therefore, we focused on those site chronologies as potential predictors for inflow reconstruction purposes.
Standard dendrochronological techniques were used to develop well-dated dendrochronological series. We detrended the ring-width measurements using negative exponential curves and straight lines of a positive or negative slope to eliminate age-related trends. This procedure was conducted in the dplR library package with the R program R Core Team 2013 [Bunn 2008]. The chronology indices were computed by ratio, dividing the ring-width measurements by the values obtained from the fitted curve. The Expressed Population Signal (EPS), Rbar, Coefficient of Efficiency, and Signal-to Noise ratio were obtained from the dplR package [36].

Development of Reconstruction Model
A correlation matrix was computed for the nine ring-width chronologies to determine potential association between chronologies and to establish if there was a similar response between the ones located in the upper subbasin section of the NZW and those distributed in neighbor watersheds. The amount of variance explained by the chronologies with a common climate signal was obtained by running a Principal Component Analysis (PCA).
Factor scores from the first and second principal component, full pool (e.g., considering all nine sites) and a stepwise regression procedure were tested for inflow reconstruction. The multiple linear regression relied on stepwise procedure, which iteratively adds and removes predictors (e.g., site chronologies) in order to find out the variables in the best performing model with the higher explained variance and lowest Akaike Information Criterion (AIC) [37]. The regression models were summarized by the adjusted R 2 , and Forests 2022, 13, 2057 6 of 20 AIC; and assessed for possible multicollinearity of predictors with the variance inflation factors (VIF < 5). Residuals for all models were inspected for normality, trend and autocorrelation. Once we selected the best model, we ran a 1000 bootstrap regression in order to provide a more accurate inference among predictors. In this case, we used a bootstrap with resampling, which results in multiple replicates of the original data. This procedure is recommended when small databases are used and data does not come from a normal distribution [37].

Ocean-Atmosphere Phenomena
The influence of ocean-atmospheric phenomena (ENSO, PDO, AMO, and NAMS) on the water yield at the upper NZW was analyzed by determining the association with different indices of each phenomenon. The influence of the ENSO was analyzed by three indices, the first one involving reconstructed values of the Southern Oscillation Index (SOI_Rec) from December to February, period 1707-1977 [39]; the second one is the Multivariate ENSO Index (MEI) regarded as the most comprehensive index for monitoring ENSO since it combines analysis of multiple meteorological and oceanographic components [40]; the MEI records were obtained from NOAA (https://psl.noaa.gov/enso/mei.ext/table.ext.html, accessed on 28 May 2021); and the third ENSO index is the Tropical Rainfall Index (TRI) for the period 1894-1995, an estimate of ENSO variation using rainfall anomalies in the central Pacific that is more stable than the Tahiti-Darwin [41]. The indices of the PDO phenomenon, considered as a long-term ocean fluctuation of the Pacific Ocean and characterized by frequencies of 20 to 30 years [42], were obtained from NOAA (https://www.ncei.noaa.gov/pub/data/cmb/ersst/v5/index/ersst.v5.pdo.dat, accessed on 28 May 2021). To verify the multiple regression model, the calibration-verification procedure was conducted on the two half subperiods calibrating on one and verifying on the other subperiod and vice versa. Once we selected the best model (highest variance explained and lowest AIC), 60% of the inflow records (34 years) was used for calibration and 40% (33 years) for verification. If both subperiods were significant for the statistical tests (correlation, reduction of error, negative first difference), then the whole period of recorded data was used to build the model for reconstruction purposes.
Given the lack of precipitation records in the upper forested area where water is yielded, we used a reanalysis procedure from NASA's NLDAS-2 model to estimate monthly precipitation data [38], and to determine the relationship between both monthly precipitation and recorded annual water inflow.

Ocean-Atmosphere Phenomena
The influence of ocean-atmospheric phenomena (ENSO, PDO, AMO, and NAMS) on the water yield at the upper NZW was analyzed by determining the association with different indices of each phenomenon. The influence of the ENSO was analyzed by three indices, the first one involving reconstructed values of the Southern Oscillation Index (SOI_Rec) from December to February, period 1707-1977 [39]; the second one is the Multivariate ENSO Index (MEI) regarded as the most comprehensive index for monitoring ENSO since it combines analysis of multiple meteorological and oceanographic components [40]; the MEI records were obtained from NOAA (https://psl.noaa.gov/enso/mei.ext/table.ext.html, accessed on 28 May 2021); and the third ENSO index is the Tropical Rainfall Index (TRI) for the period 1894-1995, an estimate of ENSO variation using rainfall anomalies in the central Pacific that is more stable than the Tahiti-Darwin [41]. The indices of the PDO phenomenon, considered as a long-term ocean fluctuation of the Pacific Ocean and characterized by frequencies of 20 to 30 years [42], were obtained from NOAA (https://www.ncei.noaa.gov/pub/data/cmb/ersst/v5/index/ersst.v5.pdo.dat, accessed on 28 May 2021).

Drought Indices
The reconstructed water inflow was compared against the instrumental and reconstructed June-August Palmer Drought Severity Index (PDSI) [10]. This drought index is based in a hydrological balance for the summer months and constitutes a good measure of the available water for rainfed agriculture, grassland productivity, and filling of the reservoirs, in such a way that the water inflow to the reservoirs and PDSI values may constitute a proxy to analyze water yield fluctuations and potential trends on the watershed.
The fluctuations in water inflow variability into the LCD were also investigated using spectral analysis techniques and were used to determine which potential circulation modes exerted the greatest influence on affecting the streamflow variability. The periods in which the reconstruction shows significant frequencies was determined with a Spectral Wavelet Analysis developed in the library biwavelet in R software [43].

Quality of Tree-Ring Chronologies
The dendrochronological parameters of the tree-ring chronologies distributed on the upper-forested sub-basins and neighbor watersheds involved in assessing their quality are shown in Table 1. Based on the series intercorrelation, average mean sensitivity, Rbar, coefficient of efficiency, EPS, and SNR values, all of the chronologies are considered suitable for climate reconstruction purposes [44].

Tree-Ring Chronologies Association and Common Climatic Response
The correlation matrix shows that the association between chronologies for the common period 1773-2002 was significant (p < 0.01) for all the comparisons. Higher correlations, however, were found when comparing chronologies developed from the same species, i.e., Psedutosuga menziessi. The association was lower between the Arnulfa (ARN) chronology developed from Picea chihuahuana and the rest of P. menziesii chronologies. Results of the Principal Component Analysis indicated the presence of two principal components with the highest variance. The first principal component (PC1) loaded most of the chronologies explaining 67% of the variance, indicating that most of the chronologies are responding to a common climatic factor (Figure 3). The correlation matrix shows that the association between chronologies for the common period 1773-2002 was significant (p < 0.01) for all the comparisons. Higher correlations, however, were found when comparing chronologies developed from the same species, i.e., Psedutosuga menziessi. The association was lower between the Arnulfa (ARN) chronology developed from Picea chihuahuana and the rest of P. menziesii chronologies. Results of the Principal Component Analysis indicated the presence of two principal components with the highest variance. The first principal component (PC1) loaded most of the chronologies explaining 67% of the variance, indicating that most of the chronologies are responding to a common climatic factor (Figure 3).  The ring-width chronologies explaining greater variance in PC1 were COC, CHI, PIT, and CBA, whereas in PC2 was ARN (Figure 4).

Relationship between Dendrochronological Series and Water Inflow to the Lazaro Cardenas Dam
To assess the relationship between the tree-ring chronologiers and common variance captured by PC1 and PC2 eigenvectors and water inflow into the LCD, a Pearson correlation analysis was applied for the period 1947-2003. Correlation values were significant for the PC eigenvalues (PC1, PC2), and for the single chronologies ( Table 2).

Relationship between Dendrochronological Series and Water Inflow to the Lazaro Cardenas Dam
To assess the relationship between the tree-ring chronologiers and common variance captured by PC1 and PC2 eigenvectors and water inflow into the LCD, a Pearson correlation analysis was applied for the period 1947-2003. Correlation values were significant for the PC eigenvalues (PC1, PC2), and for the single chronologies (Table 2). The PC1 factor scores and the ARN chronology showed the highest correlation (0.75, 0.82, repectively) with water inflow records, conversely, poor associations (0.24, 0.15, and 0.29) were found for PIT, SAN and SDT, respectively, and were discarded from further analysis.
The regression models accounted for between 66% and 70% of the variance of flow, for the pooled and stepwise models, respectively. However, not all sites were significant for the full pool model, for instance, only three sites were relevant. In concordance with these results, the stepwise procedure selected the same three sites as the most relevant explaining the water inflow (ARN, CBA, and COC). The statistics for the tested models are presented in Table 3; however, the best reconstruction model was obtained with the stepwise procedure since it showed the highest adjusted R 2 and lowest AIC (Table 3).These results are confirmed with the 1000 bootstrap regression runs providing an accurate inference about predictors (complemental material). Based on these results, we decided to use the ARN, CBA, and COC chronologies for developing a reconstruction model for the water inflow to the LCD. The residual analysis indicated normality (p = 0.71) and homoscedastic error variance (p = 0.49) for the stepwise model. Furthermore, all sites included in the selected model showed a VIF lower than five (VIF < 5) which assures the no collinearity between response variables. For calibration purposes, we used subperiods 1947-1980, while the verification was from 1981-2003, and the total period 1948-2003 (Table 4). All the statistical tests for the calibration-verification procedure were significant (Table 5). 1 r = Pearson correlation coefficient (a measure of the relationship between two random quantitative variables and is independent of the measured scale). 2 RE = Reduction of Error there is no specific procedure to perform a statistic test for this parameter, but any positive result is an indication that the reconstruction contributes to unique paleoclimatic information [46]. 3 ST = Sign test relates the number of agreements and disagreements between the observed and reconstructed values and applies a statistic for the association known as t-test [47]. 4 Root mean squared error * Values are significant (p < 0.05).
Given that the total period and the half subperiods of recorded data were significant for the correlation analysis, we used the whole period  to develop a reconstruction model of the water inflow to the LCD. The dependent variable in this model was the annual water inflow to the LCD, whereas the independent variables were the annual ring-width index values of the ARN, CBA, and COC tree ring sites.
The regression equation was: where:

width index (standard version) for the ARN, CBA and COC sites, respectively
A graphical comparison between the water inflow records and the reconstructed inflow shows a good agreement between both variables with a correlation value of 0.85 (p < 0.01, n = 56). Most of the reconstructed data was in agreement with the observed inflow, but some specific years showed great differences, i.e., 1951, 1966, and 1974 ( Figure 5).

Water Inflow Reconstruction
Given that most of the statistics were significant for the verification analysis, the stepwise model developed with the gauged records 1947-2003 was used to reconstruct the water yield at the upper NZW on the total length of the chronology ( Figure 6). This information allows a more solid analysis of the water inflow to the LCD and the water outflow for irrigation purposes.
The reconstruction shows the presence of extreme hydroclimatic events affecting the amount of runoff reaching the LCD. Thus, extreme low inflow events (Mean − 1.0 S.D.) occurred from 1771 to 1773, followed by minimum inflow volumes in 1787 and 1789. The period from 1821 to 1826 had a mean inflow of 578 × 10 6 m 3 ; however, the most extended period with extreme low inflow occurred from 1951 to 1963 (12 years) with an average of 595 × 10 6 m 3 . On the other hand, high inflow events (Mean + 1.0 S.D.) were recorded in years 1754, 1759, 1791-1803, 1944, 1968, and 1990-1992. However, when we obtained descriptive statistics for 50-years subperiods, we found that the greater mean inflow occurred for the 1903-1952 subperiod and the minimum for the 1853-2002 subperiod, although the means were not statistically different. (Table 6). A graphical comparison between the water inflow records and the reconstructed inflow shows a good agreement between both variables with a correlation value of 0.85 (p < 0.01, n = 56). Most of the reconstructed data was in agreement with the observed inflow, but some specific years showed great differences, i.e., 1951, 1966, and 1974 ( Figure 5).

Water Inflow Reconstruction
Given that most of the statistics were significant for the verification analysis, the stepwise model developed with the gauged records 1947-2003 was used to reconstruct the water yield at the upper NZW on the total length of the chronology ( Figure 6). This information allows a more solid analysis of the water inflow to the LCD and the water outflow for irrigation purposes.   When an inflow analysis was conducted for 50-year subperiods, we did not find significant differences in the mean ( Table 6).
Based on the reconstructed water inflow to the LCD, and considering the algorithm proposed by the research hydrologists of IMTA to define the annual water allocation for irrigation in the Comarca Lagunera region, the dominant water outflow was <1100 × 10 6 m 3 with a probability of occurrence of 78% (36% + 42%). A close result was obtained with the recorded outflow when a water outflow <1100 × 10 6 m 3 was released in 74% (49% + 25%) of the years (Table 7).

Influence of Atmospheric Circulation Patterns
The reconstructed water inflow shows the influence of several ocean-atmosphere phenomena affecting its high-and low frequency variability. The ENSO phenomenon, throughout its diverse indices (SOI, MEI, TRI) and based on a Pearson correlation analysis, had a significant influence on the water inflow variability (−0.24, 0.25, and 0.28, respectively). Given that water inflow corresponds to the annual runoff reaching the LCD, the effect of this phenomenon at the annual level is lower as when is compared at the seasonal level, e.g., ENSO exerts a greater influence on the precipitation variability of the winterspring season. Other phenomena such as AMO, NAMS, and PDO also affect the water inflow behavior of the Nazas watershed (Table 8).

Influence of Drough Indices
The North American Monsoon System (NAMS) is the most important phenomenon affecting summer precipitation along the WSM and the one that most explains the interannual variability of the runoff volumes reaching the LCD. This finding is corroborated by a significant relationship with this phenomenon (r = 0.48, p < 0.01). On the other hand, given that the PDSI for the summer season (June-August) represents a hydrological balance for that period it constitutes a good indicator of the NAMS influence on this region [10]. The association between the inflow reconstruction and the instrumental June-August PDSI was of 0.60 (p < 0.01), resulting in greater or reduced runoff volumes reaching the reservoirs in the NZW (Figure 7).
Forests 2022, 13, x FOR PEER REVIEW 15 of 23 [10]. The association between the inflow reconstruction and the instrumental June-August PDSI was of 0.60 (p < 0.01), resulting in greater or reduced runoff volumes reaching the reservoirs in the NZW (Figure 7). The association between the monthly precipitation data obtained by the NLDAS-2 model and the reconstructed streamflow data indicated that the accumulated precipitation from July to October had a correlation value of 0.42 (p < 0. 05,, explaining 18% of the annual water inflow to the LCD.

Dominant Frequencies in the Water Inflow Reconstruction
The reconstructed water inflow for the last 251 years shows significant wavelengths influenced by dominant interannual climatic conditions at 21.2 and 64.98 years, which may be attributed to the influence of atmospheric circulatory modes (Figure 8). Relationship between the reconstructed water inflow to the LCD and reconstructed June-August PDSI values for the Nazas watershed [14].
The association between the monthly precipitation data obtained by the NLDAS-2 model and the reconstructed streamflow data indicated that the accumulated precipitation from July to October had a correlation value of 0.42 (p < 0.05, 1981-2003), explaining 18% of the annual water inflow to the LCD.

Dominant Frequencies in the Water Inflow Reconstruction
The reconstructed water inflow for the last 251 years shows significant wavelengths influenced by dominant interannual climatic conditions at 21.2 and 64.98 years, which may be attributed to the influence of atmospheric circulatory modes (Figure 8).

Discussion
The assemblage of tree-ring chronologies developed at the headwaters of the NZW and contiguous watersheds draining toward the western and eastern slopes of the WSM represent one of the most complete dendrochronological networks ever developed for a watershed in northern Mexico. However, a limitation for extracting the climatic signal of these time series and to extend climatic information back in time is the limited availability of weather and gauged streamflow records needed to establish a proper relationship between the dendrochronological series and the hydroclimatic information [21].
Even though the chronologies developed in the subbasins of the NZW showed a significant response with the water inflow records to the LCD, a higher correlation was observed when combining different chronologies present in and outside of the upperforested subbasins of the NZW. The greatest response, however, was obtained when using a stepwise regression procedure to derive a mixed model involving three of the chronologies (Arnulfa, Cerro Barajas, and Cocono) located at the headwaters of the NZW. High-elevation sites (2800 to 3000 m asl) on the western slopes of the WSM receive an annual precipitation around a 1000 mm, coming mostly on the summer months when NAMS has its greater impact on this region. In contrast, on the leeward side of the WSM, there is a decrease in precipitation due to a lower moisture content of dominant winds, which drop their moisture content on the windward side of the WSM. This condition affects forest productivity and runoff volumes draining toward the LCD that are lower in comparison to the ones recorded on the western slopes of the WSM [1].
The annual radial growth of tree species growing outside the limits of a particular watershed is often involved in dendrochronological studies to explain the behavior of climate and streamflow variability in locations far away from the place where the trees grew [50]. This response is explained by the fact that climatic factors (precipitation, evapotranspiration) limiting water availability for tree growth are the same variables controlling streamflow [51]. Streamflow reconstructions in the Colorado River basin developed by Woodhouse et al. [51] included tree-ring chronologies distributed all over the basin, but far away from gauge stations, and explained between 72 to 81% of the interannual streamflow variability. Similar studies have been conducted in basins of Canada [52], Argentina [53], Mexico [21,29,54], and other parts of the world [55,56].
In this study, of the nine chronologies involved, those with a better response with the water inflow records were Pseudotsuga menziesii and Picea chihuahuana located at the headwaters of the NZW. The development of a network of chronologies for streamflow reconstruction purposes is suitable in cases where the water yield produced in the whole watershed has a significant contribution to the streamflow volumes recorded at a gauge located at some point on the main stream. However, in semiarid watersheds where water yield come from specific locations, developing a few climate-sensitive chronologies may capture most of the climate variability required to explain hydrological records and to model their past behavior.
Given the high-climate sensitivity of P. chihuahuana to water yield in the NZW, by improving the number of chronologies of this species, we could develop a better statistically supported model to extend back in time inflow reconstructions in this watershed and to improve the results provided in this study.
The interannual and multiannual precipitation variability in basins located along he WSM is significantly influenced by ENSO, which favors a greater winter-spring rainfall in the warm phase (El Niño) and produces droughts in the cold phase (La Niña) [17][18][19]21,23,32].
Rainfall occurring in the winter-spring season favors conifers' earlywood growth and total ring-width in northern Mexico increasing biomass production [57]. However, cool season precipitation and streamflow represent less than 20% of the total annual precipitation and less than 10 % of the annual streamflow volume in the NZW. Intense Niño years, however, may produce extraordinary streamflow volumes in the winter-spring season as in years of 1958, 1968, 2008, 2010, and 1991-1992, when additional streamflow of the LCD was released for structural safety reasons and to prevent social and economic damages to settlements located downstream of the reservoir [58].
The Influence of intense ENSO events in the summer season is weak, except when this phenomenon is in phase with the PDO, thus wet conditions in summer are favored when the PDO is in its negative phase, and vice versa [5,7,59]. In this study, we did not find a significant association with ENSO indices and summer water inflow records to the LCD, except with the mean July-October PDO (r = 0.30, p < 0.05) ( Table 8).
In the irrigation area of the Comarca Lagunera, the presence of warm El Niño events during the cool season allows for an increase in irrigated area. In contrast, the extreme dry period of the 1950s reduced the agriculture area to less than 20,000 ha year −1 [3]. After the 1950s drought, there was a period when a greater water allocation was delivered to irrigate an area of near 100,000 ha, particularly in the years 1974, 1988-1989, and 1991-1992 ( Figure 1).
The NAMS is the most important ocean-atmosphere phenomenon affecting summer precipitation along the WSM in the states of Durango, Chihuahua, and Sonora, Mexico, known as the Monsoon Core Region [4,26,32]. The water yield from watersheds on this region is incorporated into streams with intermittent or perennial flows, and stored in reservoirs for agriculture, power generation, or industrial use. Alternatively, some of the streams drain freely to the ocean, creating riparian corridors of high biodiversity that demand high water volume for their ecological stability [16].
The significant relationship between the accumulated precipitation from July 5 to September 15 in the Monsoon Core Region (r = 0.48, p < 0.0005, n = 48) ( Table 8) and reconstructed water inflow to the LCD, support the high association between the runoff volumes and accumulated precipitation. This relationship is confirmed by the significant association found between reconstructed June-August PDSI values and water inflow records (r = 0.60, p < 0.01, n = 79) (Figure 7), indicating that the hydrological balance in this period constitutes a significant parameter of the annual volume of streamflow that could be recorded in the LCD as inflow.
Although the inflow reconstruction was statistically well supported, some particular years showed significant differences between gauged and reconstructed water inflow; thus, 1974 had a recorded inflow of 3259 × 10 6 m 3 , whereas the reconstructed values was 1781 × 10 6 m 3 , a difference of 1478 × 10 6 m 3 (83%). In contrast, other years showed greater reconstructed values than the recorded inflow, such as 1966 with a reconstructed inflow of 2475 × 10 6 m 3 and a recorded value of 1277 × 10 6 m 3 , a difference of 1197 × 10 6 m 3 (93.8%) ( Figure 5).
The lack of agreement between recorded and reconstructed inflow for specific years may be explained by heavy rains exceeding the soil infiltration capacity, when rainfall occurs in a period of tree inactivity, or at the end of the growing season [46]. In the second case, when reconstructed inflow is greater than the recorded, previous soil moisture conditions or rainfall occurring right in the growing season, even if is not of great magnitude, it may be sufficient to influence a greater radial growth, conditions that may distort the relationship between both variables [46].
According to the Power Spectrum Analysis, significant peaks on the water inflow reconstruction were detected at cycles of 21 and 65 years, similar results to the ones detected in a streamflow reconstruction for the Nazas watershed, and which correspond to the 50-to-100-year wavelengths reported for north-central Mexico by Stahle et al. [60]. Dominant wavelengths of 20 to 30 years may be related to the PDO influence [48], given that a low but significant correlation was found between reconstructed water inflow records and PDO indices for the NZW.
Wavelet Power Spectral Analysis showed significant periods less than 10 years along the reconstructed period, falling in the ENSO domain ( Figure 9). The ENSO influence on streamflow reconstructions has often been reported for other streamflow reconstructions worldwide. Thus, significant peaks of 3, 6,22,47,115,200,260, and 300 years were found for a streamflow reconstruction of the Yellow River in China [61]. Frequencies of 22, 110, and 115 years were reported for streamflow reconstruction in North America [62].
The influence of ENSO (SOI, TRI, MEI) on the reconstructed annual inflow volume had correlation values in the range of 0.24 to 0.28 (p < 0.01), particularly in the cool season; this significant but low correlation is due to the fact that most of the annual inflow in this and neighbor watersheds is yielded in the summer season as influenced by NAMS characterized by low interannual variability, and, therefore, it is difficult to model their behavior [27]. In comparison, rainfall and streamflow variability is greater in the winter-spring season as modulated by ENSO conditions that has some potential to be forecasted [63]. In this study, long-term streamflow reconstructions are relevant for a better understanding of the impact of atmospheric circulation modes, to analyze their influence on the interannual and multiannual climate variability, to determine potential trends of water availability, and to use that information for establishing sustainable management strategies of water resources.
Inflow volumes <2355.2 × 10 6 m 3 implying outflow volumes <1100 × 10 6 m 3 have dominated the allocated water outflow for irrigation in the DDR 01. This volume, however, is more persistent when considering the water inflow reconstruction, which provides a more extended record. These findings suggest the need for developing a management plan strategy for drought periods as long as 12 consecutive years with water yields <595 × 10 6 m 3 . These extended dry periods have already been documented in this water inflow reconstruction and may re-occur in the next few decades. Based on this information, it is important to analyze in detail the influence of the NAMS, the main circulatory phenomenon affecting the summer water yield in watersheds located in the WSM.
Climate change models from the IPCC suggest an increase in aridity for northern Mexico in the next few decades [6], conditions that may already be in the process of reducing the water availability for irrigation and other uses [64] The models suggest drier than normal conditions in dry years and increases in overland flow in wetter years, which may be enhanced by deforestation and other land-use changes taking place on this region [1,12]. Conservation strategies focused on the upper forested subbasins where water is yield in the NZW and based on the payment of environmental services to the owners of those areas should be important to implement to ameliorate the potential impacts forecasted by climate change models.

Conclusions
Developing an assemblage of tree-ring chronologies outside the limits of a particular watershed is important to capture the regional climate signal and to establish a significant relationship with gauge streamflow records for hydroclimatic reconstructions. In this study, selected ring-width series from trees located at the upper-forested section of the NZW and neighbor watersheds had a common climatic signal and were used to analyze the historical water inflow to the LCD.
The water inflow reconstruction represents a significant contribution to understand the interannual and long-term variability of the runoff reaching the LCD, the main reservoir that captures the streamflow yield at the upper NZW used for irrigation in the DDR017.
Ocean-atmospheric phenomena had a significant influence on the amount of rainfall and streamflow yield at the upper-forested NZW. The warm phase of ENSO increased precipitation and runoff in the cool season, particularly for some of the wettest years detected in the recorded and reconstructed data. Winter-spring precipitation is around 10% of the annual rainfall, but significantly correlates with the annual radial increase in conifer species in this watershed.
The influence of ENSO in summer precipitation was non-significant, but in extreme ENSO events, its influence was extended to the summer season, producing extreme droughts or pluvials. In the WSM, where the main moisture source comes from NAMS, summer precipitation represents over 90% of the annual rainfall and 92% of the runoff in watersheds. Other ocean-atmospheric phenomena, e.g., AMO and the PDO, may have a synergic effect on summer rainfall when they are in phase with ENSO.
Historical water allocation based on gauge records had a dominant range of <1100 × 10 6 m 3 in 74% of the years examined (water inflow < 2355 × 10 6 m 3 ), and this percentage increased to 78% of the years examined when considering the reconstructed water inflow for the last 251 years. Water administrators of the DDR 017 should take proper measures to anticipate and ameliorate droughts lasting over 10 consecutive years to plan for suitable water management strategies that could prevent the social and economic calamities that these events provoke.