The Use of Historical Data and Global Climate Models to Assess Historical and Future Surface Water and Groundwater Availability in the Trinity River Basin in Texas

THE TEXAS WATER JOURNAL is an online, peer-reviewed, and indexed journal devoted to the timely consideration of Texas water resources management, research, and policy issues. The journal provides in-depth analysis of Texas water resources management and policies from a multidisciplinary perspective that integrates science, engineering, law, planning, and other disciplines. It also provides updates on key state legislation and policy changes by Texas administrative agencies. For more information on on the Texas Water Journal as well as our policies and submission guidelines, please visit texaswaterjournal.org. As a 501(c)(3) nonprofit organization, the Texas Water Journal needs your support to provide Texas with an open-accessed, peer-reviewed publication that focuses on Texas water. Please consider donating.


INTRODUCTION
In 2019, the U.S. Congress provided the U.S. Geological Survey (USGS) Water Availability and Use Science Program with resources to implement Integrated Water Availability Assessments (IWAAs). The purposes of the IWAAs are to provide nationally consistent assessments of water availability and identify factors that limit water availability. The IWAAs are designed to meet the following six objectives: (1) provide accurate assessments of available water resources, (2) determine the quantity of water available for human and ecological needs, (3) quantify long-term trends in water availability, (4) provide assessments of changes in water availability, (5) explore factors that limit water availability, and (6) forecast water availability for economic development, energy production or conservation, and environmental or other in-stream uses (USGS 2021a). The Trinity River Basin (TRB) in Texas was selected as one of 10 basins to support development of IWAAs with state and local partners using cooperative matching funds (USGS 2021b).
The TRB is a major source of water for large, rapidly growing metropolitan areas in Texas. Maintaining the quantity and quality of water resources is vital to meeting the water demands of the greater Dallas-Fort Worth metropolitan area and downstream metropolitan areas such as the Houston metropolitan area, which are among the fastest growing cities in the United States (U.S. Census Bureau 2020a, 2020b. In response to current (2022) water demands and increasing demand requirements associated with population growth, an improved understanding of current and future water availability is needed to help resource managers efficiently manage water resources in the TRB and plan for future needs. To gain a better understanding of water availability, the USGS entered into a cooperative agreement with the following entities that manage water resources in the TRB: Trinity River Authority (TRA), City of Dallas, North Texas Municipal Water District, and Tarrant Regional Water District.
Texas Water Journal, Volume 14, Number 1 The Use of Historical Data and Global Climate Models to Assess Historical and Future 36

Purpose and Scope
The purpose of this study was to assess TRB water availability, quantify long-term trends, and forecast future TRB water availability. This was done using surface-water and groundwater models to evaluate future conditions under different global climate scenarios. This assessment builds on a recently completed study that evaluated historical long-term trends in streamflow and other hydrologic properties for the TRB as well as other water-supply basins in Texas . In this study, those historical long-term trends were used to forecast climate scenarios, and the results were used to inform a surface-water model (Precipitation-Runoff Modeling System [PRMS] model) and groundwater model (MODFLOW numerical groundwater-flow model). The water-availability assessment was done by using historical and recently collected data to assess possible changes in future TRB water availability. The surface-water and groundwater models were used to evaluate future conditions under different existing global climate models (GCMs) through 2099 and 2087, respectively. Using historical data and GCMs provided a means for comparison of various climate projections in the TRB.
The GCMs that were used were published as part of the Localized Constructed Analogs (LOCA) downscaling Coupled Model Intercomparison Project Phase 5 (CMIP5; LLNL 2021). To help account for the inherent uncertainty associated with GCMs, 30 different GCMs were chosen to simulate the ranges of possible values for climatic variables such as precipitation and temperature for surface-water and groundwater predictive models. An ensemble mean was computed for each climatic variable from these ranges of possible values.

Study Area
The Trinity River's headwaters are in north-central Texas, west of the Dallas-Fort Worth metropolitan area. From there, the river flows approximately 550 miles southeast into the Gulf of Mexico, east of Houston, Texas (TWDB 2019). South of the Dallas-Fort Worth metropolitan area, four main tributaries join to form the main stem of the Trinity River: Clear Fork Trinity River, West Fork Trinity River, Elm Fork Trinity River, and East Fork Trinity River (Figure 1; TRA 2021). The TRB covers 17,913 square miles and is the largest river basin contained entirely in Texas (TWDB 2019). The Texas Water Development Board's (TWDB) regional water planning groups (RWPGs) divide Texas into different regions for water-management purposes (TRA 2021). Most of the TRB (81%) is included in either RWPG C (hereinafter referred to as Region C) or RWPG H (hereinafter referred to as Region H). Region C includes the upstream part of the TRB, whereas Region H includes the downstream part of the TRB. The Trinity River provides water to more than half the population of Texas. As of July 2016, the populations of regions C and H were 7.23 and 6.80 million, respectively, and are projected to increase to 14.0 and 11.7 million, respectively, by 2070 (TRA 2021). Municipal water demands accounted for 90% of the total use in Region C in 2016and 55% in Region H in 2015(TRA 2021. About 90% of the water supply in Region C is from surface water, mostly from reservoirs, and about 71% of the water supply in Region H is from surface water (TRA 2021).
There are 32 reservoirs in the TRB, with a total of about 7.0 million acre-feet of conservation storage (TRA 2021;TWDB 2019). As of 2022, the USGS operates 24 lake and reservoir water-surface elevation stations in the TRB (USGS 2022). Of the 32 reservoirs, 14 were analyzed in Harwell et al. (2020) for historical long-term trends, and these 14 reservoirs represent 74% of the total storage in the TRB.  divided the TRB into five sections; section 1 was the most downstream, and section 5 was the most upstream section. Within Harwell et al.'s (2020) five defined sections, the mean annual precipitation from 1900 through 2017 was 51. 28, 44.87, 38.55, 37.38, and 34.53 inches for sections 1, 2, 3, 4, and 5, respectively ( Figure 2). Annual precipitation and annual reservoir surface evaporation indicate that for most of Texas, evaporation exceeded the mean annual precipitation: annual evaporation during 1954-2013 averaged 55.1 inches, while the mean annual precipitation during 1940-2014 averaged 39.4 inches (Wurbs and Zang 2014;Wurbs 2021;TWDB 2021a).
Although surface water accounts for 90% of Region C's water supply, groundwater is an important source of municipal water supply in some of Region C's rural areas (TRA 2021). Within Region H, groundwater accounts for about 28% of water supply (TRA 2021). One groundwater source, the Trinity River alluvium aquifer (TRAA), underlies the Trinity River and its adjacent stream corridor and tributaries ( Figure 1). The TRAA consists of alluvium and terrace deposits of gravel, sand, silt, and clay of Quaternary age (Hanko and Brikowski 2009;USGS 2014) and covers about 5,265 square miles, or 29% of the TRB. Groundwater-surface-water interactions take place between the alluvium and terrace geologic units that contain the TRAA and the overlying Trinity River and its tributaries, indicating that groundwater resources can be affected by changes in streamflow. Although TWDB does not recognize the TRAA as a major or minor aquifer in Texas, it is recognized as a viable aquifer by Groundwater Management Area 14, Region H, and the Bluebonnet Groundwater Conservation District (Groundwater Management Area 14 et al. 2016;Region H Water Planning Group et al. 2020;Williams 2010). Although few data related to the TRAA were available as of 2021 for model input and calibration (TWDB 2021b), modeling of the groundwater in the alluvium was included in this study to better understand future water availability.  Figure 1. Trinity River Basin study area, Trinity River alluvium aquifer extent, and Regional Water Planning Group extents.

Figure 1.
Trinity River Basin study area, Trinity River Alluvium aquifer (TRAA) extent, and regional water planning group extents.  (Hay 2019) and showing sections used in long-term trend analysis . Surface Water and Groundwater Availability in the Trinity River Basin in Texas

METHODS
The following methods were used to better understand water availability: a statistical analysis of hydrologic trends and the development of surface-water and groundwater models for forecasting purposes. A surface-water model of the TRB was developed by using the PRMS (Leavesley et al. 1983). A numerical groundwater-flow model of the TRAA was developed by using MODFLOW-NWT (Niswonger et al. 2011), which is a modified form of the MODFLOW-2005 groundwater-flow model (Harbaugh 2005). The surface-water and groundwater models were then used for forecasting climate scenarios. Historical precipitation and air temperature values, as well as GCMs, were used to simulate current and projected water-budget components. The methods for each approach and the methodology for selecting the GCMs are discussed below.

Long-term Trends
Using the same methods used in Harwell et al. (2020), Kendall's tau was used to detect upward or downward trends in precipitation, groundwater levels, and streamflow. Kendall's tau is a rank-based correlation coefficient that measures the strength of the monotonic relationship between two variables. The relationships between precipitation and streamflow and between streamflow and storage were also assessed using Kendall's tau. Multiple regression equations with periodic functions were developed to test the statistical significance of any changes in annual mean air temperature over time at the 95% confidence level (p-value ≤ 0.05; Helsel et al. 2020).
The approach used by Harwell et al. (2020) to analyze for trends in different sections of the TRB was also used in this study. The previously mentioned five sections defined by Harwell et al. (2020) were created by making roughly equal divisions of the aggregated counties that overlap the TRB. These sections were used for analyzing historical long-term trends ( Figure 2). This accounted for possible latitudinal and longitudinal climate differences across the TRB with respect to precipitation trends. An area-weighted daily mean precipitation total was computed for each section from the daily precipitation data (NOAA 2021). The area-weighted daily mean precipitation totals were analyzed for temporal trends in precipitation over monthly, seasonal, and annual time steps ). The trends from Harwell et al. (2020) informed scenarios for the surface-water and groundwater modeling.
Where Harwell et al. (2020) detected statistically significant monotonic historical long-term trends in annual or seasonal precipitation, this study used the Theil slope estimate to calculate the seasonal or annual change in precipitation quantities to estimate additional future reservoir water volume gains or losses. The Theil slope, a nonparametric estimate of a regression slope (Helsel et al. 2020), was also used to quantify projected water-budget components from modeling results.
Nonparametric tests were used to facilitate statistical comparisons of datasets that might differ from a normal distribution (Helsel et al. 2020). The Mann-Whitney rank sum test was used to test for differences in decadal data. The Mann-Whitney test indicates whether one group-or decade in this case-tends to produce larger observations than a second group (Helsel et al. 2020). No assumptions were made about the distributions of the data in either group.

Surface-water Model: Precipitation-Runoff Modeling System
PRMS is a deterministic, distributed-parameter, physical-process-based modeling system developed to evaluate watershed-scale hydrologic responses to various combinations of climate variables (Markstrom et al. 2015). Hydrologic simulations are done on a daily time step with daily input data, and the model outputs are designated by the user as daily, monthly, or annual time steps.

Model Inputs and Parameters
Hydrologic simulations for the TRB were done by using past and projected climate data as inputs to the PRMS model. Climate input variables-daily precipitation, daily minimum air temperature, and daily maximum air temperature-were used in the hydrologic simulations to evaluate annual hydrologic response to changes in climate variables from 2018 to 2099. Outputs of the model include annual water budget variables such as precipitation, actual evapotranspiration, surface runoff, and groundwater recharge. In the PRMS model, surface runoff is generated when the precipitation rate exceeds the infiltration rate of soil that may not be saturated; this type of surface runoff is referred to as "Hortonian runoff" (Horton 1933).
To account for the complexity of the hydrologic cycle, several inputs and parameters are required to compute hydrologic simulations in PRMS. A full list of required components is provided in Markstrom et al. (2015). In conjunction with the PRMS software, the surface-water model used in this study includes the USGS National Hydrologic Model (NHM) data infrastructure, which is designed "to fill the gap between the detailed local models used in engineering hydrology and global land-surface models." (p. 193 in Regan et al. 2019). The NHM data infrastructure was configured for use with PRMS and allows for extraction of one or more watersheds or hydrologic response unit (HRU). This data infrastructure provides a standardized modeling platform for model distribution, comparability, and interoperability; a consistent geospatial structure; and default parameter values (Regan et al. 2018).
NHM-PRMS can be used to compute hydrologic-simulation results of the temporal and spatial distribution of water availability and storage across the continental Unites States using national-scale datasets. These datasets include hydrography, The Use of Historical Data and Global Climate Models to Assess Historical and Future 40 solar radiation, evapotranspiration, geology, soils, land cover, topography, snow-covered area, and snow-water equivalent. NHM-PRMS is currently configured with the Daymet Version 3 dataset (Thornton et al. 2016), which includes daily values of precipitation, minimum air temperature, and maximum air temperature from January 1980 through December 2016. A national PRMS model, referred to as the "NHM-PRMS, by HRU Calibrated Version" was recently published (Hay 2019) and was used in this study. This surface-water model includes PRMS parameter values calculated using the calibration procedure from the NHM-PRMS, by HRU Calibrated Version. The multiple-objective, stepwise, automated calibration procedure was used to identify the optimal set of parameters for each HRU using historical climate data from 1980 to 2016.
The historical climate inputs used in the surface-water model weredaily precipitation, daily minimum air temperature, and daily maximum air temperature on a sub-watershed or HRU scale. As shown in Figure 2, the model consists of 1,192 HRUs and 620 stream segments and has an area of 11,471,544 acres (Yesildirek et al. 2023). The model also includes daily streamflow inputs of 71 streamgaging stations from 1980 through 2016, selected for the national PRMS model. Parameters comprise values calculated geospatially by HRU or stream segment with a monthly time step or for the duration of the model-simulation period, depending on the parameter type.
PRMS hydrologic outputs are computed using methods based on physical laws and/or empirical relations. Climate outputs for the model are calculated by HRU and then geographically weighted within the model to provide an average for the TRB. The actual evapotranspiration output in the model is the computed rate of water loss, which reflects the availability of water to satisfy potential evapotranspiration; specifics of the computation are presented in Leavesley et al. (1983). To calculate Hortonian runoff, the "srunoff_smidx" module was used (Regan et al. 2018). In PRMS, recharge is the current available water in the soil recharge zone; details for calculating recharge are presented in Leavesley et al. (1983).

Model Scenarios
The surface-water model runs consist of 32 scenarios for the period from 2018 to 2099 (Table 1). The scenarios used were 30 downscaled LOCA GCMs (LLNL 2021), one forward run using climate data from 1980 through 2016 for the TRB (TRB-fwd), and one trend run applying historical long-term trends from Harwell et al. (2020) to TRB-fwd climate data (TRB-trend). Climate data used in the scenarios were daily precipitation, daily minimum air temperature, and daily maximum air temperature calculated by HRU.
The TRB-fwd scenario was generated by repeating 1980-2016 Daymet data from Thornton et al. (2016) starting in 2017 until the year 2099. Because precipitation progressively increases from the northern extent of the TRB to the southern extent, the TRB-trend scenario was generated using TRB-fwd precipitation and then applying statistically significant historical monthly precipitation trends (Table 1) from Harwell et al. (2020) by section ( Figure 2) from 2018 to 2099. Additionally, temperature inputs for the TRB-trend scenario used a statistically significant upward trend in historical basin air temperature of 0.02 °F per year from Harwell et al. (2020) applied to TRB-fwd minimum and maximum air temperatures.

Section 1 Section 2 Section 3 Section 4 Section 5
January Table 1. Statistically significant historical monthly precipitation trends (in inches/year) for the Trinity River Basin-trend scenario from Harwell et. al. (2020). --indicates no value was reported (months with no reported value did not have a statistically significant trend). Surface Water and Groundwater Availability in the Trinity River Basin in Texas

Numerical Groundwater Model
The numerical groundwater-flow model (McDowell et al. 2023) uses MODFLOW-NWT (Niswonger et al. 2011) to simulate steady-state and transient groundwater flow, recharge, and discharge across the TRB. The full model grid consists of 909 rows, 788 columns, and two layers for a total of 716,292 cells per layer, of which 54,551 are active in each layer, with a total area of approximately 3,369,961 acres. Model grid-cell dimensions, which were selected to align with the 1,000-meter USGS National Hydrogeologic Grid (Clark et al. 2018), are 500 meters in both directions, which equates to slightly less than 62 acres per model cell. The historical groundwater simulation period starts in 1949 (designated as a steady state simulation to obtain initial conditions) and continues through 2018. The future scenarios start in 2018 (designated as a steady state simulation to obtain initial conditions) and continue through 2087.
The surficial extent of the model was set to the surficial exposure extent of the alluvium and terrace deposits containing the TRAA (USGS 2014). Unique model cell values were selected for Quaternary alluvium and terrace geologic units underlying the Trinity River and its tributaries; drain cells (model cells where water leaves the model to the Trinity River and its tributaries); each major lake in the TRB; and all major and minor aquifers underlying the alluvium and terrace deposits containing the TRAA.
The top layer of the model represents the TRAA, and the bottom layer represents underlying geologic units that contain major and minor aquifers (as defined by TWDB [2021b), each with unique hydraulic conductivity values. Grid cells in the bottom layer (all set to a thickness of 250 meters) located in areas between the TWDB-defined extents of the major and minor aquifers-Carrizo, Gulf Coast, Nacatoch, Queen City, Sparta, Trinity, and Woodbine aquifers-were assigned an average hydraulic conductivity value based on the lithology of those areas.
The surface of the model was developed to represent the land surface based on an approximately 30-meter digital elevation model from the USGS National Elevation Dataset (USGS 2020), which was resampled to the 500-meter model cell size. The base of the TRAA was delineated by using Railroad Commission of Texas geophysical logs (Railroad Commission of Texas 2021), TWDB groundwater database geophysical logs (TWDB 2021c), TWDB Brackish Resources Aquifer Characterization System geophysical logs (TWDB 2021d), and TWDB drillers' reports (TWDB 2021e) from wells and boreholes within the boundaries of the Quaternary alluvium and terrace deposits. Keywords used to select drillers' reports containing the alluvium included "alluvium," "alluvial," "sand," "silt," "clay," and "gravel." Reports with these keywords, along with the geophysical logs, were used to make picks on the base of the aquifer from underlying units. The base of wells that were labeled as completed in the alluvium were selected as the aquifer base, as well formations above which sand, gravel, silt, and clay were found. Aquifer thickness values ranged from a minimum of 5 meters (set manually to avoid convergence issues with thin cells) to a maximum of slightly less than 44 meters, based on values provided by the Brazos River alluvium aquifer conceptual model (Ewing et al. 2016). The thickness of the Brazos River alluvium aquifer is likely similar to that of the TRAA because their depositional histories are similar and the formations that contain the aquifers are of similar age, size, and lithology. Because there were little available data to characterize the thickness of the TRAA throughout its extent, the uncertainty associated with the assigned thickness was large, and it is likely that the alluvium is thinner than depicted in parts of the model.
MODFLOW-NWT packages used in the TRAA numerical groundwater-flow model include discretization, basic, upstream weighting, drain, general-head boundary, well, head-observation, output control, and recharge. Detailed descriptions of these packages are presented in the MODFLOW-NWT documentation (Niswonger et al. 2011). All these packages, except for recharge and well, are held constant across all simulated scenarios. The discretization package is used to specify model settings such as layers, rows, columns, cell sizes, and time discretization. The basic package is also used to specify model settings, which include defining active and inactive cells and the starting hydraulic heads for all model cells. The upstream weighting package is used to define storage properties, such as specific storage and specific yield, and flow properties, such as horizontal and vertical hydraulic conductivity. The well package is used to define locations, volumes, and times for groundwater pumpage, or withdrawals. The drain package is used to simulate groundwater flowing out of the aquifer as outflows contributing to surface water. The drain package only accommodates surface-water outflow; it does not allow water to return to the aquifer. In contrast to the drain package, the river package accommodates flow both into and out of the aquifer. Streamflow gain-loss data are commonly used for assessing flows into and out of the aquifer. Because streamflow gain-loss data that could be used to assess groundwater-surface water exchanges between the alluvium aquifer and the Trinity River were not available, the drain package was used instead of the river package. The general-head boundary package is also used to simulate head-dependent flux boundaries. In this case, it was used to handle low-elevation areas with convergence issues near the seaward extent of the model north of Trinity Bay. The head-observation package is used to specify observations of hydraulic head so observed groundwater levels in wells can be compared with the simulated values. Finally, the recharge package is used to input a specified flux distributed across the top layer of the model, in units of meters per day (m/day).

Model Values
To simplify the model, uniform values were used for both specific storage (the volume of groundwater released from one unit volume of the aquifer under one unit decline in hydraulic head) and specific yield (the ratio of the volume of water that a saturated aquifer will yield by gravity relative to the total volume of the aquifer [Johnson 1967]) for both layers. Similarly, to model-thickness values, the specific storage and specific yield values used in the Brazos River Alluvium aquifer conceptual model (Ewing et al. 2016) were also used to guide decisions for these model values. Specific storage was set to 0.0001 inverse meters for the top layer and 1.0x10-7 inverse meters for the bottom layer, whereas specific yield was set to 0.15 for the top layer and 0.01 for the bottom layer. Hydraulic conductivity, a coefficient describing the capacity of a rock to transmit water (Fetter 1994), was set in units of meters per day. The value used for much of the top layer (including Quaternary alluvium and Quaternary terrace geologic units) was 100 m/day, whereas the bottom layer varied 8.4-40 m/day, with a value of 27.4 m/ day for the non-aquifer areas. The value of horizontal hydraulic conductivity for all lakes in the top layer of the model was 1,500 m/day-a high value to allow water to pass through the lake cells because of the complex nature of inter-basin transfer and regulation in the TRB. Vertical hydraulic conductivity for both layers was set to 1% of horizontal hydraulic conductivity, ranging 1-15 m/day for the top layer and 0.08-0.4 m/day for the bottom layer.

Primary Model Control: Recharge
The main control on water going into the groundwater model is recharge, which affects the output to surface water from aquifer storage. Recharge was calculated using the USGS Soil-Water-Balance (SWB) code (Westenbroek et al. 2010) for the historical model recharge calculation and SWB code version 2.0 (Westenbroek et al. 2018) for all climate scenario recharge calculations. SWB 2.0 was used for the climate scenarios because it has been refactored to allow use of Network Common Data Form (NetCDF) version 4 input files, which is the native format of the climate data sets used (LLNL 2021).
SWB uses a modified Thornthwaite and Mather (1957) soil-water-balance method on a gridded data structure to compute the daily volume of net infiltration; net infiltration is assumed to take place any time the soil-moisture value exceeds the total available water for the cell. Inputs for SWB are daily climate data (precipitation, minimum air temperature, and maximum air temperature), elevation, flow direction (generated from the elevation grid), land cover, and soil type. Climate data for the historical model run were acquired from the National Oceanic and Atmospheric Administration's historical daily dataset (NOAA 2021), and climate data for the future scenarios was acquired from Lawrence Livermore National Laboratory (LLNL;LLNL 2021). Land-cover types from the National Land Cover Database (MRLC 2016) were used to assign runoff curve numbers and plant root-zone depths, values that respectively control the surface-water runoff and rate of infiltration through the soil (Westenbroek et al. 2010). Four Gridded Soil Survey Geographic Database (gSSURGO) hydrologic soil groups, categorized from A (high infiltration, low overland flow) to D (low infiltration, high overland flow; NRCS 2021), were used to characterize cells by available water content. All inputs were resampled to fit the 500-meter model grid cells. SWB results were filtered to remove the fifth percentile outliers from each annual result for all scenarios and otherwise left as-is, without an automated calibration process.
Annual recharge values across all scenarios range from about 8,000 acre-feet/year (0.03 inches/year) to about 2,400,000 acre-feet/year (8.55 inches/year) with a mean value of about 620,000 acre-feet/year (2.21 inches/year).

Uncertainty and Sensitivity
Uncertainty in groundwater modeling is assessed in a variety of ways. In this study, hydraulic head values measured in wells were compared to simulated hydraulic head values at the same time (for the historical period) and location in the model. Both observed and simulated values of hydraulic head were plotted and compared to a 1:1 line that represents a perfect fit ( Figure 3) to evaluate the overall simulation-to-observation performance. The figure shows that the simulated values tend to be less than the observed values and, on average, the simulated hydraulic heads are 5.01 meters (6.3%) lower than the observed hydraulic heads (which average just under 80 meters). Over 80% of all simulated hydraulic heads were between 10 meters lower and 5 meters higher than the observed hydraulic heads. Additionally, the differences between the GCMs and the historical model run were compared to understand variance and the range of the simulated values. In contrast to a traditional calibration approach, the outputs of both SWB and MODFLOW were used to bound the uncertainty through the modeling process.
Varying specific yield ratio in the top model layer was also used to assess model sensitivity-understanding how this parameter affects water-budget components is important in understanding model sensitivity. Specific yield values were adjusted from low (0.1) to high (0.2) and compared to what was used in the base model (0.15). When adjusted to low specific yield, the average annual volumetric rate of water leaving to drains of the GCM ensemble mean(the annual mean of the 30 GCMs) increased from about 596,000 acre-feet to about 598,000 acre-feet, an increase of about 0.3%. In contrast, the average annual volume of water going to storage in the aquifer decreased from about 7,000 acre-feet to about 5,000 acre-feet, Surface Water and Groundwater Availability in the Trinity River Basin in Texas a decrease of about 29%. The change of simulated hydraulic heads when adjusted to low specific yield was a negligible increase (less than a hundredth of a percent). When adjusted to high specific yield, the average annual volumetric rate of water leaving to drains of the GCM ensemble mean also decreased by about 2,000 acre-feet/year, whereas the average annual volumetric rate of water going to storage in the aquifer increased by about 2,000 acre-feet/year. The change of simulated hydraulic heads when adjusted to high specific yield was a negligible decrease (less than a hundredth of a percent).

Model Scenarios
The groundwater model runs consist of 32 scenarios for the period from 2018 to 2087 (Table 2) in addition to the historical base run. The scenarios are 30 downscaled LOCA GCMs (LLNL 2021), one forward run using climate data from 1949 through 2017 (TRAA-fwd), and one trend run applying historical long-term trends from Harwell et al. (2020) to TRAAfwd climate data (TRAA-trend). The TRAA-fwd scenario was generated by repeating historical climate data (NOAA 2021) starting in 2018 until the year 2087. To develop the TRAAtrend scenario, local historical climate data (NOAA 2021) were averaged over the period of record from 1949 to 2017 and extrapolated out for future data use starting in 2018 after applying an upward linear rate of change of 0.06 inches/year from Harwell et al. (2020). Additionally, historical air temperature trends of 0.02 °F per year from Harwell et al. (2020) for both daily minimum air temperature and daily maximum air temperature were applied to account for projecting rates of change in the future.

Global Climate Model Scenarios
Thirty LOCA GCM scenarios (LLNL 2021) were selected to evaluate possible future basin conditions as climate inputs to the surface-water and groundwater models developed for this study. Fifteen unique LOCA GCMs were used, each with representative concentration pathways (RCPs) of 4.5 and 8.5 (Table 2). An RCP of 4.5 simulates a scenario that represents moderate global emissions, whereas an RCP of 8.5 simulates a scenario that represents high global emissions. Where the temperature data for each RCP differed, once factored into the models, the variance was minimal because precipitation was the primary driver of model outputs. The USGS Geo Data Portal (Blodgett et al. 2011) was the source of the surface-water model GCM data, and the LLNL (LLNL 2021) was the source for the groundwater model GCM data. Climate data from   List of Trinity River Basin (TRB) model scenarios used for both the surface-water model and the Trinity River alluvium aquifer (TRAA) groundwater model. A representative concentration pathway (RCP) of 4.5 simulates a scenario that represents moderate global emissions, whereas an RCP of 8.5 simulates a scenario that represents high global emissions. The climate data for all scenarios with an RCP value were acquired from the Lawrence Livermore National Laboratory's downscaled climate projections dataset (LLNL 2021). Surface Water and Groundwater Availability in the Trinity River Basin in Texas the GCM scenarios that were used as inputs in both the surface-water model and the TRAA groundwater model include daily precipitation, minimum air temperature, and maximum air temperature. A majority of the GCM data extend to the year 2099; therefore the choice was made to run the surface-water model scenarios to the year 2099. The simulation period for the TRAA groundwater model extends to 2087 to allow for roughly equivalent time periods for both historical data and future projections.
To minimize the number of GCMs in this study and use the most suitable climate data, GCMs were selected based on Venkataraman et al. (2016), which evaluated how well GCM historical data fit the different Texas climate divisions. Metrics used to evaluate the GCMs include mean absolute error, normalized standard deviation, and Kendall's tau. Additionally, because of the variability of the GCM climate data, the model outputs in this study are reported using the annual mean of the 30 GCMs, referred to as GCM ensemble mean. The GCM ensemble mean is used because it reduces the dispersion of the results as compared to individual GCM scenarios.

RESULTS
Historical long-term trends, surface-water model results, and groundwater model results were analyzed using Kendall's tau, Theil slopes, and Mann-Whitney rank-sum tests. When possible, comparisons were made with projected water availability estimates in the 2016 RWPG water plan Region H Water Planning Group et al. 2015).

Historical Long-term Trend Analysis
From the results of previous work by Harwell et al. (2020), precipitation trend analyses on an annual time step in the TRB indicated upward trends in most sections ( Figure 2). Data from eight of the 36 stations selected in Harwell et al. (2020, p. 4) and analyzed for annual streamflow trends indicated upward trends, and all eight stations are in the upper sections (sections 4 and 5) of the study area. None of the data from stations in the lower sections indicated trends in annual streamflow. Data from 16 of the 36 stations indicated upward trends in annual minimum streamflow. All the trends in annual peak streamflow were in the sections that include the Dallas-Fort Worth metropolitan area. Data from two monitoring stations-one USGS streamgage and one U.S. Army Corps of Engineers (USACE) simulated-inflow station-indicated upward trends in annual peak streamflow, and data from one other USGS streamgage indicated a downward trend in annual peak streamflow. Simulated-inflow stations are reservoir stations with inflow data calculated as a mass balance over a 24-hour period. These data were provided by USACE for analyses in Harwell et al. (2020).
Refer to Harwell et al. (2020) for the mass balance equation used by USACE to simulate reservoir inflow.
Of the different river basins included in Harwell et al. (2020), the TRB has the second largest potential flood storage volume at 8,947,349 acre-feet available in the numerous reservoirs built between 1890 and 2013. Potential flood storage volume is defined as the difference between maximum storage volume and normal storage volume. A positive association between potential flood storage volume and annual streamflow was detected at 11 monitoring stations in the TRB, indicating that annual streamflow increases as potential flood storage increases. Data from seven of the 11 monitoring stations also indicated upward trends in annual streamflow. The ratio of streamflow volume to precipitation volume (percent of total water that falls on a watershed that results in streamflow) from analysis in Harwell et al. (2020) was used to estimate the amount of runoff volume to reservoirs in response to upward trends in precipitation in the historical record. Precipitation data in Harwell et al. (2020) included the period from 1900 through 2017. Streamflow volume data included variable time periods ranging from 1869 through 2017, and periods of record for each station are included in Harwell et al. (2020). The purpose of this analysis was to determine how much additional surface water might be available in reservoirs within Region C ( Figure  1) in the future if the upward trends in precipitation reported in Harwell et al. (2020) were to continue. For this study, the Theil slope estimate was calculated for all years and seasons in sections 3, 4, and 5 with statistically significant upward trends in precipitation. The mean annual slope was 0.06 inches/year, and the mean seasonal slope was 0.02 inches/year for each of the three seasons. Figure 4 shows the annual increase in volume by season to the 14 reservoirs analyzed in Harwell et al. (2020) within Region C using the aforementioned ratios of streamflow volume to precipitation volume. Harwell et al. (2020) defined three seasons for the purpose of analysis: season 1 (November, December, January, and February), season 2 (March, April, May, and June), and season 3 (July, August, September, and October). All 14 reservoirs are expected to increase in volume during season 1, with a total annual increase in volume of 3,440 acre-feet/year (Figure 4). Five of the 14 reservoirs are expected to increase in volume during season 2, with a total annual increase in volume of 1,338 acre-feet/year. Lastly, three of the 14 reservoirs are expected to increase in volume during season 3, with a total annual increase in volume of 708 acrefeet/year. Therefore, the estimated total annual increase in volume to the 14 reservoirs from upward trends in precipitation is 5,486 acre-feet, or 0.12% of the projected water availability in 2070 of 4,444,916 acre-feet/year in regions C and H (Table  3). Projected water availability in 2070 is calculated from the values in Table 3 Table 3. Current and estimated future water supply availability (acre-feet/year) as reported in the 2016 regional water planning group water plans Region H Water Planning Group et al. 2015).

PRMS Surface-water Model
The surface-water model results for the scenarios include water budget variables of precipitation as an input and actual evapotranspiration, Hortonian runoff, and recharge as outputs. Figure 5 shows annual means of the water budget variables from 2018 to 2099 for all the model scenarios: the 30 GCMs, TRB-fwd, TRB-trend, and the GCM ensemble mean, which is the annual mean of all the GCM scenarios. Precipitation is a large driver for model outputs, as seen in the Figure  5. As precipitation increases, Hortonian runoff and recharge values increase.
From 2018 to 2099, TRB-fwd and TRB-trend mean annual precipitation values were 42.11 and 43.74 inches/year, respectively, a difference of 1.63 inches/year (Table 4)       a Kendall's tau value of -0.2314 (p-value ≤ 0.01), indicating a downward trend in actual evapotranspiration. The computed Theil slope estimate for the period from 2018 to 2099 is -0.0127, which equates to a downward trend in actual evapotranspiration of about 0.01 inches/year. This is likely a result of less water input to the system as precipitation, resulting in less surface water available for evapotranspiration.
From 2018 to 2099, TRB-fwd and TRB-trend mean Hortonian runoff values were 1.86 inches/year and 1.94 inches/year, respectively, a difference of 0.08 inches/year. Neither the TRBfwd nor the TRB-trend mean annual Hortonian runoff values indicate statistically significant time trends. The GCM ensemble mean had a mean annual Hortonian runoff value of 1.53 inches/year. The time trend of the GCM ensemble mean data is statistically significant with a Kendall's tau value of -0.2096 (p-value ≤ 0.01), indicating a downward trend in Hortonian runoff. The computed Theil slope estimate for the period from 2018 to 2099 is -0.0013, which equates to a downward trend in Hortonian runoff of about 0.0013 inches/year. TRB-fwd and TRB-trend mean annual recharge values were 5.03 inches/year and 5.2 inches/year, respectively, a difference of 0.17 inches/year. Neither the TRB-fwd nor the TRB-trend mean annual recharge values indicate statistically significant time trends. The GCM ensemble mean for 2018 to 2099 yielded a mean annual recharge value of 3.64 inches/year. The time trend obtained from the surface-water model for the GCM ensemble mean is statistically significant with a Kendall's tau value of −0.4263 (p-value ≤ 0.01), indicating a downward trend in recharge. The computed Theil slope estimate for the period from 2018 to 2099 is −0.0092, so the downward trend in recharge is about 0.0092 inches/year. Decadal plots of the GCM ensemble data are included to show the distributions and the variability in the GCM data by decade (Figures 6-9). The violin plots (Hintze and Nelson 1998) were constructed to visualize the distribution and probability density of the GCM ensemble data (Figures 6-9) and include 300 data points per decade, representing 10 years of annual variable data for each of the 30 GCM scenarios. Variables plotted include precipitation (Figure 6), actual evapotranspiration (Figure 7), Hortonian runoff (Figure 8), and recharge ( Figure 9). According to results of the Mann-Whitney rank-sum test (p-value ≤ 0.05), the 2090 decadal precipitation data are not equivalent to the 2020, 2030, and 2040 decadal precipitation data. As expected from the downward precipitation trend of the GCM mean annual precipitation data previously discussed, the 2090 decadal precipitation data are significantly less than the 2020, 2030, and 2040 decadal precipitation data ( Figure 6). The differences between the 2090 decadal precipitation data and the 2050, 2060, 2070, and 2080 decadal precipitation data are not statistically significant.
Decadal differences for evapotranspiration, Hortonian runoff, and recharge are similar to those of precipitation.
Mann-Whitney test results indicate that the 2090 decadal values are significantly less than the earlier decades of 2020, 2030, and 2040 for the other components of the water budget, except for evapotranspiration. For evapotranspiration, there is no difference between the 2090 and 2030 data. However, 2090 evapotranspiration values are less than 2020 and 2040 values. These results are expected given the downward trends in the ensemble means for all water budget components.

Numerical Groundwater Model Inflows, Outflows, and Storage
The simulated output of the numerical groundwater model for the TRAA-forward, TRAA-trend, and GCM scenario runs includes volumetric water budget components in acre-feet/year for the inputs and outputs of the model (Figure 10). These water budget components are recharge, drains, storage change, cumulative storage change, head-dependent boundaries, and wells. Except for groundwater storage change, positive values indicate groundwater going into the model (groundwater inflows), whereas negative values indicate groundwater flowing out of the model (groundwater outflows).
Recharge to the aquifer-the primary inflow to the MOD-FLOW model-is largely controlled by precipitation but is also dependent on other variables, such as land cover and soil type. Groundwater storage change represents the change in groundwater being stored in the aquifer at any given time. When groundwater storage change is negative, net groundwater is entering the aquifer, and when groundwater storage change is positive, net groundwater is leaving the aquifer. Cumulative groundwater storage change is simply the cumulative change in aquifer storage by year, representing the change in simulated volume of groundwater stored in the aquifer in any given year. For reference, at the end of the base model transient period, the model estimates there are about 4.3 million acre-feet in storage for the TRAA. Drains signify the volume of water flowing out of the aquifer and include all simulated reaches of rivers and streams. Head-dependent boundaries show groundwater flowing out of the aquifer (and sometimes entering the system) at the seaward extent of the model near Trinity Bay, and wells show groundwater being pumped out of the system.
Because of the many different approaches and objectives when creating the GCMs (Table 2), the GCM input data exhibit a high level of variance, which is reflected in the output of the various GCM scenarios and is why the GCM ensemble mean is analyzed. Annual simulated values for recharge from all GCM scenarios range from about 8,000 acre-feet to about 2,385,000 acre-feet, whereas the GCM ensemble mean ranges from about 377,000 acre-feet/year to about 743,000 acre-feet/ year (    The GCM ensemble mean, TRAA-fwd run, and TRAAtrend run water budget components are compared in Figure  10. Because of the linear nature of the trend extrapolation for the trend run (which was based on the 1900-2017 period of record), it is likely a high estimate of the range of possible outcomes. Conversely, the GCM ensemble mean indicates that the overall water availability for the future will likely be within the upper and lower bounds of water availability determined in the historical simulation period .
Decadal analysis of the water budget also provides insight into future groundwater conditions in the TRB. Violin plots of the primary water budget components (recharge, drains, and storage change) were created to better understand the distribution of values across decades (Figures 11-13). Figure 11 shows the recharge component of the groundwater model water budget by decade, including every year in each decade for all 30 GCM scenarios. Because recharge is the main control on volume of groundwater in the model, more recharge to the model means more groundwater is available to leave the model via drains or to be stored in the aquifer. Simulated median annual recharge (Figure 11) ranges from about 632,000 acre-feet/ year in the 2020 decade to about 462,000 acre-feet/year in the 2070 decade. Simulated median annual drain values ( Figure  12) range from about 633,000 acre-feet/year out of the system in the 2040 decade to about 517,000 acre-feet/year out of the system in the 2070 decade. Simulated median annual storage change (Figure 13), while highly variable, ranges from about 16,000 acre-feet/year being added to aquifer storage in the 2020 decade to about 40,000 acre-feet/year flowing out of the aquifer in the 2050 decade.
Negative cumulative storage change values mean that water is going into the aquifer-an increase in groundwater storage in the alluvial aquifer. There is a high degree of variability across different climate scenarios, as some scenarios show groundwater will be increasing in the aquifer whereas others show groundwater will be decreasing. Typically, when there is more recharge to the aquifer for a given scenario and year, cumulative storage change values will be more negative, representing an increase in groundwater storage and vice versa. Figure 14 shows that the GCM ensemble mean has an initial increase of groundwater available (negative values) in the alluvial aquifer, followed by a gradual decrease over the rest of the projected time period.

Limitations
All water budgets, including the results of this study, have uncertainties associated with them because of simplifying assumptions within the models. Uncertainties for the results of this study stem from measurement and modeling errors, as well as natural variability in precipitation patterns, evapotranspiration, soil and vegetation properties, and diurnal, seasonal, and long-term climate trends (Healy et al. 2007). Model limitations associated with PRMS include the following four factors: (1) groundwater or reservoir withdrawals are not simulated, (2) interbasin transfers are not simulated, (3) land use changes are not simulated, and (4) the effect of frozen ground on runoff is not simulated (Bjerklie et al. 2015). Additional uncertainties exist specifically for the TRAA. Because the TRAA is not formally classified as an aquifer, few data have been collected to characterize its geologic or hydraulic properties, and no gainloss studies have been done to assess groundwater-surface water exchanges between the alluvium aquifer and the Trinity River. As such, this groundwater model should be used for basin-wide water budget assessments and not localized regions within the TRB. Population growth, land-use changes, and other likely anthropogenic effects that could appreciably reduce water availability in the future were not considered in the simulations described herein.
Application and interpretation of outputs from GCM simulations require an understanding of some basic considerations and limitations of the results as summarized by Taylor et al. (2012). These include unforced variability, bias correction, downscaling, and multi-model ensemble (Krinner et al. 2020;Soriano et at. 2019). GCMs typically have grid cells that are approximately 100 kilometers by 100 kilometers. Climate change models are being tasked to provide climate change effects on increasingly smaller spatial scales. To accomplish this task, downscaling techniques have been developed to take GCM output and provide meaningful information at scales smaller than the size of the GCM's grid cells (Taylor et al. 2012). LOCA is a statistical downscaling technique that uses history to add improved fine-scale detail to GCMs (Pierce et al. 2014). LOCA errors tend to pattern those of random variability in sampling as opposed to errors showing spatial pattern biases. However, like all statistical models, LOCA is based on historical data and thus assumes that spatial relationships and local and average climate fields will not change for future climate scenarios.
Research has shown that the use of multi-model ensemble means and the ensemble range (spread of minimum and maximum across many GCMs) will provide better projections of future changes in climate and represent the most conservative approach given the challenges of predicting complicated systems like climate (Venkataraman et al. 2016). Therefore, longterm climate trends from an ensemble mean of many GCMs Surface Water and Groundwater Availability in the Trinity River Basin in Texas  will provide a broad representation of future climate conditions, not an accurate description of the timing and magnitude of individual future events (Venkataraman et al. 2016).

DISCUSSION
The TRB is predominately within two regions of the TWDB RWPGs: Region C is in the upstream part of the TRB, and Region H is in the downstream part of the TRB (Figure 1). Both RWPGs inform stakeholders of the water supply and demand by providing and updating a water plan every 5 years. Although water plans for these two regions are available for 2021 (Freese and Nichols, Inc. et al. 2020;Region H Water Planning Group et al. 2020), for the purpose of this study, values provided in the 2016 water plans were used Region H Water Planning Group et al. 2015). Using land use change projections and predicted population changes in the area, TWDB and the RWPGs produce water supply availability values for the regions for the next 50 years in decadal increments. Local entities use the information provided in the RWPG water plans to make management decisions regarding water use. Although the RWPGs' water plans are thorough and comprehensive, one limitation is that they do not consider future climate change as part of their projected water budget. This study aimed to provide an initial framework for stakeholders to evaluate the vulnerability of the TRB to a changing climate through the use of GCMs. Howev-er, the scope of this study is limited by simplified models that do not account for the complex regulation in the upper part of the TRB or the TRB's projected anthropogenic changes, which include substantial population growth and land-use changes, particularly in and near the Dallas-Fort Worth and Houston metropolitan areas. As such, the results of this study and the findings of the RWPG water plans cannot be compared directly but could be studied in conjunction to better understand TRB water availability.
The mean annual increase in precipitation to Region C reservoirs from historical data (1900-2017) is 0.06 inches/year, and the mean seasonal increase is 0.02 inches/year for each of the three seasons. The estimated total annual volume increase from upward precipitation trends reported in Harwell et al. (2020) to 14 Region C reservoirs is 5,486 acre-feet/year  inches/year, or 0.12% of regions C and H's projected water availability in 2070 of 4,444,916 acre-feet/year (Table 3). Projected water availability in 2070 is calculated from the values in Table 3 as the sum of the total volumes including groundwater in regions C and H and subtracting the sum of the groundwater volumes in both regions.
However, according to the surface-water model analysis, the GCM ensemble mean annual precipitation indicates a downward trend in precipitation of about 0.03 inches/year (2018-2099), resulting in a downward trend in Hortonian runoff. The surface-water model GCM ensemble mean indicates a downward trend in Hortonian runoff of about 0.0013 inches/year. Simulated results of the surface-water model GCM ensemble mean indicate a downward trend in actual evapotranspiration of about 0.01 inches/year, likely a result of less water input to the system as precipitation and therefore less surface water available for evapotranspiration. Lastly, the surface-water model GCM ensemble mean indicates a downward trend in recharge of −0.0092 inches/year or approximately 8,764 acrefeet/year to the TRB.
Similar to the surface-water model results, the TRAA groundwater model analysis of the GCM ensemble mean indicates a downward trend of about 2,376 acre-feet/year of recharge to the TRAA, or 0.30% of the projected 2070 groundwater availability in regions C and H ( Table 3). The downward trend in recharge is a result of the downward trend in precipitation from the GCM ensemble mean. The TRAA groundwater model GCM ensemble mean indicates an upward trend in the amount of groundwater flowing out of the aquifer to rivers and streams at a rate of about 1,877 acre-feet/year. The TRAA model GCM ensemble mean also indicates an upward trend of 3,655 acre-feet/year in cumulative storage change, and the amount of groundwater in the aquifer is decreasing despite an initial increase in groundwater storage, as depicted by the annual GCM ensemble mean data in figure 14.
The GCM ensemble mean and trend scenario results show long-term stability in the water budget for both surface water and groundwater for the study period. The estimated total annual increase in volume to 14 Region C reservoirs from upward trends in precipitation is 5,486 acre-feet/year . However, the surface-water model GCM ensemble mean annual precipitation indicates a downward trend in precipitation of about 0.03 inches/year (2018-2099), resulting in a downward trend in Hortonian runoff of approximately 0.0013 inches/year. The GCM ensemble mean for the surface-water model (11,471,544 acres) indicates a downward trend in recharge of about 8,764 acre-feet/year. The GCM ensemble mean for the TRAA groundwater model (3,369,961 acres) also indicates an upward trend (decrease of groundwater storage) of 3,655 acre-feet/year in cumulative storage change, and the amount of groundwater in the aquifer is decreasing despite an initial increase in groundwater storage. The results of this analysis show that the overall change in future water availability attributable to climate change is small, assuming the average of the ensembles is the best predictor of the future. The scientific consensus is that water availability in the future will be more variable compared to the past because of the likelihood of longer and more severe droughts punctuated by more intense storms (Nielsen-Gammon et al. 2021).
Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.