Characterization of the Water Supply of the Rio Grande Project Based on Rio Grande Compact Reports 1940–2020

Water sharing agreements between upstream and downstream entities are customary in many countries in the world. Issues arise when present obligations are based on hydrological conditions of the past that no longer represent the present. That has, and continues to be, the cause of disagreement and disputes between parties.

The focus of this study is the 1938 Rio Grande Compact (the Compact) water supply; however, to study it, one must learn its history as well as that of the Rio Grande Project (RGP). Fig. S1 shows the map of the geographical location encompassed by the Rio Grande Compact.

In response to claims by Mexico in 1890s of the US overdiverting water and increasing shortages in the neighboring country (American Journal of International Law 1907; Shupe et al. 1989; King and Maitland 2003), the US Department of the Interior began researching ways by which the water from the Rio Grande could be better managed (Littlefield 1999). An act of Congress in 1905 approving the construction and operation of a project called the RGP (Autobee 1994) was the result of research and a proposal from the US Bureau of Reclamation to build a dam where the current Elephant Butte Reservoir is located.

The project delivers water to two irrigation districts, one in Texas and one in New Mexico, in addition to helping fulfill the requirements of the 1906 treaty with Mexico (Littlefield 1999). The storage capacity of Elephant Butte Reservoir, about 2 million acre-feet (MAF) (about $2.5\mathrm{billion}{\mathrm{m}}^3$) buffers wide swings in Rio Grande water supply from upstream to provide a more reliable water supply for the RGP downstream. Caballo Dam, built later, is intended to provide additional regulating and flood control capability.

The Compact is an interstate agreement signed in 1938 between the states of Colorado, New Mexico, and Texas to evenly apportion the waters of the Rio Grande Basin. The Compact was ratified by the states in 1938, and by Congress in 1939, to settle disputes going back to the late 1800s (Hill 1974). A paramount purpose of the Compact was, and continues to be, to safeguard the RGP and its operations under the conditions that existed at the time the Compact was signed (Rio Grande Compact Commission 1939, p. 5, 2020).

Over the years, the three signatory states of the Compact have had disagreements (Mix et al. 2012) and from time to time sued each other, but the disagreements have been resolved. However, the current conflict began in 2011 when New Mexico sued the federal government and the two irrigation districts in the US that are part of the RGP, alleging that an agreement made in 2008 between the US Bureau of Reclamation and the RGP districts was shortcutting New Mexico of the water to which the Compact gave it a right. In 2013, Texas took the case to the Supreme Court, accusing New Mexico farmers and cities below Elephant Butte Reservoir of extracting groundwater to their detriment, thus shortening their share of the water. Given that the Rio Grande is the main water supply of the RGP and the shortage in deliveries is at the center of the current court case, this study characterized the relationship between Compact obligations and actual deliveries by doing a statistical examination of the Rio Grande Compact operations and accounting from 1940 to 2020. The authors used annual reports issued by the Compact Commission to describe the patterns of credits and debits relative to the available water supply between 1940 and 2020 to estimate when Colorado and New Mexico are more likely to have a credit or a debit.

This technical note contributes to the extensive body of literature of the Rio Grande Compact in a way that no other study has done, specifically by characterizing debits and credits for Colorado and for New Mexico. In 1974, Reynolds and Mutz (1974) studied water deliveries under the Rio Grande Compact. They plotted the deliveries of Colorado and New Mexico but did not statistically characterize the credits or debits. They pointed out something that was confirmed with this technical note, which is that the release from project storage was less than 790,000 acre-feet a year for all but 4 of the 33 years since the ratification of the Compact and that New Mexico and Colorado were in debit status most of that time. In 2011, Sandoval-Solis and McKinney (2011) conducted a risk analysis of the 1944 treaty between the United States and Mexico for the Rio Grande/Rio Bravo Basin. They contrasted the expected treaty obligations with actual historic performance and used measures of reliability, resilience, and vulnerability for their analysis. The topic and approach of their study is like this technical note in the sense that it is also an attempt to characterize differences between the expected and the actual to enlighten the path ahead. They concluded that under drought conditions, the system was not able to recover as expected, just like the one under the Rio Grande Compact, and that the deficits were bigger than what was estimated resulting in increased difficulties for signatories of the treaty as well as beneficiaries of the water agreement in both countries. In 2012, Mix et al. (2012) characterized the changes in streamflow at Lobatos gauge in Colorado to determine whether they were statistically significant between periods, but did not characterize the behavior of credits and debits of Colorado or New Mexico like this technical note does. They also conducted a water balance accounting to reveal volumetric changes to identify when those require attention. Besides the references mentioned, the authors did not find other published work that has statistically characterized the pattern of credits and debits for Colorado and New Mexico under the Rio Grande Compact of 1938, or other rivers in the US.

