Numerical groundwater models for Edwards Aquifer systems Open Access

Water-resource managers rely on groundwater flow models to characterize groundwater flow conditions in the Edwards (Balcones Fault Zone) Aquifer and to serve as the basis for predicting effects of water-resource management decisions. There are recognized limitations and shortcomings in all numerical groundwater models, such as uncertainty about the underlying conceptual model and its implementation in the ensuing numerical model (De Marsily, 1986; Anderson and Woessner, 2002).

As described throughout this memoir, the Edwards Aquifer is regional in scale with properties and characteristics that vary spatially and, often, temporally. Capturing these extreme heterogeneities imposes a large cost with regards to computation time and data input requirements. At the same time, aquifer management concerns are usually defined at the local scale, requiring sufficiently high resolution to account for local-scale characteristics of the aquifer. For these reasons and others, no single numerical model has been developed to simulate flow for the entirety of the Edwards Aquifer. This chapter discusses the major numerical models developed for individual segments and subsets of the Edwards Aquifer, including the San Antonio segment, the Barton Springs segment, Kinney County, and the Devils River watershed. In addition, there have been regional-scale models, referred to as groundwater availability models, that span large areas of the Edwards Aquifer and related aquifers, such as the Edwards-Trinity and Trinity Aquifers. Although lumped parameter (Wanakule and Anaya, 1993; Loáiciga et al., 2000) and solute transport (Lindgren et al., 2011) models of the Edwards Aquifer have been developed, this survey focuses on distributed groundwater flow models.

Several regional-scale groundwater flow models cover parts of the Edwards Aquifer. Models developed after the year 2000 have been mostly generated as part of the Texas Water Development Board Groundwater Availability Model program (Mace et al., 2008). These models target portions of the Edwards and Edwards-Trinity Aquifers, including: the northern segment (Jones, 2003), the Barton Springs segment (Scanlon et al., 2001, 2003; Smith and Hunt, 2004; Winterle et al., 2009; Hutchison and Hill, 2011), the San Antonio segment (Klemt et al., 1979; Maclay and Land, 1988; Lindgren et al., 2004; Lindgren, 2006; Fratesi et al., 2015; Liu et al., 2017), Kinney County (Hutchison et al., 2011a), and Val Verde County watersheds (Eco-Kai Environmental, Inc., and Hutchison, 2014; Green et al., 2015, 2016; Toll et al., 2017). In addition, there are groundwater availability models for the Edwards-Trinity Aquifer (Anaya and Jones, 2004, 2009; Young et al., 2009; Hutchison et al., 2011b) and the Hill Country Trinity Aquifer (Jones et al., 2009; Mace et al., 2000), although the focus of these models includes the Trinity Aquifer (see Sharp et al., this volume, their figures 2 and 3). There was also a U.S. Geological Survey Regional System Aquifer Analysis model that targeted the Edwards Plateau, but it covered the western portions of the Edwards Aquifer (Kuniansky and Holligan, 1994; Kuniansky, 1995). Models of interest are those models wherein the scales allow for in-depth analysis of local flow regimes. Models listed in Table 1 are discussed by area in the following sections.

Klemt et al. (1979) developed the first comprehensive groundwater model of the San Antonio segment of the Edwards Aquifer using the U.S. Geological Survey Groundwater Simulation Program GWSIM (Fig. 1), based on the groundwater simulation program developed by Prickett and Lonnquist (1971). Their model domain extended from a groundwater divide located near Brackettville, Texas, in Kinney County on the west to a groundwater divide located near the Blanco River in Hays County on the east (Fig. 1). The northern or updip boundary was the upgradient edge of the Edwards Aquifer recharge zone and the southern or downdip boundary was designated as the downdip edge of freshwater in the Edwards Aquifer, which was defined as less than 1000 mg/L dissolved solids. Thorkildsen and McElhaney (1992) reevaluated the Klemt et al. model using refined water levels and spring flows. The Thorkildsen and McElhaney model domain was the same as that in the Klemt et al. model.

Maclay and Land (1988) developed an independent groundwater flow model that relied on a groundwater flow simulator developed by Trescott et al. (1976). The southern or downdip boundary of their model domain was essentially the same as that designated in the Klemt et al. model (Fig. 1); however, the northern or updip boundary was modified based on a refined interpretation of the recharge zone (Puente, 1978). In particular, the northern limit of the recharge zone was extended to include more area in eastern Kinney County and most of Uvalde County. The rest of the recharge zone in Bexar, Comal, and Hays Counties remained the same as that defined by Klemt et al. A regional-scale, finite-element model that included the Edwards, Trinity, and Edwards-Trinity Aquifers was developed to provide a large-scale assessment of the three aquifers (Kuniansky and Holligan, 1994; Kuniansky, 1995).

The Edwards Aquifer Authority commissioned a groundwater flow model of the Edwards Aquifer that was completed in 2004 (Lindgren et al., 2004). The northern or updip boundary of the model domain was essentially the same as that defined by Maclay and Land (1988). The location of the southern or downdip boundary at the freshwater–saline water interface was revised from 1000 mg/L to 10,000 mg/L, thereby extending the model domain farther south (Fig. 1). MODFLOW-96 and MODFLOW-2000 (Harbaugh and McDonald, 1996; Harbaugh et al., 2000) were both used to simulate flow.

Fundamental similarities among the Klemt et al., Maclay and Land, and Lindgren et al. models were treatment of the northern or updip model boundary and the way in which recharge was incorporated. None of these models included the Edwards Aquifer contributing zone. In all three models, recharge from the contributing zone was input into the upstream boundary of the Edwards Aquifer recharge zone. Recharge rates were based on surface-flow measurements and estimates provided by the U.S. Geological Survey (Puente, 1978).

The 2004 MODFLOW model (Lindgren et al., 2004) incorporated focused, high-transmissivity zones to represent conduit flow, as conceptualized in the Edwards Aquifer by Hovorka et al. (2004) and Worthington (2004). The 2004 version of the MODFLOW model was used for a series of assessments of recharge and recirculation scenarios (Todd Engineers, 2004, 2005, 2008). Lindgren (2006) modified the 2004 MODFLOW model by replacing the high-transmissivity representations of conduits with a diffuse-flow conceptualization. The 2006 version of the MODFLOW model was used for additional simulations of the Edwards Aquifer (Lindgren, 2006). The 2004 and the 2006 versions of the model were similar in both their residual calibration errors and their ability to predict spring-flow and water-level responses to variations in recharge and pumping.

In 2011, the Edwards Aquifer Authority undertook two initiatives to reduce uncertainty in models used to perform waterresource management analyses. One initiative was undertaken to advance the 2004 MODFLOW model through a series of recalibration exercises (Liu et al., 2017). The second initiative was to develop a second groundwater flow model that was conceptually independent of the 2004 model (Fratesi et al., 2015). The objective of a second model was not to replace the 2004 model, but to provide an independent numerical tool against which to compare model predictions. The alternative model (Fratesi et al., 2015) employs a finite-element formulation instead of the finitedifference formulation used in Lindgren et al. (2004) and Liu et al. (2017). FEFLOW (Diersch, 2014) was selected as the finiteelement modeling software package.

The FEFLOW model domain was specified to allow for model boundaries to be no-flow, to the degree possible, with the exception of springs (Fratesi et al., 2015). The FEFLOW model domain included all groundwater and surface-water basins that contribute water to the San Antonio segment of the Edwards Aquifer (Fig. 1). Thus, this domain includes the entire contributing zone of the Edwards Aquifer. The major river basins in the contributing zone were characterized as hydrologically independent. By doing this, surface-water and groundwater flow from each basin to adjoining basins was minimized. This characterization honored the conceptual model developed for the contributing zone in which surface-water and groundwater flow in each basin was mostly restricted to each basin (Fratesi et al., 2015). This conceptualization also allowed the precipitation/recharge model to be calibrated for each watershed.

The fundamental difference between the MODFLOW and FEFLOW models is the manner in which recharge is input into each model. By including recharge calculated directly from precipitation over both the contributing and recharge zones, the time lag between the time of precipitation and the time at which a hydraulic signal is transmitted through the aquifer is captured in the Fratesi et al. model. Conversely, recharge is a specified input parameter in the Lindgren et al. (2004) and Liu et al. (2017) models. Although other factors differ in the FEFLOW and MODFLOW models, it is the manner in which recharge is incorporated that establishes the FEFLOW model as conceptually independent from the MODFLOW model.