The study considered the entire period of record for both series, credits and debits, and then introduced qualifiers to make distinctions between hydrological conditions of an uncertain and changing environment (NMBGMR 2022). The conditions were called high (H) and low (Lo), early (E) and late (La). The (H) and (Lo) resulted from plotting a cumulative distribution for scheduled delivery for each state (SchDCO and SchDNM) to determine when natural breaks in the distribution occurred, reflecting wet periods and drought periods. The (E) and (La) indicate the period before and after the spill that occurred in 1985 (Rio Grande Compact Report 1939), which coincides with long-term conditions of the region toward a more arid climate (NMBGMR 2022). With the entire period of record and the four qualifiers, the authors calculated descriptive statistics for all the series, visualized the distribution of the series and their relationships, and conducted correlation analyses. The purpose of using the four qualifiers was to capture differences between specific hydrologic conditions of the flows and not necessarily of the system. In turn, the point of using a statistical approach is to focus on index and delivery flows to avoid the infinite complexity in the hydrology of the system.

A time series model using the autoregressive integrated moving average (ARIMA) process was performed to understand the temporal characteristics of each series of credit/debit for Colorado and New Mexico. The ARIMA model can be denoted ARIMA $(p,d,q)$, where $p$ is the order of an autoregressive model, $d$ is the degree of differencing, and $q$ is the order of the moving average model. The authors then fitted the transfer function model (TFM) to identify hydrologic factors that determine patterns of credit/debit series. The TFM combines two approaches, time series with causal, to estimate the output series from past values of more than one interconnected input time series (Hasanah et al. 2013). This study used a TFM with annual credit/debit as the output series ($Y_t$) and scheduled delivery as the input series ($X_t$). For more details, refer to Supplemental Materials. Once the models that fit the input and output series were identified, the second step was to do prewhitening for both series. Prewhitening uses ARIMA models for the input series and consists of “transformation processes which are correlated to the behavior of white noise which are not correlated” (Hasanah et al. 2013, p. 319). The prewhitening process detects the cross correlation between the input and output series using the plot.

The model selection was primarily based on the Akaike information criterion (AIC), which is a mathematical method for evaluating how well a model fits the data from which it was generated. The model with the smallest AIC value was chosen as the best model. The root-mean-square error (RMSE) and the mean absolute percentage error (MAPE) were considered to determine which model was most accurate in prediction. For details, refer to Supplemental Materials. All statistical analyses were conducted using the software R (version 4.1.0).

This technical note contributes to build the broad body of literature regarding how transfer function models have been used to date in the field of hydrology because none of the publications reviewed used a transfer function model to describe the future prediction value of the annual credit/debit time series (output series or $Y_t$) based on the past values from one or more interconnected time series (input series or $X_t$ scheduled delivery and previous year credit/debit). One documented use of TFM was in 1985 by Maidment et al. (1985), when the model was used to represent the relationship between daily rainfall and urban water use. TFM has also been used in regionalization at the finer scale for snow interception, infiltration root zone, surface runoff evapotranspiration, slow and fast interflow, baseflow, routing, and potential evapotranspiration (Samaniego et al. 2010). TFM has helped estimate travel time of recharge water through the deep vadose zone by describing the relationship between the variation in rainfall and evapotranspiration, and the changes in the level of the groundwater table (Mattern and Vanclooster 2010). In turn, Labat and Mangin (2015) applied the TFM to trace the flow of a contaminant via a tracer test in karstic systems. Also in 2015, Baillieux et al. (2015) used the TFM in the study of the evolution of groundwater quality at aquifer outlets by considering historical land use in the area under study. In the paper they document five studies that used TFM in research related to groundwater quality at an aquifer outlet. Continuing with the groundwater topic, Jimenez-Martinez et al. (2016) derived a TFM between effective rainfall and groundwater levels to predict groundwater levels for rainfall for high, medium, and low emissions scenarios from climate projections for the United Kingdom. TFM was used by Pezij et al. (2020) to characterize soil moisture dynamics, which they stated was the first time TFM was used for that purpose, with precipitation and evapotranspiration as inputs. In their paper they provide a literature reference related to the use of TFM in hydrological applications from 1980 to 2014. A study conducted by Jeong et al. (2018) used TFM for estimating water table fluctuations of unsaturated infiltration fluxes at two monitoring locations. Other examples of the use of TFM in hydrology include a comparison of the use of TFM to determine the impact of barometric pressure on well water levels (Sun and Xiang 2020), the estimation of non-point-source solute travel times along an unsaturated zone (Bancheri et al. 2021), and the quantification of the variability in discharge response of a confined aquifer due to variation in rainfall (Chang et al. 2021). There are hundreds of other articles that describe the use of TFM in the field of hydrology; however, the authors were not able to identify another study that uses TFM for the same purpose they used it.