Liu et al.’s (2017) updated model of the San Antonio segment of the Edwards Aquifer made use of improved data availability to implement improvements such as the use of actual pumping locations and better estimates of pumping rates, and the addition of significantly more water-level observation well locations for use in model calibration. The model domain for the Liu et al. (2017) model included the recharge and confined zones of the Edwards Aquifer, but not the contributing zone. Structural changes to the model included removal of the Barton Springs segment, addition of barrier features to represent the Knippa Gap area and Haby’s Crossing fault, modifications to hydraulic conductivity zones, and additional drain discharge locations.

These structural changes and subsequent recalibration of the MODFLOW model significantly reduced the errors between modeled and observed water levels in certain areas—particularly in the recharge zone and the Knippa Gap area (Liu et al., 2017). Since the updated model (Liu et al., 2017) and the original model (Lindgren et al., 2004) were calibrated for different time periods, it was difficult to quantify the overall improvement in terms of the ability of the model to match observations. The updated model had less error in water-level and spring-flow predictions compared to the original model for the period of 2001–2009, but the original model by Lindgren et al. (2004) was not calibrated to that period. The updated model was specifically calibrated to that period and therefore should be expected to perform better. When the updated model was run forward for the 4 yr validation period of 2012–2015, the match to observed water levels and spring flows was generally satisfactory, but not quite as close as for the period of 2001–2011 to which it was calibrated.

When the updated model was run for 1947–1958 to simulate the drought-of-record, the match to observations at index well J-17 and Comal Springs was qualitatively good, as was the original 2004 version. The updated model slightly underestimated the lowest flows at San Marcos Springs compared to the original model, but sensitivity analyses showed that the model fit to San Marcos Springs could be significantly improved by a modest adjustment to recharge in that area of the model. The updated model consistently underestimated water levels at index well J-27 in Uvalde County, but the observation data for that well during that time are questionable due to structural flaws in the well casing that resulted in the well straddling several aquifers (Green, 2013).

Overall, the model by Liu et al. (2017) represents an incremental improvement in the iterative modeling process with better representation of pumping and inclusion of new structural features that clearly reduce water-level errors in those areas of the model. A limitation of the ability of the Liu et al. (2017) model (or any model) to match observations is the uncertainty in its inputs. While uncertainty will always exist regarding the distribution of hydrologic properties in the subsurface, the Liu et al. (2017) model does adequately match the magnitude and timing of observed hydrologic responses to transient changes in recharge and pumping rates throughout the modeled region. Uncertainty in estimates of monthly recharge is possibly the most significant contributor to model error (Liu et al., 2017). The FEFLOW model by Fratesi et al. (2015) offers the opportunity to simulate recharge due to precipitation applied to the contributing and recharge zones.

Scanlon et al. (2001, 2003) developed the original groundwater availability model for the Barton Springs segment of the Edwards (Balcones Fault Zone) Aquifer using MODFLOW-96 (Harbaugh and McDonald, 1996; see also Sharp et al., this volume, their figure 2). The model was calibrated for steady-state conditions using a recharge rate based on the average spring discharge of 1.56 m³/s plus pumpage of 0.14 m³/s for 1989. The model was calibrated to transient conditions for the 10 yr period from 1989 to 1998. This model generally reproduced spring-flow observations for the transient period. However, simulated spring flows during high-flow periods tended to be larger than observed, and simulated water-level hydrographs for observation wells fluctuated more strongly in response to changes in recharge than the observed hydrographs. Simulated heads in the southwest portion of the model tended to be too high. Additionally, simulated spring flows for extreme drought conditions tended to be lower than observed. To address this problem, the Barton Springs Edwards Aquifer Conservation District developed an alternative calibration of the model for simulating extreme low-flow conditions (Smith and Hunt, 2004).

A dual-continuum model, MODFLOW-DCM version 2.0 (Painter et al., 2006, 2007), was developed to simulate diffuse-conduit flow in the Barton Springs segment, built on the groundwater availability model by Scanlon et al. (2001, 2003). Winterle et al. (2009) used MODFLOW-DCM to extend the model domain to the south and west. These updates necessitated recalibration for steady-state and transient conditions. The steady-state model was calibrated to match observations of hydraulic heads in 74 wells during July and August 1999 using an estimated average spring discharge of 1.56 m³/s and total pumpage of 0.14 m³/s. The transient MODFLOW-DCM model was calibrated to match spring discharge measurements at Barton Springs and hydraulic head responses in three observation wells for the 10 yr period from January 1989 to December 1998. In addition to the model updates and calibrations, a simple algorithm was developed for estimating recharge input to the model based on average monthly precipitation over the recharge and contributing zones.

The root-mean-square residual error for the MODFLOW-DCM (Winterle et al., 2009) steady-state model was reduced to approximately half that of the original groundwater availability model (Scanlon et al., 2001, 2003). Using recharge input developed from the precipitation-to-recharge algorithm, the transient model was able to reproduce the pattern of hydraulic head fluctuations in observation wells, as well as the pattern of observed discharge at Barton Springs. The calibrated transient model performed well in matching the spring-flow conditions and the cumulative spring flow for the 10 yr simulation period.

Hutchison and Hill (2011) recalibrated the groundwater availability model for the Barton Springs segment of the Edwards Aquifer that was built by Scanlon et al. (2001, 2003). The existing MODFLOW-96 (Harbaugh and McDonald, 1996) packages used by Scanlon et al. (2001, 2003) were converted to MODFLOW-2000 (Harbaugh et al., 2000). MODFLOW-2000 was used with the Geometric Multigrid (GMG) solver (Wilson and Naff, 2004). The resulting alternative groundwater availability model was developed to ascertain the amount of pumping that would result in specified spring flows of 0.31, 0.25, 0.20, 0.14, and 0.08 m³/s under drought-of-record conditions.

The 2001 groundwater availability model for the Barton Springs segment of the Edwards (Balcones Fault Zone) Aquifer had been calibrated based on data from 1989 to 1998 (Scanlon et al., 2001, 2003). Thus, the calibration period did not include the historic drought-of-record that lasted from 1950 through 1956, when the minimum observed discharge of 0.31 m³/s (Slade et al., 1986) occurred for Barton Springs. To account for the drought-of-record, Hutchison and Hill (2011) recalibrated the model for the Barton Springs segment of the Edwards (Balcones Fault Zone) Aquifer for the period January 1943 to December 2004 using 744 estimated or measured discharges for Barton Springs (Slade et al., 1986). Simulated discharges at Barton Springs using the alternative groundwater availability model provided satisfactory agreement with the minimum observed discharge of 0.31 m³/s that occurred in July and August of 1956 during the historic drought-of-record.

Although previous models (i.e., Klemt et al., 1979; Maclay and Land, 1988; Lindgren et al., 2004; Anaya and Jones, 2004, 2009; Young et al., 2009; Fratesi et al., 2015; Liu et al., 2017) covered at least part of the Edwards Aquifer in Kinney County, only the model by Hutchison et al. (2011a) was developed exclusively for the Kinney County segment of the Edwards Aquifer (Fig. 2). The Hutchison et al. (2011a) model, developed using MODFLOW-2000 (Harbaugh et al., 2000), covered Kinney County and its surrounding counties with a refined hydrogeologic framework. This model was developed by the Texas Water Development Board and is considered an alternative groundwater availability model.

The Hutchison et al. (2011a) model was developed to evaluate the groundwater resources of the Kinney County segment of the Edwards Aquifer and to simulate the effects of potential groundwater withdrawal at wells on springs and river flows. Specifically, the objective of the Kinney County model was to ascertain how pumping levels impacted flows at the three largest springs (Las Moras, Mud, and Pinto Springs) under potential future groundwater-use conditions. The model contained one steady-state stress period (stress period 1) that was intended to initiate a transient simulation, and 56 transient annual stress periods. The transient stress periods for the years 1950 through 2005 included the drought-of-record. Model boundaries were reviewed against observed data and modified, as necessary. A FORTRAN preprocessor was developed to estimate net groundwater recharge distribution based on a revised method originally developed by Maxey and Eakin (1949). Groundwater pumping records that were assigned to Kinney County in the alternative groundwater availability model (Hutchison et al., 2011a) were significantly greater than historical records used for Kinney County in the regional Edwards-Trinity Aquifer groundwater availability model (Anaya and Jones, 2004, 2009; Young et al., 2009). Models calibrated to pumping levels much greater than actual pumping rates can lead to estimation of elevated hydraulic property values or improper boundary conditions.