Researchers characterized the scheduled delivery, actual delivery, and their difference, the annual credits and debits, as well as the accrued debits. Refer to Supplemental Materials for details. For details about the difference in the levels of scheduled and actual deliveries for early and for late for Colorado and New Mexico, refer to Supplemental Materials. They also characterized early and late for the credits and debits and the series for Colorado were easier to describe statistically than those for New Mexico; refer to Supplemental Materials.

To visualize the distribution of the series and their relationships for the high and low, and early and late considerations for the Colorado and New Mexico series, the authors created histograms (refer to Supplemental Materials). Continuing with the study of the behavior of annual credits and debits for Colorado and New Mexico, the authors tested if the variables were autocorrelated with past values of themselves by using the autocorrelation function (ACF) (refer to Supplemental Materials). For Colorado, the ACF revealed that the current series is positively correlated to two past series with Lag 1 and Lag 2 (Fig. S6). The correlation plot for Colorado (Fig. S7) shows that the series credit/debit (CrDebCO) is significantly correlated with SchDCO and its first lag (SchDCO.1), as well as its own two lags (CrDebCO.1 and CrDebCO.2). There is also evidence of significant correlation between Lag 1 of CrDebCO (CrDebCO.1) and lag 2 of scheduled delivery (SchDCO.2).

For New Mexico, the ACF shows that there is no evidence of autocorrelation of the series with its past values because none of the lags considered are outside of the bands (Fig. S8). The correlation plot (Fig. S9) between schedule delivery and credit/debit for New Mexico shows that the series credit/debit for New Mexico (CrDebNM) is not significantly correlated with SchDNM or its lags (SchDNM.1 and SchDNM.2), nor is it significantly correlated to its own lags (CrDebNM.1, CrDebNM.2). There is only evidence of significant correlation between scheduled delivery Lag 1 (SchDNM.1) and its Lag 2 (SchDNM.2).

Given the multiple relationships of the credit/debit and scheduled delivery time series for Colorado, time series analysis using the transfer function model was used. The credit/debit for Colorado, or output variable, was fit to the ARIMA(1,0,1). In contrast, the model for the Scheduled Delivery, or input variable, was fit with the ARIMA(0,0,0), which implies the error term is uncorrelated across time. ARIMA(1,0,1) means the credit/debit for Colorado is a function of the predictors of the lagged value of credit/debit by one period [$p=1$ in $\mathrm{ARIMA}(p,d,q)$], the model error, and lagged error term ($q=1$). Because $d=0$, a time series does not need to be differenced [$Y_t$ itself, not considering $Y_t-Y_{(t-1)}$]. ARIMA(0,0,0) means the $Y_t$ is stationary, neither autocorrelated with the previous series nor previous error terms.

The results of the prewhitening (Fig. S10) for the input and output series for Colorado indicated that the trend in credit/debit for Colorado lagged one year behind the scheduled delivery, which is consistent with the time series plots in Fig. 2. Two initial models (Model 1 and Model 2) with $b=1$ and $s=0$ were run with different parameters $r=1$ and $r=2$. The cross correlation between prewhitened input and output series determined the parameters $b=1$ and $s=0$. The $r$ value is 1 if it has an exponential plot pattern, or 2 if it has sine wave pattern (Hasanah et al. 2013, p. 320) (refer to Supplemental Materials).

Using the model selection criteria explained in the methods section and in Supplemental Materials, Model 1 was selected and defined for Colorado credit/debit series ($Y_t$) using the input series of the scheduled delivery ($X_t$)where ${\eta }_t$ has the autoregressive and moving average form with 0.8674 and $-0.6138$ of the AR(1) and MA(1) parameter estimates. For details about parameters estimated of the transfer function for Colorado, refer to Table S6.

Two additional variables were added, one at a time, to Model 1—high and low, and early and late—to determine whether there were any changes in parameter estimation. The same three model evaluation criteria were used. The parameter HighLow is statistically significant as predictor of credit/debit behavior in Colorado.

For New Mexico, the authors followed the same procedure, identification of time series model for input and output series from the autocorrelations function plot. The credit/debit for New Mexico as the output series could be fitted to the model ARIMA(0,0,0) with white noise, while the model for the scheduled delivery, or the input series, could be fitted to an ARIMA(0,0,1). ARIMA(0,0,1) means that $Y_t$ is a function of combinations of error and lagged error by one period ($q=1$), but it is not a function of past series for $Y$ because $p=0$.

Consistent with results previously obtained when examining the autocorrelation in the New Mexico credit/debit series (Figs. S8 and S11), the plot did not show a clear pattern to help determine parameters $b$, $s$, and $r$ to run the TFM; therefore, the transfer function model was not used.

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon request (table of raw data and copies of tables from Rio Grande Compact Commission reports 1940–2020).

The authors would also like to express their gratitude and recognize the support from the New Mexico Water Resources Research Institute, whose 2022 Student Water Research Grant NMWRRI-SG-2022 program provided the funding for the publication of this article.