A groundwater flow model was developed in 2014 for the Edwards-Trinity Aquifer within Val Verde County (Eco-Kai Environmental, Inc., and Hutchison, 2014). The objective of the model was to determine correlations and potential effects of groundwater pumping on local spring flows, lake elevations, and groundwater levels in Val Verde County. The Kinney County groundwater model (Hutchison et al., 2011a) also encompassed Val Verde County and was used as the foundation for the Val Verde County model. In particular, the geologic framework and many of the boundary conditions of the Kinney County model were used as the foundation of the Val Verde County model (Eco-Kai Environmental, Inc., and Hutchison, 2014). Although both models were constructed using half-mile elements, the Kinney County model was developed using annual stress periods, while the Val Verde County model was developed using monthly stress periods from 1968 to 2013.

Two of the watersheds that discharge to the Rio Grande in Val Verde County, i.e., the lower Pecos River watershed (Green et al., 2016) and the Devils River watershed (Green et al., 2015; Toll et al., 2017), have been recently modeled (Fig. 3). Parsing out watersheds within the Edwards-Trinity Aquifer was justified by the hydrologic separation exhibited by watersheds in the Edwards Plateau, particularly where the aquifer is phreatic. Similar hydrologic separation of phreatic aquifers was observed in the contributing zone of the Edwards Aquifer (Fratesi et al., 2015).

Fast conduit flow and slow diffuse flow in the karstic Edwards-Trinity Aquifer were replicated in the lower Pecos River (Green et al., 2015) and the Devils River (Green et al., 2015; Toll et al., 2017) watershed models using highly transmissive preferential flow paths aligned with the river channel pathways (Green et al., 2014). The Devils River model was eventually converted to a coupled surface-water/groundwater model to accommodate recharge (Toll et al., 2017). The coupled surface-water/groundwater model successfully replicated both flashy flow and low base flow in the Devils River at major springs and where it discharges to Amistad Reservoir.

Regional and local groundwater models have been developed to help manage the resources of the Edwards Aquifer. Given the interdependent hydrologic connections between the Trinity and Edwards Aquifers, integration of the Trinity and Edwards Aquifers within the same model is commonly undertaken and frequently necessary. This allows for boundary conditions that are more easily designated and justified. Regional-scale groundwater models, such as the Texas Water Development Board–designated groundwater availability models, have limited use in simulating local-scale pumping/drawdown scenarios; however, these models have been useful in the determination of the modeled available groundwater based on desired future conditions (Mace et al., 2008).

Simulating local-scale pumping/drawdown scenarios in the Edwards Aquifer and reducing the limitations and often unacceptable levels of uncertainty of existing models will likely require smaller-scale models or regional-scale models that have sufficiently high resolution to account for local-scale characteristics of the aquifer. As discussed in this chapter, improved simulation of the Edwards Aquifer has been achieved as regional-scale groundwater flow models are improved and local-scale models are developed. Additional future advances in numerical groundwater models of the Edwards Aquifer and related systems will benefit from improved and expanded input and property data, refined conceptual models of the targeted aquifer systems, and increased temporal and spatial resolution in the numerical renderings of the physical systems.

Limitations in existing models remain. Advances are needed to resolve or reduce uncertainties to improve the capability of groundwater models used to manage the resources of the Edwards Aquifer. One area of uncertainty is the southern or downdip boundary of the San Antonio segment of the Edwards Aquifer. The Klemt et al. (1979) and Maclay and Land (1988) models designated the 1000 ppm southern boundary as a no-flow boundary. Lindgren et al. (2004) and Fratesi et al. (2015) still modeled the southern boundary as no-flow, but at the 10,000 ppm isopleth. The 1000–10,000 ppm brackish-water zone downdip of the so-called bad-water line might be an area of horizontally convergent flow within the Edwards Group (Hoff and Dutton, 2017). There are two issues here: whether a no-flow boundary is correct, and what error propagates from ignoring variable-density flow updip of the 10,000 ppm boundary (Brakefield et al., 2015). If the conclusion made by Hoff and Dutton (2017) is correct, convergent flow implies significant vertical flow. Existing two-dimensional Edwards Aquifer models, however, assume no vertical or cross-formational flow, or they allow or imply upward-directed cross-formational flow from the Trinity Aquifer. Hoff and Dutton (2017) speculated that even a small percentage error in a model’s volume balance, owing to cross-formational flow not being represented, might have an impact on calibration of recharge and of aquifer properties.

Characterization of interformational flow between the Trinity and Edwards Aquifers is a second critical area of uncertainty. Early estimates of Trinity-Edwards Aquifer interformational flow of 66.4 × 10⁶ M³/yr (53,800 acre-ft/yr) (Lowry, 1955) and 132.0 × 10⁶ M³/yr (107,000 acre-ft/yr) (Bader et al., 1993) included only the Cibolo Creek watershed. Interformational flow from the Trinity Aquifer to the Edwards Aquifer was not included in the model by Klemt et al. (1979). Subsequent models by Maclay and Land (1988) and Lindgren et al. (2004) did include inflow from the Trinity Aquifer as a source of groundwater. The domain of the model by Kuniansky and Holligan (1994) and Kuniansky (1995) incorporated the Edwards-Trinity, Trinity, and Edwards Aquifers, and so interflow was inherently included in the model. Maclay (1995) identified two areas of groundwater inflow along the updip limit of the San Antonio segment of the unconfined Edwards Aquifer, one area in northeastern Medina County and the other in Comal County (Maclay and Land, 1988). The Maclay and Land (1988) model did not indicate significant inflow from the Trinity Aquifer to the Edwards Aquifer in either Kinney or Uvalde Counties.

Lindgren et al. (2004) calculated that inflow through the northern and northwestern model boundaries contributes 6.5% of total recharge to the Edwards Aquifer. Mace et al. (2000) used the Hill Country Trinity Aquifer groundwater availability model to estimate that 72.8 × 10⁶ M³/yr (59,000 acre-ft/yr) recharged the Edwards Aquifer from the Trinity Aquifer as interformational flow, based on conditions representative of 1975. The Hill Country portion of the Trinity Aquifer only extends to the Dry Frio/Frio River watersheds to the west, excluding the West Nueces/Nueces River watersheds. The Hill Country Trinity Aquifer model by Jones et al. (2011) calculated that ~60% of discharge by the Trinity Aquifer is to streams, springs, and reservoirs, and 35% discharges through cross-formational flow to the Edwards (Balcones Fault Zone) Aquifer. Fratesi et al. (2015) refined the estimates for interformational flow by explicitly including the contributing zone of the Edwards Aquifer in their model; nonetheless, significant uncertainty remains in quantifying interformational flow between the Trinity and Edwards Aquifers. In recognition of this limitation, the Edwards Aquifer Authority has established a multiyear program to resolve this area of uncertainty.

Most numerical models of the Edwards Aquifer and related systems have time steps or stress periods from 1 mo to 1 yr. The greatest limitation when attempting to reduce the length of time steps or stress periods is the lack of pumping data or precipitation information with sufficient temporal or spatial resolution. One aspect of this limitation that will remain, even if current pumping and precipitation data resolution are improved, is the lack of historical data with similar improved resolution. Thus, improvement in simulating historical periods such as the drought-of-record in south-central Texas (i.e., 1950s) will always retain a certain level of uncertainty due to lack of quality data.

Perhaps the greatest source of uncertainty in models of regions with semiarid climates (i.e., Edwards-Trinity Aquifer) is recharge. Whereas models whose domains receive average annual precipitation rates in excess of 500 mL/yr can be expected to experience distributed or direct recharge, areas that receive precipitation <~400 mL/yr typically rely on runoff that is focused to river and stream channels for recharge (Green et al., 2014). Quantification of recharge for areas where precipitation is limited and highly variable, such as the Edwards Plateau, introduces significant uncertainty into numerical groundwater models. Advances to improve recharge estimates and measurements using regional and local investigations and the inclusion of remotely sensed data could greatly improve numerical groundwater model accuracy.

We wish to acknowledge technical reviews by Alan Dutton, Yongli Gao, and Jack Sharp. Their comments are greatly appreciated.