UNITED STATES

MEXICO

Model boundary

Active 

model-area

Lake

 Travis

Blanco River

Canyon 

Lake

Creek

Medina

Lake

Pedernales

River

San

Marcos

River

Nueces

River

Frio

River

West

Nueces

River

Medina

River

San Antonio River

Cibolo

Creek

Leona

River

Guadalupe River

Sabinal River

Colorado River

River Frio

RIO BRAVO

RIO GRANDE

TEXAS

LLANO COUNTY

BURNET

COUNTY

BLANCO

COUNTY

WILLIAMSON

COUNTY

BASTROP

COUNTY

TRAVIS

COUNTY

HAYS COUNTY

COMAL

COUNTY

KENDALL

COUNTY

GILLESPIE

COUNTY

EDWARDS COUNTY KERR COUNTY 

REAL

COUNTY 

BANDERA

COUNTY

BEXAR COUNTY

GUADALUPE

COUNTY

CALDWELL

COUNTY

GONZALES

COUNTY

WILSON COUNTY

MEDINA COUNTY UVALDE COUNTY

KINNEY

COUNTY 

MAVERICK

COUNTY 

ZAVALA COUNTY FRIO COUNTY

ATASCOSA

COUNTY

LA SALLE COUNTY

DIMMIT COUNTY 

TEXAS

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing density Increasing density Increasing density

PILOT POINT MULTIPLIER PARAMETERS

Specific storage

Specific storage

Horizontal hydraulic 

conductivity

Vertical hydraulic 

conductivity

Horizontal hydraulic 

conductivity

Initial dissolved-solids

concentration

Initial dissolved-solids

concentration

Initial dissolved-solids

concentration

U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2015–5081

Prepared in cooperation with the San Antonio Water System 

Updated Numerical Model with Uncertainty Assessment of 

1950–56 Drought Conditions on Brackish-Water Movement 

within the Edwards Aquifer, San Antonio, Texas
Cover illustration: modified versions of figures 1 and 22.
Updated Numerical Model with 

Uncertainty Assessment of 1950–56 

Drought Conditions on Brackish-Water 

Movement within the Edwards Aquifer, 

San Antonio, Texas

By Linzy K. Brakefield, Jeremy T. White, Natalie A. Houston, and 

Jonathan V. Thomas

Prepared in cooperation with the San Antonio Water System

Scientific Investigations Report 2015–5081

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

U.S. Geological Survey, Reston, Virginia: 2015

For more information on the USGS—the Federal source for science about the Earth, its natural and living 

resources, natural hazards, and the environment—visit http://www.usgs.gov or call 1–888–ASK–USGS.

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod/.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials 

as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.

Suggested citation:

Brakefield, L.K., White, J.T., Houston, N.A., and Thomas, J.V., 2015, Updated numerical model with uncertainty 

assessment of 1950–56 drought conditions on brackish-water movement within the Edwards aquifer, San Antonio, Texas: U.S. Geological Survey Scientific Investigations Report 2015–5081, 54 p., http://dx.doi.org/10.3133/

sir20155081.

ISSN 2328-0328 (online)
iii

Acknowledgments

The authors thank Richard Lindgren (U.S. Geological Survey, retired) for his insight and knowledge of previously developed MODFLOW models of the Edwards aquifer and his willingness to 

provide assistance.
iv

Contents

Acknowledgments........................................................................................................................................iii

Abstract...........................................................................................................................................................1

Introduction.....................................................................................................................................................2

Purpose and Scope ..............................................................................................................................4

Background and Previous Investigations.........................................................................................4

Updated Conceptualization and Numerical Model of the Edwards Aquifer........................................5

Updated Conceptualization of Hydrogeology for Numerical Model ............................................6

Review of Groundwater-Flow System...............................................................................................7

Water-Quality Zones....................................................................................................................7

Aquifer Hydraulic and Transport Properties ...........................................................................9

Model Development ...........................................................................................................................11

History-Matching Phase....................................................................................................................11

Discretization..............................................................................................................................13

Hydrologic Boundary Conditions ............................................................................................13

Lateral Model Boundaries ..............................................................................................13

Groundwater Recharge ...................................................................................................17

Groundwater Discharge ..................................................................................................17

Well Discharge.........................................................................................................19

Spring Discharge .....................................................................................................19

Initial Conditions ........................................................................................................................19

Initial Heads.......................................................................................................................19

Initial Dissolved-Solids Concentrations........................................................................19

Requirements for SEAWAT.......................................................................................................22

Calibration through History Matching....................................................................................22

Parameterization ..............................................................................................................23

Hydraulic and Transport Parameters ...................................................................23

Model Sources and Sinks ......................................................................................27

Initial Dissolved-Solids Concentrations...............................................................32

Historical Observation Dataset ......................................................................................32

Groundwater Head Observations .........................................................................35

Spring Discharge Observations ............................................................................35

Dissolved-Solids Concentration............................................................................35

Constraints on Total Mass of Dissolved Solids at Production Wells ..............35

Goodness of Fit..................................................................................................................38

Parameter Uncertainty—Schur’s Complement for Linear-Based Conditional Uncertainty 

Propagation ...........................................................................................................................38

Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement within 

the Edwards Aquifer .....................................................................................................................42

Predictive Phase ................................................................................................................................43

Drought Recharge Conditions 1950–56 ..................................................................................43

Withdrawals by Wells ...............................................................................................................43
v

Prediction Uncertainty—Schur’s Complement for Linear-Based Conditional Uncertainty 

Propagation............................................................................................................................43

Predictive Uncertainty Results.........................................................................................................45

Model Limitations ........................................................................................................................................46

Summary........................................................................................................................................................47

References Cited..........................................................................................................................................49

Figures

1 Map showing location of model area, Edwards aquifer segments, depositional 

provinces, groundwater divides, and physiographic regions, San Antonio region, 

Texas...............................................................................................................................................3

2. Diagrammatic northwest-to-southeast section showing hydrogeologic framework 

and generalized groundwater-flow directions, San Antonio region, Texas .......................6

3. Graph showing correlation of Cretaceous stratigraphic units and hydrogeologic 

units, relative permeability and model layers in the Edwards aquifer model area, 

San Antonio region, Texas ..........................................................................................................8

4. Map showing recharge, freshwater, brackish-water transition, and saline zones 

of the Edwards aquifer, San Antonio region, Texas..............................................................10

5. Map showing active model area and boundary conditions, San Antonio region, 

Texas.............................................................................................................................................12

6. Maps showing faults and altitude of top of Edwards aquifer active model area, 

and thickness of Edwards aquifer active model area, San Antonio region, Texas .........14

7. Maps showing altitude of top, and thickness of the middle part of the 

Edwards aquifer (model layer 4), San Antonio region, Texas .............................................15

8. Map showing number of active model layers per cell, San Antonio region, Texas. .......16

9. Map showing lines of equal dissolved-solids concentrations, dissolved-solids 

concentration locations of measured or estimated values based on data collected 

from 1930 through 2009, and distribution of interpolated dissolved-solids concentrations, San Antonio region, Texas..............................................................................................18

10. Map showing average annual groundwater withdrawals during 1999–2009, and 

the production wells of predictive interest near the brackish-water transition zone, 

San Antonio region, Texas ........................................................................................................20

11. Map showing locations of pilot points where hydraulic properties are adjusted 

during history matching, San Antonio region, Texas............................................................24

12. Maps showing the history-matched best-fit distributions of hydraulic conductivities 

of the Edwards aquifer, San Antonio region, Texas. Horizontal hydraulic conductivity 

of the upper part of the Edwards aquifer. Vertical hydraulic conductivity of the middle 

part of the Edwards aquifer. Horizontal hydraulic conductivity of the lower part of the 

Edwards aquifer..........................................................................................................................26

13. Maps showing the history-matched best-fit distributions of specific storage of the 

Edwards aquifer, San Antonio region, Texas for the upper part of the Edwards 

aquifer, and the lower part of the Edwards aquifer..............................................................28

14. Maps showing the history-matched best-fit distributions of effective porosity of the 

Edwards aquifer, San Antonio region, Texas for the upper part of the Edwards 

aquifer, and the lower part of the Edwards aquifer..............................................................29
vi

15. Map showing the history-matched best-fit distributions of recharge multipliers for 

the surface-water basins overlying the Edwards aquifer, San Antonio region, Texas. 

Recharge multipliers for the wetter months of April–June, September–November. 

Recharge multipliers for the drier months of January–March, July, August, and 

December.....................................................................................................................................30

16. Graph showing the history-matched best-fit annual recharge multipliers, 

1999–2009 .....................................................................................................................................32

17. Maps showing inflow of groundwater from the Trinity aquifer (flux) into the overlying 

Edwards aquifer, San Antonio region, Texas. Initial flux. History-matched bestfit flux ............................................................................................................................................33

18. Maps showing the history-matched best-fit distributions of initial dissolved-solids 

concentrations of the Edwards aquifer, San Antonio region, Texas for the upper 

part of the Edwards aquifer, the middle part of the Edwards aquifer, and the 

lower part of the Edwards aquifer...........................................................................................34

19. Map showing locations of observation data used for history matching (groundwater 

withdrawals, springs, dissolved-solids concentrations, and groundwater heads), 

San Antonio region, Texas, 1999–2009 ....................................................................................36

20. Graphs showing observed values compared to simulated values and residuals for 

groundwater heads, spring discharges, and dissolved-solids concentrations, San 

Antonio region, Texas.................................................................................................................37

21. Graphs showing simulated and observed time-series of Comal Springs discharge, 

San Marcos Springs discharge, and Bexar County index well J-17 groundwater 

head, San Antonio region, Texas .............................................................................................39

22. Graphs showing distributions of prior (before history matching) and posterior (after 

history matching) parameter uncertainty...............................................................................40

23. Graphs showing groundwater flows (fluxes) and cumulative volumes for the historymatched model (1999–2009) and the predictive model (1950–56), San Antonio region, 

Texas. Positive values represent groundwater gains, and negative values represent 

groundwater losses....................................................................................................................44

24. Graph showing the upper bound of the 95-percent credible interval for the prior 

and posterior predictive uncertainty of dissolved-solids concentration change 

during the 1950–56 predictive model period at production wells of interest....................45

Tables

1. The history-matched best-fit Drain (DRN) package altitude and conductance parameters that control spring discharge from the model, San Antonio region, Texas ............32
vii

Conversion Factors

Inch/Pound to International System of Units

Multiply By To obtain

Length

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Volume

cubic foot (ft3) 0.02832 cubic meter (m3) 

acre-foot (acre-ft) 1,233 cubic meter (m3)

Flow rate

acre-foot per year (acre-ft/yr) 1,233 cubic meter per year (m3/yr)

cubic foot per day (ft3/d) 0.02832 cubic meter per day (m3/d)

Density

pound per cubic foot (lb/ft3) 16.02 kilogram per cubic meter (kg/m3)

Hydraulic conductivity

foot per day (ft/d) 0.3048 meter per day (m/d)

DRN (Drain Package) conductance

foot squared per day (ft2/d) 0.09290 meter squared per day (m2/d) 

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as 

°F = (1.8 × °C) + 32.

Supplemental Information

Specific conductance is given in microsiemens per centimeter at 25 degrees Celsius (μS/cm at 

25 °C).

Concentrations of chemical constituents in water are given in milligrams per liter (mg/L).

Datums

Vertical coordinate information is referenced to the National Geodetic Vertical Datum of 1929 

(NGVD 29).

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83).

Altitude, as used in this report, refers to distance above the vertical datum.
viii

Abbreviations

BTN Basic Transport Package (in MT3DMS)

DRN Drain Package (in MODFLOW)

EAA Edwards Aquifer Authority

GIS geographic information system

HFB Horizontal-Flow Barrier Package (in MODFLOW)

MODFLOW A modular three-dimensional finite-difference ground-water flow model

MT3DMS Modular 3-D multi-species transport model for simulation of advection, 

dispersion, and chemical reactions of contaminants in groundwater systems

MUD Fort Clark Municipal Utility District

NWIS National Water Information System

PEST Model-independent Parameter ESTimation code 

RCH Recharge Package (in MODFLOW)

RDM regional dense member

SAWS San Antonio Water System

SEAWAT A computer program for simulation of three-dimensional variable-density 

groundwater flow

SP spontaneous potential

SSM Source/Sink Mixing Package (in MT3DMS)

SSP static self-potential

TCEQ Texas Commission on Environmental Quality

TWDB Texas Water Development Board

USGS U.S. Geological Survey

WEL Well Package (in MODFLOW)
Abstract

In 2010, the U.S. Geological Survey, in cooperation with 

the San Antonio Water System, began a study to assess the 

brackish-water movement within the Edwards aquifer (more 

specifically the potential for brackish-water encroachment into 

wells near the interface between the freshwater and brackishwater transition zones, referred to in this report as the transition-zone interface) and effects on spring discharge at Comal 

and San Marcos Springs under drought conditions using a 

numerical model. The quantitative targets of this study are 

to predict the effects of higher-than-average groundwater withdrawals from wells and drought-of-record rainfall conditions 

of 1950–56 on (1) dissolved-solids concentration changes at 

production wells near the transition-zone interface, (2) total 

spring discharge at Comal and San Marcos Springs, and 

(3) the groundwater head (head) at Bexar County index 

well J-17. The predictions of interest, and the parameters 

implemented into the model, were evaluated to quantify their 

uncertainty so the results of the predictions could be presented 

in terms of a 95-percent credible interval.

The model area covers the San Antonio and Barton 

Springs segments of the Edwards aquifer; the history-matching effort was focused on the San Antonio segment. A previously developed diffuse-flow model of the Edwards aquifer, 

which forms the basis for the model in this assessment, is 

primarily based on a conceptualization in which flow in the 

aquifer is predominately through a network of numerous small 

fractures and openings. Primary updates to this model include 

an extension of the active area downdip, a conversion to an 

8-layer SEAWAT variable-density flow and transport model to 

simulate dissolved-solids concentration effects on water density, history matching to 1999–2009 conditions, and parameter 

estimation in a highly parameterized context using automated 

methods in PEST (a model-independent Parameter ESTimation code). 

In addition to the best-fit parameter values derived from 

history matching, the uncertainty of model parameters was 

also estimated by using linear uncertainty analysis. Comparison of “prior” (before history matching) and “posterior” 

(after history matching) variances of parameters indicate that 

the information within the observation dataset used for history matching informs many parameters. The concentration 

threshold parameters were well-informed by the observation 

dataset as their posterior distributions were much narrower 

than their prior distributions. The transition-zone scaling 

parameters of hydraulic conductivity, effective porosity, and 

specific storage were all informed by the observation dataset, 

as evidenced by the difference between the prior and posterior 

variances. Saline-zone scaling parameters, alternatively, were 

not informed by the observation dataset for effective porosity 

and specific storage. Resulting posterior drier-month, wettermonth, and annual recharge multiplier parameter variances are 

important to understanding how well recharge is estimated and 

implemented within the model. The shifts of the posterior distributions left and right indicate that there were zones where 

less or more water was needed in the model. The widths of the 

distributions were not decreased substantially, indicating that 

many of the best-fit recharge parameters are not statistically 

different from the initial values specified in the history-matching effort. Recharge from rainfall is the driving force behind 

groundwater flow and heads in the aquifer; therefore, an 

increase in understanding of this process would benefit model 

development by potentially decreasing the uncertainty of this 

parameter. The history-matching effort was most helpful in 

informing the parameters in the model that control discharge 

at springs, namely, the spring orifice (drain) altitude and drain 

conductance parameters for each spring.

The uncertainty assessment of the predictive model (a 

hypothetical recurrence of 1950–56 drought conditions and 

higher-than-average groundwater withdrawals from wells) 

provided insights into the potential effects of these conditions 

on dissolved-solids concentration changes at production wells 

near the transition-zone interface, discharges at Comal and San 

Marcos Springs, and heads at Bexar County index well J-17. 

Results at the 25 production wells near the transition-zone 

Updated Numerical Model with Uncertainty Assessment 

of 1950–56 Drought Conditions on Brackish-Water 

Movement within the Edwards Aquifer, San Antonio, 

Texas 

By Linzy K. Brakefield, Jeremy T. White, Natalie A. Houston, and Jonathan V. Thomas
2 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

interface indicate that the uncertainty of model input parameters based on expert knowledge yielded an upper bound of 

the 95-percent credible interval of dissolved-solids concentrations that exceeds the secondary drinking water standards of 

1,000 milligrams per liter (mg/L) of the Texas Commission 

on Environmental Quality (TCEQ) for many wells. However, 

the history-matching process provided key information to 

inform prediction-sensitive model parameters and therefore, 

contributed to a substantial decrease of the upper bound of the 

95-percent credible interval to below the secondary drinking 

water standards. Reductions in dissolved-solids concentration 

changes were on the order of 400 mg/L to 1,300 mg/L. The 

reduction in uncertainty in regards to this prediction implies 

that this prediction of dissolved-solids concentration change 

can be made with some certainty using this current model 

and that those parameters that control this prediction are 

informed by the observation dataset. Even though predictive 

uncertainty was reduced for this prediction, dissolved-solids 

concentration changes were still greater than zero, indicating 

a minimal increase in concentration at these 25 production 

wells during the 7-year simulation period is likely. However, 

this minimal concentration increase indicates a small potential for movement of the brackish-water transition zone near 

these wells during the 7-year simulation period of drought-ofrecord (1950–56) rainfall conditions with higher-than-average 

groundwater withdrawals by wells.

Predictive results of total spring discharge during the 

7-year period, as well as head predictions at Bexar County 

index well J-17, were much different than the dissolved-solids 

concentration change results at the production wells. These 

upper bounds are an order of magnitude larger than the actual 

prediction which implies that (1) the predictions of total 

spring discharge at Comal and San Marcos Springs and head 

at Bexar County index well J-17 made with this model are not 

reliable, and (2) parameters that control these predictions are 

not informed well by the observation dataset during historymatching, even though the history-matching process yielded 

parameters to reproduce spring discharges and heads at these 

locations during the history-matching period. Furthermore, 

because spring discharges at these two springs and heads at 

Bexar County index well J-17 represent more of a cumulative effect of upstream conditions over a larger distance (and 

longer time), many more parameters (with their own uncertainties) are potentially controlling these predictions than the 

prediction of dissolved-solids concentration change at the 

prediction wells, and therefore contributing to a large posterior 

uncertainty. 

Introduction

The Edwards aquifer is one of the most prolific aquifers 

in the world and supplies water to users in several counties in 

south-central Texas (Maclay, 1995). The largest municipality 

that uses the Edwards aquifer for water supply is San Antonio, 

Tex., the seventh most populous city in the U.S. (U.S. Census Bureau, 2014). Residents in south-central Texas rely on 

water from the Edwards aquifer for residential, recreational, 

industrial, and agricultural uses (Lindgren and others, 2004). 

In the San Antonio and Austin areas, the Edwards aquifer 

was deemed a sole-source aquifer by the U.S. Environmental 

Protection Agency (1975, 1988, and 2015), and several endangered and threatened species are sustained by groundwater 

discharged at Comal, San Marcos, and Barton Springs (fig. 1). 

The Edwards aquifer consists of regionally extensive 

carbonate rocks that dip from the outcrop region within the 

Edwards Plateau and the Balcones fault zone toward the south 

and southeast beneath the Gulf Coastal Plain (fig. 1). The 

southernmost extent of the freshwater zone of the Edwards 

aquifer is defined by a dissolved-solids concentration of 

1,000 milligrams per liter (mg/L), a concentration that serves 

as the interface between the freshwater zone of the aquifer 

and the brackish-water transition zone that consists of water 

with dissolved-solids concentrations ranging from 1,000 to 

10,000 mg/L. The brackish-water transition zone discussed in 

this report has been referred to as the freshwater/saline-water 

transition zone in other reports (Lambert and others, 2010; 

Thomas and others, 2012). The dissolved-solids concentrations of 1,000 to 10,000 mg/L are within the range termed 

“brackish” by Freeze and Cherry (1979). In this report, the 

interface between the freshwater and brackish-water transition 

zones is referred to as the transition-zone interface, and the 

brackish-water transition zone is referred to as the transition 

zone. The part of the Edwards aquifer downdip of the transition zone is referred to as the saline zone, which is characterized by water with dissolved-solids concentrations in excess of 

10,000 mg/L. 

Because the Edwards aquifer is an important water 

resource for the area, there are concerns about a repeat, or 

exceedance of, the drought-of-record (1950–56) conditions. 

Recent (1999–2009) groundwater withdrawal rates in the 

Edwards aquifer are much greater than they were during the 

1950s (Lindgren and others, 2004) and could potentially cause 

movement of the transition zone, cause production wells to 

draw nonpotable water, or result in sustained reductions in 

spring discharges from Comal and San Marcos Springs 

(fig. 1). Since 2011, drought conditions have existed in southcentral Texas (Texas Water Development Board, 2014), further 

underscoring the need to better understand drought effects 

on potential brackish-water movement, groundwater levels 

(heads), and spring discharges. 

Many studies have characterized the transition zone to 

better understand this zone between the freshwater and saline 

water zones of the Edwards aquifer (for example, Garza, 

1962; Groschen, 1994; Groschen and Buszka, 1997; Pavlicek and others, 1987; Perez, 1986; Schultz, 1992; Schultz, 

1993; Schultz, 1994). Several wells constructed in transects 

across the transition zone have been the subject of studies by 

the U.S. Geological Survey (USGS) in cooperation with the 
Introduction 3

UNITED STATES

MEXICO

0 10 20 30 KILOMETERS

0 10 20 30 MILES

29°

30°

31°

98°

99°

100°

Base from U.S. Geological Survey and

City of Austin digital data, 2010

Albers equal area projection 

Texas Centric Mapping System

North American Datum of 1983

MAVERICK BASIN

DEVILS RIVER

TREND

SAN MARCOS

San Pedro PLATFORM

Springs

San Antonio Springs

Hueco Springs

Comal

Springs

Las Moras Springs

Leona Springs

San Marcos

Springs

Barton

Springs 

Lake

 Travis

Blanco River

Canyon 

Lake

Creek

Medina

Lake

Pedernales

River

San

Marcos

River

Nueces

River

Frio

River

West

Nueces

River

Medina

River

San Antonio River

Cibolo

Creek

Leona

River

Guadalupe River

Sabinal River

Colorado River

River Frio

RIO BRAVO

RIO GRANDE

TEXAS

COAHUILA

LLANO COUNTY

BURNET

COUNTY

BLANCO

COUNTY

WILLIAMSON

COUNTY

BASTROP COUNTY

TRAVIS

COUNTY

HAYS COUNTY

COMAL

COUNTY

KENDALL

COUNTY

GILLESPIE

COUNTY

EDWARDS COUNTY KERR COUNTY 

REAL

COUNTY 

BANDERA

COUNTY

BEXAR COUNTY

GUADALUPE

COUNTY

CALDWELL COUNTY

GONZALES

COUNTY

WILSON COUNTY

MEDINA COUNTY UVALDE COUNTY

KINNEY

COUNTY 

MAVERICK

COUNTY 

ZAVALA COUNTY FRIO COUNTY

ATASCOSA

COUNTY

LA SALLE COUNTY

DIMMIT COUNTY 

Austin

Buda

New Braunfels

San Marcos

Brackettville

Uvalde

Hondo

San 

Antonio 

Kyle

G U L F C O A S T A L P L A I N

E D W A R D S P L A T E A U

B A L C O N E S F A U L T Z O N E

STUART CITY

REEF TREND

TEXAS

Model boundary 

Active

model area

Bexar County index well (J–17)

EXPLANATION

Model area

 Active regional-scale

 model area

 San Antonio segment of

 the Edwards aquifer

 Barton Springs segment

 of the Edwards aquifer

 Inactive regional-scale

 model area

Model boundary

Active model-area boundary

Depositional province boundary (modified from

 Maclay, 1995, fig. 2)

Boundary of Balcones fault zone

Balcones escarpment (modified from

 Abbott and Woodruff, 1986, fig. 1)

Recharge zone boundary (modified from 

 Puente, 1978)

Groundwater divide

Groundwater divide (extrapolated)

Spring—Tail points in direction of flow

Figure 1 Location of model area, Edwards aquifer segments, depositional provinces, groundwater divides, and physiographic 

regions, San Antonio region, Texas (modified from Lindgren and others [2011a, fig. 1]).
4 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

San Antonio Water System (SAWS). The objectives of these 

studies have included gaining a better understanding of the 

hydraulic dynamics between the freshwater and transition 

zones, and insights into the structure and properties of the 

aquifer in the transition zone at the local scale (Lambert and 

others, 2009; Lambert and others, 2010; Thomas and others, 

2012). 

In 2010, the USGS, in cooperation with SAWS, began 

a study to assess the brackish-water movement within the 

Edwards aquifer (more specifically the potential for brackish-water encroachment into wells near the transition-zone 

interface) and effects on discharge at Comal and San Marcos 

Springs under drought conditions using a numerical model. 

Twenty-five production wells near the transition-zone interface were selected to evaluate encroachment of brackish 

water updip from the transition zone into the freshwater zone. 

Discharges at Comal and San Marcos Springs are a concern 

because of the potential ecological effects on riparian habitat 

that sustains several endangered species; the vitality of the 

riparian habitat is dependent upon flows from these springs. 

Groundwater heads (heads) at the Bexar County index well 

J-17 (hereinafter referred to as “ well J-17”) are an indicator 

of groundwater availability in the Edwards aquifer; heads at 

this well have been measured since 1932 (Edwards Aquifer 

Authority, 2015) and a statistical correlation between heads at 

this well and discharge at Comal Springs has been established 

(Miller and Long, 2006). Therefore, it was also of interest to 

determine effects from drought conditions on heads at this 

monitoring well.

Purpose and Scope

This report describes a variable-density groundwater 

flow and transport numerical model constructed by using 

the computer program SEAWAT version 4. Modeling was 

completed in two phases: the history-matching phase and 

the predictive phase. In addition to modeling, an uncertainty 

assessment was done to quantify the combined effects of 

rainfall amounts identical to those of the drought of record 

(1950–56) and recent (1999–2009) groundwater withdrawals 

from the Edwards aquifer on the potential movement of the 

transition zone in the San Antonio segment of the Edwards 

aquifer, potential effects on discharge at Comal and San Marcos Springs, and potential effects on heads at well J-17. The 

model that was developed is an updated numerical model of 

the San Antonio and Barton Springs segments of the Edwards 

aquifer that was based on the diffuse-flow model of Lindgren 

(2006). Major updates to this model, as described in this 

report in the history-matching phase, include simulation of 

the Edwards aquifer farther downdip than in previous models, 

the addition of more model layers to better represent aquifer 

heterogeneities and anisotropy, and an updated simulation 

(1999–2009) for history matching. The model was history 

matched in a highly parameterized context using PEST, a 

model-independent Parameter ESTimation code that allows 

for model calibration, parameter estimation, and calculation 

of parameter and predictive uncertainties (Doherty, 2005; 

Doherty and Hunt, 2010; Doherty and others, 2010b). The 

focus of the predictive model created in the predictive phase 

(hypothetical scenario of drought-of-record rainfall combined 

with recent groundwater withdrawals), and therefore the 

history-matching effort, was on the freshwater and transition zones of the confined part of the San Antonio segment 

of the Edwards aquifer (fig. 1). Linear uncertainty analysis 

was completed by using PEST to quantify the uncertainty of 

model parameters and predictions. The uncertainty assessment 

focused on (1) dissolved-solids concentrations changes in 

target production wells near the transition-zone interface, 

(2) spring discharge at Comal and San Marcos Springs, and 

(3) heads at well J-17. 

Background and Previous Investigations

Several numerical groundwater-flow models have been 

constructed for the San Antonio and Barton Springs segments 

of the Edwards aquifer. Lindgren and others (2009) described 

and evaluated six regional-scale groundwater-flow models. 

One model discussed in Lindgren and others (2009) was the 

one-layer (vertically averaged) diffuse-flow model (Lindgren, 

2006). This model forms the basis of the model developed in 

this study. The diffuse-flow model is based primarily on the 

conceptualization in which water in the aquifer flows predominantly through a network of numerous small fractures 

and openings. The diffuse-flow model, as well as the previous 

one-layer conduit-flow model (Lindgren and others, 2004), 

were modified by Lindgren and others (2011b) to have a finer 

horizontal discretization, an additional second layer, and to 

simulate conditions of 2001–2003. The focus of the previous 

MODFLOW modeling efforts were to calibrate these models to steady-state and transient periods to better understand 

groundwater budgets, regional flow, and transmissivity and 

storativity distributions. A local-scale MODFLOW model was 

also developed to investigate the contributing zone to a public 

supply well in San Antonio (Lindgren and others, 2011a) 

using particle-tracking. All of these models assume a uniform, 

freshwater density for the Edwards aquifer, that extends to the 

southern boundary (approximated at the 10,000-mg/L contour 

of equal dissolved-solids concentration [hereinafter referred 

to as the “10,000-mg/L dissolved-solids contour”]), and do 

not account for variable-density conditions that exist in the 

aquifer. Furthermore, the regional-scale models are vertically 

simplified representations of the complex anisotropic and 

heterogeneous sloping aquifer. The only previously published 

model simulating potential transition-zone encroachment is a 

one-dimensional solute-transport model (Perez, 1986). Results 

from the simulations and sensitivity analysis completed by 

Perez indicated that to accurately simulate encroachment of 

saline water into the freshwater zone, more information was 
Updated Conceptualization and Numerical Model of the Edwards Aquifer 5

needed to characterize groundwater heads, transmissivity, 

porosity, and dispersivity in the transition zone separating the 

freshwater and saline water zones. Since the one-dimensional 

simplified model of Perez (1986), improved modeling tools, 

increased computational power, additional data, and more scientific knowledge exist, allowing for more rigorous simulation 

of variable-density flow in the freshwater and transition zones.

Several field-based assessments have been completed 

to better characterize the transition zone by using water-level 

(head), lithologic, temperature, and geophysical log data collected at wells that traverse the transition zone (transect wells) 

(Lambert and others, 2010; Lambert and others, 2009; Thomas 

and others, 2012). Analysis of equivalent freshwater-head 

gradients at these transect wells indicate that, although some 

potential for flow between the freshwater and transition zones 

exists, this potential is highly variable across the region. Furthermore, the transition zone is stable laterally and vertically 

and in time, and is likely to remain stable (Lambert and others, 

2010). A previous local-scale field investigation by Groschen 

(1994) was focused on gaining a better understanding of the 

hydrology of the aquifer near the downdip extent of freshwater 

and also was focused on gaining a better understanding of the 

evolution of water chemistry. Groschen (1994) also presented 

a thorough literature review of current understanding of the 

transition zone and a summary about saline-zone groundwater 

flow and geochemical alteration of the aquifer from freshwater 

flows and dedolomitization [the dissolution of dolomite and 

net precipitation of calcite (Deike, 1990)] through time.

Contours of dissolved-solids concentrations in the 

Edwards aquifer (Schultz 1992, 1993, 1994; A.L. Schultz, 

private consultant, written commun., 2000) were mapped (two 

dimensions) based on geophysical-log estimates of dissolvedsolids concentrations and water-quality data. To date, this 

map is the most widely accepted understanding of the complicated three-dimensional distribution of dissolved solids in the 

Edwards aquifer, although the map contours are a two-dimensional interpretation and do not depict changes in dissolved 

solids at depth. The mapped 10,000-mg/L dissolved-solids 

contour (A.L. Schultz, private consultant, written commun., 

2000) has been used as the southern boundary in models simulating freshwater flows (Lindgren and others, 2004; Lindgren, 

2006, Lindgren and others, 2011b).

Groschen and Buszka (1997) focused on geochemical 

characterization of the saline zone of the Edwards aquifer. The 

primary findings based on geochemical and isotope analysis 

were that water of the saline zone can be separated into two 

distinct hydrological and geochemical regimes: the shallower 

regime is mostly meteoric water that has been recharged 

within recent geologic time and the deeper regime is hydrologically stagnant when compared to the updip regime. These 

two regimes exist, in part, because faults function as barriers 

to downdip flow of recharge water and likely, to updip flow of 

saline water (Groschen and Buszka, 1997). 

Updated Conceptualization and 

Numerical Model of the Edwards 

Aquifer

The conceptualization of the Edwards aquifer and implementation of this conceptualization into a numerical model 

was based on previous modeling efforts (Lindgren and others, 

2004; Lindgren, 2006; Lindgren and others, 2011b). In previous USGS models, the Edwards aquifer was simulated on a 

regional scale by using a model with a single layer (Lindgren 

and others, 2004; Lindgren, 2006) or a model with two layers 

(Lindgren and others, 2011b). Previous models simulated a 

freshwater zone and brackish-water zone of groundwater containing as much as 10,000 mg/L of dissolved solids, although a 

uniform freshwater density was assumed throughout the model 

domain (Lindgren and others, 2004; Lindgren, 2006; Lindgren 

and others, 2011b). Because variations in density can potentially affect groundwater flows, the uniform freshwater-density 

diffuse-flow model (Lindgren, 2006) was converted into an 

eight-layer coupled variable-density flow and transport model. 

The USGS’s SEAWAT version 4 computer program (Langevin 

and others, 2008) was used to develop the updated model. 

This updated conceptualization and model implementation 

was completed to include a larger active downdip extent of the 

model, include simulation of density-driven flow and transport 

processes (specifically dissolved-solids concentration effects 

on density) and to better represent the vertical and horizontal 

anisotropy and heterogeneity of the aquifer. 

Uncertainty in the parameters implemented into the 

model and the model predictions of interest were quantified 

by using a linear uncertainty analysis approach (Doherty 

and others, 2010b) with 95-percent credible intervals. The 

linear uncertainty analysis used in this assessment is based 

on Bayes’ equation (Doherty and others, 2010b). The uncertainty analysis yields estimates of credible intervals instead 

of the more common confidence intervals (used in frequentist statistics). Gelman and others (2014, p. 3) explain that a 

credible interval (also known as Bayesian interval or Bayesian probability interval) “for an unknown quantity of interest can be directly regarded as having a high-probability of 

containing the unknown quantity, in contrast to a frequentist 

(confidence) interval, which may strictly be interpreted only in 

relation to a sequence of similar inferences that may be made 

in repeated practice.” In this study, the term Gaussian distribution is a probability density function and the uncertainty is 

measured through estimation of 95-percent credible intervals. 

The prior is a multivariate Gaussian distribution that describes 

parameter uncertainty defined solely on the basis of the innate 

variability of parameters (Doherty and others, 2010b). The 

prior must be quantified by the modeler (Doherty and others, 2010b) based on their expert knowledge of the system 

(Doherty, 2005; Fienen and others, 2013). The posterior is a 

multivariate Gaussian distribution, which combines the prior
6 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

with information contained in the observation dataset used 

for history matching. History matching is often referred to as 

model calibration, but the term history matching has been used 

extensively in reservoir modeling in the petroleum industry 

(Oliver and others, 2008) and is chosen for this report because 

it more “clearly conveys the essence of the modeling process 

than does the term calibration (Konikow and Bredehoeft, 

1992).” Gelman describes the posterior as the “conditional 

probability distribution of the unobserved quantities of ultimate interest, given the observed data” (Gelman and others, 

2014). Model parameters and predictions of interest given as 

credible intervals, rather than one estimate derived from the 

“best fit” of the model parameters to the data, where “best fit” 

refers to the minimum of the weighted least-squares objective 

function, are more defensible and more meaningful. Furthermore, uncertainty estimates (credible intervals) are indicative 

of a model’s reliability in making a prediction (Refsgaard and 

others, 2007).

Updated Conceptualization of Hydrogeology for 

Numerical Model

Details on the geology and hydrogeologic setting of the 

freshwater extent have been extensively discussed in several 

modeling and geology reports (Lindgren and others, 2004; 

Lindgren, 2006; Maclay, 1995; Rose, 1972). Lower Cretaceous limestone and dolomite compose the karstic Edwards 

aquifer (Hovorka and others, 1996; Rose, 1972), which 

extends downdip of the freshwater zone within the Balcones 

Fault Zone to the Stuart City Reef Trend (fig. 1) (Groschen 

and Buszka, 1997; Oetting and others, 1996; Waite, 2008) at 

depths greater than 10,000 feet below the National Geodetic 

Vertical Datum of 1929 (Rose, 1972). The unconfined outcrop 

of the Edwards aquifer in the northern part of the study area is 

the recharge zone (figs. 1–2). South and southeast of the outcrop, the aquifer becomes confined and is overlain by younger, 

lower-permeability rocks (figs. 1–2).

NW

Spring

Not to scale

SE

Recharge

Discharge

Spring

Supply

EXPLANATION well

Upper confining unit

Edwards aquifer

 Unconfined aquifer

 Confined aquifer

Trinity aquifer

Limestone

Sand

Underlying rock units

Approximate freshwater/saline-

 water interface

Generalized conduit system

Generalized direction of 

 groundwater flow

Fault—Arrows indicate direction

 of movement

Figure 2. Diagrammatic northwest-to-southeast section showing hydrogeologic framework and 

generalized groundwater-flow directions, San Antonio region, Texas (modified from Barker and Ardis, 

1996, plates 1 and 3).
Updated Conceptualization and Numerical Model of the Edwards Aquifer 7

The Edwards aquifer is characterized by three different depositional provinces (figs. 1 and 3): the San Marcos 

Platform in the eastern part of the study area, the Maverick 

Basin in the western part, and the Devils River Trend in the 

west-central part (Lozo and Smith, 1964; Maclay, 1995; 

Rose, 1972) (figs. 1 and 3). A high degree of spatial variability of aquifer properties exists within each depositional 

province. Within the San Marcos Platform, the regional dense 

member (RDM) (fig. 3) at the base of the Person Formation 

of the Edwards Group is a unit between the Person Formation and the Kainer Formation (Rose, 1972) characterized by 

very small permeability relative to other units of the Edwards 

aquifer. In some areas, the RDM may separate overlying freshwater from the underlying brackish water (Groschen, 1994). 

The depth of the RDM was used to subdivide the single layer 

MODFLOW models into two model layers (Lindgren and 

others, 2011b). In the western part of the Maverick Basin, the 

middle of the low-permeability McKnight Formation (fig. 3) 

was determined to be a unit that is similar to the RDM within 

that depositional province. Stratigraphic determinations of the 

top altitude and thicknesses of the RDM and midpoint of the 

McKnight Formation were made across the freshwater extent 

of the Edwards aquifer to map a thin unit that is conceptually 

referred to as the middle part of the Edwards aquifer (fig. 3) 

in this report. Within the Devils River Trend, a unit similar to 

the middle of the McKnight Formation and the RDM does not 

exist; thus, the midway point of the Devils River Limestone 

was used as a marker for the top altitude of the middle part of 

the Edwards aquifer in this area. The middle part then divides 

the aquifer into an upper part of the Edwards aquifer and a 

lower part of the Edwards aquifer. The upper part is representative of the Salmon Peak Formation and upper McKnight 

Formation in the Maverick Basin, the top of the Devils River 

Limestone in the Devils River Trend, and the Georgetown and 

Person Formations in the San Marcos Platform. The lower part 

of the Edwards aquifer is representative of the bottom of the 

McKnight Formation and the West Nueces Formation in the 

Maverick Basin, the bottom of the Devils River Limestone in 

the Devils River Trend, and the Kainer Formation in the San 

Marcos Platform (fig. 3). The Edwards aquifer was informally 

divided and named into three parts to aid in model conceptualization, implementation, parameterization, and discussion, and 

is not intended for use as naming conventions for any other 

purposes outside of this report.

Previously published information on top and bottom 

altitudes of the Edwards aquifer freshwater extent, as well as 

the extensive faults, is summarized by Lindgren and others 

(2004). The top and bottom altitudes of the Edwards aquifer model (Lindgren and others, 2004) were also used in the 

diffuse-flow model of Lindgren (2006), the specific model 

updated for this report. In addition to the top and bottom 

altitudes of the model extent of Lindgren and others (2004) 

for the model described in this report, interpretations of the 

tops and bases of hydrostratigraphic units were made based on 

cross-sectional information presented in Rose (1972), which 

was supplemented by data from geophysical logs (Railroad 

Commission of Texas, 2012). Interpretations of the top and 

base of the middle part of the Edwards aquifer were also made 

from the same downdip cross-sections and geophysical logs 

to extend this vertical aquifer subdivision (fig. 3) downdip as 

well. Geographic information system (GIS)-generated raster 

grids of top and bottom altitude information were interpolated 

for the revised downdip area using a natural neighbor interpolation method in ArcMap (Esri, 2015b). Updip of the 10,000-

mg/L dissolved-solids contour (A.L. Schultz, private consultant, written commun., 2000), the top and bottom altitudes and 

faults are the same as those of Lindgren and others (2004) and 

Lindgren (2006). Downdip of the 10,000-mg/L dissolved-solids contour, the raster grids were used to populate the updated 

extent of the model grid.

Review of Groundwater-Flow System

Groundwater recharge to the Edwards aquifer is from 

rainfall infiltration and stream losses within the recharge zone 

of the aquifer (which is the primary source of recharge), and 

to a lesser extent, groundwater flow (flux) from the underlying 

Trinity aquifer (figs. 2–3) (Lindgren and others, 2004). From 

the recharge zone, groundwater moves in directions controlled 

mainly by regional and local faulting (Maclay and Land, 

1988), but the general trend is for the water to move downdip into the confined part of the aquifer (fig. 2). In general, 

groundwater moves through fractured, highly transmissive 

rocks and is deflected by faults towards the east and northeast to discharge at many springs (fig. 1), of which Comal 

Springs and San Marcos Springs are the largest (Lindgren 

and others, 2004). Groundwater recharge from the Trinity 

aquifer is derived from the underlying Glen Rose Limestone 

(figs. 2–3), which is the uppermost unit of the Trinity aquifer 

(Lindgren and others, 2004). Historically, naturally occurring 

springs throughout the region were the primary locations of 

groundwater discharge. As groundwater withdrawals by wells 

increased throughout the freshwater extent of the aquifer to 

meet the needs of a growing population, industry, and agriculture (Lindgren and others, 2004), groundwater began discharging at both springs and wells. 

Water-Quality Zones

The freshwater zone of the Edwards aquifer is delineated 

to the south by a zone of brackish water (transition zone) 

where the aquifer transitions from fresh to saline water 

(fig. 4). The saline zone is located downdip of the transition zone. The transition zone minimizes the downdip flow 

of freshwater from the Edwards aquifer. Typically it has 

been defined on maps as the zone exceeding the 1,000-mg/L 

contour of equal dissolved-solids concentration (hereinafter 

referred to as the “1,000-mg/L dissolved-solids contour”). 

The 1,000-mg/L dissolved-solids contour is commonly 

referred to as the “bad-water line,” which coincides with the 

updip boundary of the transition zone (Schultz, 1993, 1994). 
8 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

MAVERICK

BASIN1

DEVILS RIVER

TREND1

SAN MARCOS

PLATFORM1

HYDROGEOLOGIC

UNITS

Anacacho

Limestone

Austin Chalk

Eagle Ford

Group

Buda

Limestone

Upper

confining

unit

Del Rio Clay

UPPER CRETACEOUS

Salmon

Peak

Formation2

Devils

River

Limestone2

Georgetown

Formation

LOWER CRETACEOUS

Cyclic and

marine members

(undivided)

Leached

member

Collapsed

member

Regional dense member

Grainstone

member

Dolomitic

member

Basal nodular member

Kirschberg

evaporite

member

Person Formation3 Kainer Formation3

McKnight

Edwards aquifer

Formation2

West Nueces

Formation2

Trinity

aquifer

Upper member

of the Glen

 Rose Limestone

Lower member

of the Glen

 Rose Limestone

Aquifer subdivision in the San Marcos Platform area4

I

II

III

IV

V

VI

VII

VIII

Erosional hiatus 

Very small

Very small

Anacacho

Limestone

Very small

Anacacho

Limestone

Very small

Moderate

Austin Chalk Moderate

Austin Chalk Moderate

Very small

Eagle Ford

Group

Very small

Eagle Ford

Group

Very small

Small

Buda

Limestone Small

Buda

Limestone Small

Very small

Del Rio Clay

Very small

Del Rio Clay

Very small

Very small

Very small

Large Small to moderate

Large

Moderate to large

Moderate to large

Very small

Moderate to large

Moderate

Large

Moderate

Moderate

Very small

Very small

Moderate

Small

Small

STRATIGRAPHIC UNITS

Descriptors "very small, small, moderate, and large" refer

 to relative permeability of stratigraphic units Edwards Group

SYSTEM

Upper

zone

Middle

zone

DEPOSITIONAL PROVINCE

1 Location shown on figure 1.

2

 Lozo and Smith (1964).

3 Modified from Rose (1972).

4

 Maclay and Small (1984).

Glen Rose

Limestone

Trinity Group

Glen Rose

Limestone

Glen Rose

Limestone

MODEL LAYERS

1

2

3

4

5

6

7

8

Modified from Maclay (1995)

Stratigraphic unit not included in the model

Edwards aquifer—Dotted where less certain

 Upper part

 Middle part

 Lower part

EXPLANATION

Figure 3. Correlation of Cretaceous stratigraphic units and hydrogeologic units, relative permeability and model layers in the 

Edwards aquifer model area, San Antonio region, Texas (modified from Maclay, 1995, fig. 11).
Updated Conceptualization and Numerical Model of the Edwards Aquifer 9

Previous models of the Edwards aquifer in the San Antonio 

region (Klemt and others, 1979; Maclay and Land, 1988; 

Thorkildsen and McElhaney, 1992) have used the 1,000-mg/L 

dissolved-solids contour as the southern boundary of the 

groundwater-flow model area. Other models have used the 

more conservative 10,000-mg/L dissolved-solids contour as 

an estimate of the southern boundary of the groundwater-flow 

model area (Lindgren and others, 2004; Lindgren, 2006, Lindgren and others, 2011b). The simulation for this study focused 

on the transition zone. Relocation of the southern boundary to 

minimize boundary effects on model results is described in the 

“Hydrologic Boundary Conditions” section of this report.

The thickness of the transition zone is highly variable. 

The contours of dissolved-solids concentrations (Schultz, 

1992, 1993, and 1994) are two-dimensional plan view 

contours, whereas the actual dissolved-solids concentration 

distribution is complex and three-dimensional (Lambert and 

others, 2010; Schultz, 1992, 1993, and 1994; Thomas and others, 2012). Furthermore, the water-quality data and geophysical log estimates of Schultz (1992, 1993, and 1994) were 

based on spatial information recorded during several decades. 

Therefore, considerable uncertainty exists as to the spatial 

and temporal distribution of dissolved-solids throughout the 

aquifer. Additionally, a disparity exists between the amount of 

data available for characterizing the freshwater, transition, and 

saline zones— a large amount of data is available for characterizing the freshwater zone, whereas few data are available 

for characterizing the transition and saline zones. Farther 

downdip from the “bad-water line,” data on heads, concentrations, and aquifer properties becomes increasingly scarce. This 

data deficit requires assumptions to be made based on expert 

knowledge of the system and its properties.

Dolomite dissolution within the karstic Edwards aquifer 

has created higher-permeability and higher-porosity preferential flow areas on the freshwater side of the “bad-water line” 

(Deike, 1990; Groschen, 1994; Hovorka and others, 1996; 

Hovorka and others, 1998). Just south of the “bad-water line,” 

porosity enhancement through dolomite dissolution may result 

from close proximity to the fresh-water part of the aquifer 

(Hovorka and others, 1998). As stated in Hovorka and others 

(1998), “enhancement of matrix permeability is seen regionally near the fresh-saline interface.” From Maclay and Small 

(1984, page 48), “a highly permeable belt of rocks exists along 

segments of the ‘bad-water’ line in areas where mixing of 

ground water of two different chemical types may increase the 

solution capacity of the water.” Dolomite dissolution may be 

enhanced where groundwater flow rates are faster; thus, wellconnected pores may have been further enlarged (Hovorka and 

others, 1998). Permeability within the saline zone is strongly 

related to porosity, and porosity estimates from the saline zone 

yield a lower average porosity than the freshwater zone (Hovorka and others, 1998). The updip edge of the transition zone 

has the highest estimates of matrix porosity and permeability. 

Structural controls, such as faults, play a major role in karst 

development because faults can focus dissolution and thereby 

contribute to solution enlargement in these locations (Hovorka 

and others, 1998). Some groundwater may potentially flow 

vertically (fig. 2) along faults (Groschen, 1994), but primarily faults control the horizontal direction of groundwater flow. 

Groundwater flow within the freshwater and transition zones is 

parallel to faults (Groschen, 1994; Lambert and others, 2010).

Aquifer Hydraulic and Transport Properties

Groundwater heads and flow are controlled by hydraulic 

conductivity (or the analogous transmissivity) and storativity 

in the aquifer, whereas effective porosity is the primary control 

on dissolved-solids transport. Hydraulic conductivity is known 

to span eight orders of magnitude in the Edwards aquifer 

(Hovorka and others, 1998). Calibrated hydraulic conductivity 

values for the diffuse-flow model ranged from 3 to 50,000 feet 

per day (ft/day) with an average of approximately 1,500 ft/day 

(Lindgren, 2006). Vertical variations of permeability exist in 

the Edwards aquifer (fig. 3) and, therefore, hydraulic conductivity exhibits vertical variability as well because it is related 

to permeability. Permeability of the saline zone is smaller than 

permeability in the freshwater zone; therefore, movement in 

the saline zone is thought to be less than in the freshwater zone 

(Groschen, 1994; Lambert and others, 2010). 

Storage within the Edwards aquifer is defined by specific 

yield in the unconfined part of the aquifer (outcrop) and specific storage in the confined part of the aquifer (subcrop where 

heads are above the top of the aquifer). For this assessment, 

the aquifer was conceptualized as confined and, therefore, specific storage—the volume of water that a unit volume of aquifer releases from storage under a unit decline in hydraulic head 

(Freeze and Cherry, 1979)—was used for the outcrop as well 

as the subcrop. This modeling assumption was made because 

the outcrop was not the focus area of the model. Specific 

storage values for the calibrated, confined part of the previous models (Lindgren and others, 2004; Lindgren, 2006) were 

on the order of 10-7 to 10-6 1/ft (per foot). Specific-storage 

parameters in the recharge zone were expected to be higher to 

account for specific yield, which typically is much larger than 

specific storage but was not simulated. Storage estimates are 

not well known in the transition and saline zones, but values 

used to calibrate the model were on the lower end of the range 

(Lindgren and others, 2004).

Effective porosity is a key property controlling transport of dissolved solids. Effective porosity is less than bulk 

porosity as it is a measure of the interconnectedness of the 

void space (Fetter, 1999). Effective porosity is not known on 

a regional scale and the spatial variability of effective porosity 

is not well understood. In a modeling context for contaminant 

transport problems, effective porosity serves as a final “tuning” parameter for matching concentration data (Zheng and 

Bennett, 2002), whereas hydraulic conductivity and storativity 

control groundwater heads and flows. Contaminant velocity 

is inversely proportional to effective porosity, which means 

velocities increase as effective porosity decreases. Effective porosity estimates are not available regionally for the 

Edwards aquifer because the aquifer is a complex karst system 
10 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

UNITED STATES

MEXICO

1,000 mg/L

10,000 mg/L

1,000 mg/L

10,000 mg/L

1

2

3 4

5

6

7

8

9

10

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras Springs

Leona Springs

San Marcos Springs

Barton Springs 

q

q

q

q

q

q

q

q

0 10 20 30 KILOMETERS

0 10 20 30 MILES

29°

30°

31°

98°

99°

100°

Base from U.S. Geological Survey and 

City of Austin digital data, 2010

Albers equal area projection, 

Texas Centric Mapping System,

North American Datum of 1983

TEXAS

COAHUILA

LLANO COUNTY

BURNET

COUNTY

BLANCO

COUNTY

WILLIAMSON

COUNTY

BASTROP COUNTY

TRAVIS

COUNTY

HAYS COUNTY

COMAL

COUNTY

KENDALL

COUNTY

GILLESPIE COUNTY

KERR

COUNTY 

EDWARDS COUNTY

REAL

COUNTY 

BANDERA

COUNTY

BEXAR

COUNTY

GUADALUPE

COUNTY

CALDWELL

COUNTY

GONZALES

COUNTY

WILSON

COUNTY

MEDINA

COUNTY

KINNEY UVALDE COUNTY

COUNTY 

MAVERICK

COUNTY 

ZAVALA COUNTY 

FRIO

COUNTY

ATASCOSA

COUNTY

LA SALLE COUNTY

DIMMIT COUNTY 

Lake

 Travis

Blanco River

Canyon 

Lake

Creek

Medina

Lake

Pedernales River

San

Marcos

River

Guadalupe

River

Nueces

River

Frio

River

Blanco

Creek

Dry

Frio River

West

Nueces

River

River

SanAntonio River

Cibolo

Creek

Leona

River

RIO BRAVO

Guadalupe River

Sabinal River

Dry Comal

Seco Creek

Hondo Creek

Colorado River

River Frio

Verd e Cr

Medina

Onion Creek

Barton Creek

Williamson

Creek

Bear Creek

Slaughter

 Creek

RIO GRANDE

NUECES RIVER BASIN

SAN ANTONIO RIVER BASIN

GUADALUPE RIVER BASIN

EXPLANATION

Regional-scale model area

 Catchment area

 Recharge zone (outcrop) of the Edwards 

 aquifer (modified from Puente, 1978)

 Freshwater zone of the Edwards aquifer

 Brackish-water (freshwater/saline-water)

 transition zone of the Edwards aquifer

 Saline zone q

1,000 mg/L

Model boundary

Line of equal dissolved solids concentration, in milligrams

 per liter (mg/L) (A.L. Schultz, written commun., 2000)

Contributing area—Basin extent shown terminates at 

 downdip extent of Edwards aquifer recharge zone

 Major river basin boundary

 Stream basin or contributing area boundary

Groundwater divide

Spring—Tail points in direction of flow

Contributing area

 Nueces-West Nueces River Basin

 Frio-Dry Frio River Basin

 Sabinal River Basin

 Area between Sabinal and Medina River Basins

 (includes Hondo and Seco Creeks)

 Medina River Basin

 Area between Medina River Basin and Cibolo-Dry Comal Creek Basins

 Cibolo and Dry Comal Creek Basins

 Guadalupe River Basin

 Blanco River Basin

 Barton Springs segment contributing area (includes Barton,

 Williamson, Slaughter, Bear, and Onion Creeks) 

1

2

3

4

5

6

7

8

9

10

Figure 4. Recharge, freshwater, brackish-water transition, and saline zones of the Edwards aquifer, San Antonio region, Texas 

(modified from Lindgren and others (2011a, fig. 2).
Updated Conceptualization and Numerical Model of the Edwards Aquifer 11

characterized by smaller matrix porosity as well as larger secondary porosity. Effective porosity is generally estimated on 

local scales from tracer tests, which may not be applicable on 

a regional scale. Bulk porosity estimates have been reported 

but how well these estimates reflect the interconnectedness of 

effective porosity is not well understood. An average of 0.18 

was given for porosity by Hovorka and others (1996), based 

on plug porosity data and geophysical logs, with estimates 

ranging from low values of 0.04 to 0.12 to high values of 

0.20 to 0.42 (Hovorka and others, 1996; Hovorka and others, 1993) with estimates up to 0.50 from cores (Hovorka and 

others, 1996). The matrix of the saline zone has pores that 

are less connected than those in the freshwater zone (Maclay 

and Small, 1984). Effective porosity estimates from previous modeling efforts range on the low end of reported values. 

Kuniansky and others (2001) used 0.015 to 0.035 for effective 

porosity to match estimated travel times derived from geochemical mixing models; the largest value for effective porosity used by Lindgren and others (2011a) for particle-tracking 

analyses was 0.01255. A uniform average effective porosity 

of 0.18 was used in both the conduit-flow model (Lindgren 

and others, 2004) and diffuse-flow model (Lindgren, 2006) to 

determine backward-tracked particle travel times for comparison to geochemical age-tracer results (Musgrove and others, 

2010). Results from the models indicated older groundwater 

ages than the geochemical results. For this assessment, a range 

of 0.025 to 0.50 (2.5 percent to 50 percent) was set to allow 

maximum flexibility in the history-matching effort. 

Model Development

The diffuse-flow model (Lindgren, 2006) included the 

San Antonio and Barton Springs segments of the Edwards 

aquifer (fig. 1). The San Antonio segment of the model was the 

focus area for history matching because this was the area of 

interest for this assessment. The San Antonio segment extends 

from near Brackettville, Tex. in Kinney County to the groundwater divide near Kyle, Tex. in Hays County, where the Barton 

Springs segment begins and extends to the Colorado River in 

Travis County. Both segments were simulated in the model 

because the groundwater divide near Kyle (fig. 1) is poorly 

defined and varies temporally (LBG-Guyton Associates, 1995) 

and, therefore, cannot be statically defined. 

The modeling was completed in two distinct phases: 

(1) the history-matching phase and (2) the predictive phase. 

Development of the history-matched model, within the 

history-matching phase, included construction of the model 

grid (both horizontal discretization and vertical layering), 

implementation of model boundary conditions and stresses, 

development of initial conditions, and calibration through 

history matching. History matching involves adjusting model 

input parameters to acquire a sufficient fit of the modelsimulated equivalents to the historical real-world observation 

dataset. The observation dataset consists of raw and processed 

data collected during a distinct time that the model will be 

adjusted to “match.” Another expectation of history matching, 

secondary to ensuring the model can reproduce groundwater 

responses within some time period, is improving the reliability 

of model predictions. History matching is intricately linked 

with conceptualization and with the prediction of interest in 

two ways: (1) the parameter estimation process within the history-matching exercise may yield results that indicate changes 

needed for the model, such as input parameters, boundary 

conditions, and initial conditions; and (2) the prediction of 

interest guides the history-matching effort as the model and 

observation dataset need to be in agreement with and supportive of the prediction. 

History-Matching Phase

 For the history-matching phase, the major changes from 

the diffuse-flow model (Lindgren, 2006) that were incorporated into the history-matched model (1999–2009), are 

summarized: 

1. More active model cells (fig. 5) were used to simulate 

the transition zone and saline zone downdip of the freshwater extent (fig. 4).

2. The one-layer model was converted to an 8-layer model 

(fig. 3) to better represent vertical aquifer heterogeneity 

and anisotropy as well as simulate variable-density flow 

and transport. Discussion of the model is in terms of the 

“upper” (layers 1–3), “middle” (layer 4), and “lower” 

(layers 5–8) parts of the Edwards aquifer based on the 

hydrostratigraphy.

3. SEAWAT version 4 was used to explicitly simulate 

variable-density flow and transport.

4. The history-matching period was updated to 1999–2009.

5. Parameter estimation, through history matching, was 

executed using highly parameterized inversion with 

PEST to calibrate to heads, discharges at major springs, 

and dissolved-solids concentrations at wells.

Next, the predictive phase needs to be planned early 

during the history-matching phase because the model design 

and parameterization for the history matching need to meet the 

intended use of the model, that is, the predictions of interest. 

For example, if the prediction of interest is to estimate drawdown around a well field under hydrologically wet conditions, 

the history-matching must contain wet conditions in the period 

of record and the well-field must be appropriately simulated 

in the model. The predictions of interest in this assessment are 

(1) changes in dissolved-solids concentration at production 

wells near the transition-zone interface, (2) total spring discharge at Comal and San Macros Springs, and (3) head at well 

J-17 during drought conditions of 1950–56. The predictions 

are mentioned here for readers to bear in mind the intended 

use of the model during the discussion of the model development and history-matching effort.
12 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement 30° 29° 98° 99° 100° Base from U.S. Geological Survey digital data, 2012 Texas Centric Mapping System Albers equal area projection North American Datum of 1983 Contributing area—Basin extent shown terminates at downdip extent of Edwards aquifer recharge zone Major river basin boundary Stream basin or contributing area boundary Spring—Tail points in direction of flow EXPLANATION Boundary condition Specified-flux boundary (layer 8) No-flow, constant dissolved-solids concentration southern boundary Model boundary Active model area Recharge zone (outcrop) of the Edwards aquifer (modified from Puente, 1978) Edwards aquifer subcrop Inactive model area

q

KIMBLE COUNTY SUTTON COUNTY

TRAVIS COUNTY GILLESPIE COUNTY

BLANCO COUNTY

HAYS COUNTY

KERR COUNTY

KENDALL COUNTY

CALDWELL

COUNTY

COMAL COUNTY

GUADALUPE COUNTY

GONZALES COUNTY

BEXAR COUNTY

MEDINA COUNTY

UVALDE COUNTY

KINNEY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA COUNTY

KARNES COUNTY FRIO COUNTY ZAVALA COUNTY MAVERICK COUNTY

REAL COUNTY BANDERA COUNTY EDWARDS COUNTY

Medina

Lake

Canyon Lake

South Llano River

Dr Dy raw

P ni ot Creek

Dry Frio River

Long Kent Creek

Johnson Creek

Guadalupe River

Pedernales River

Miller Creek

Blanco River

Slaughter Creek

Colorado River

On oi n Creek

Cedar Creek

York Creek

San Marcos River

Plum Creek

Picosa Creek

l E C m er ke

Chaparrosa Creek

Turkey Creek

Palo Blanco Creek

Elm Creek

Nueces River

Blanco Creek

Sabinal River

Seco Creek

Leona River

Frio River

dno H Co keer

San Miguel Creek

Black Creek

Medina River

Atascosa River

Lucas Cree k

Hondo Creek

San Antonio Creek

Cibolo Creek

Ecleto Creek

Sandies Cree k

E

ml Creek

Elm Creek

San 

Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras Springs

Leona Springs

San Marcos Springs

Barton Springs 

q

q

q q

q

q

q

q

TEXAS

Area enlarged

Modelextent 0 10 20 KILOMETERS

0 10 20 MILES

Figure 5. Active model area and boundary conditions, San Antonio region, Texas.
Updated Conceptualization and Numerical Model of the Edwards Aquifer 13

Discretization

The horizontal grid discretization from the diffuse-flow 

model (Lindgren, 2006) was retained in the SEAWAT model. 

The uniformly spaced finite-difference grid has 370 rows and 

700 columns. The dimensions of the grid cells are uniformly 

0.25 mile (1,320 feet) along rows and columns. The grid was 

oriented to align with the regional direction of groundwater 

flow (Lindgren and others, 2004) and the major faults (fig. 6) 

in the Balcones fault zone, with preference given to the direction of groundwater flow and orientation of major faults near 

Comal and San Marcos Springs. Accordingly, the grid was 

rotated 35 degrees counterclockwise from horizontal.

The grid was refined to include eight model layers to better represent vertical heterogeneity and groundwater flow and 

because SEAWAT requires a greater level of vertical discretization than MODFLOW. Also, the aquifer was divided based 

on the hydrostratigraphic units defined for the Edwards aquifer 

(fig. 3) into the upper (layers 1-3), middle (layer 4), and lower 

(layers 5-8) parts of the Edwards aquifer as discussed in 

“Updated Conceptualization of Hydrogeology for Numerical 

Model” section. Subdividing the Edwards aquifer into upper, 

middle, and lower aquifer parts allows for better representation of vertical groundwater flow, which is controlled by 

the vertical hydraulic conductivity of the middle part of the 

Edwards aquifer. Conceptually, the middle part of the Edwards 

aquifer might impede, or allow, vertical flow between the 

upper and lower parts of the Edwards aquifer, but this vertical flow is dependent on any information in the observation 

dataset regarding vertical head differences. Because conceptualization of the middle part of the Edwards aquifer was 

designed to control vertical flows between the subdivisions, 

it was assumed that horizontal variations in specific storage 

and effective porosity could be ignored and represented with 

uniform aquifer properties. Following the approach of Lindgren and others (2011b), isotropic conditions were assumed 

(vertical hydraulic conductivity of each model cell was equal 

to the horizontal hydraulic conductivity). Multiple layers were 

used for the upper and lower parts of the Edwards aquifer 

(three layers and four layers, respectively [fig. 3]) to allow for 

better simulation of variable-density flow and transport (that 

is, discretize the transition zone). However, layers within each 

aquifer subdivision contained the same property values. 

The top (top of layer 1) and bottom (bottom of layer 8) of 

the model are consistent with previous models of the Edwards 

aquifer (Lindgren and others, 2004; Lindgren, 2006; Lindgren and others, 2011b) for the active model area inherited 

(obtained from) from these models. The bottom of the model 

was calculated by subtracting the thickness of the Edwards 

aquifer (fig. 6b) from the top of the aquifer (fig. 6a). The 

model in this assessment simulates a larger downdip extent 

of the aquifer than previous models; therefore, the top and 

bottom of the new active extent of the model is based on the 

interpolated top and bottom of the Edwards aquifer as discussed in the “Updated Conceptualization of Hydrogeology 

for Numerical Model” section. Stratigraphic determinations 

of the top of layer 4 were made and interpolated (fig. 7a). The 

thickness of the upper part of the Edwards aquifer was calculated by subtracting the top of layer 4 from the top of layer 

1 (fig. 6a). The upper part of the Edwards aquifer was then 

subdivided into three layers of equal thickness. The thickness 

of layer 4 (fig. 7b), as based on stratigraphic determinations, 

was subtracted from the top of layer 4 to determine the top of 

layer 5. The resulting lower part of the Edwards aquifer (top of 

layer 5 to the bottom of layer 8) was subdivided into 4 equal 

parts to form layers 5 through 8. 

As the Edwards aquifer is a sloping aquifer with an 

outcrop recharge zone at the surface corresponding to layers 5 

through 8 which lie underneath a confined subcrop part of the 

aquifer farther downdip (layers 1 through 4), the active extents 

of each layer are different depending upon the water saturation 

of each cell and if the cell lies within the subcrop or outcrop. 

Model layer 8 has the largest number of active cells, as this 

layer includes the bottom of the Edwards aquifer. Shallower 

depths in the aquifer have smaller number of active cells for 

each layer. Model layer 1 has the least number of active model 

cells. The number of active layers for each row and column 

in the model (fig. 8) was determined based on the active cells 

in the diffuse-flow model (Lindgren and others, 2004) and 

layer bottoms and was modified where needed to aid model 

convergence. A cell was made inactive if the head fell below 

the bottom of the cell. Other inactive areas of the model grid 

are inactive because the Edwards aquifer does not exist at that 

location (north of the recharge zone) or it is not within the area 

of interest (deep in saline zone) (figs. 4 and 8).

The model was temporally discretized to have monthly 

stress periods from January 1999 through December 2009, for 

a total of 132 stress periods. Therefore, the model simulates 

monthly changes in recharge, Trinity aquifer flux, groundwater 

withdrawals, and spring discharges.

Hydrologic Boundary Conditions

Boundary conditions for the model are set along the 

lateral model boundaries of the active model area (fig. 5) and 

are specified where water is entering or leaving the aquifer. 

Boundary conditions include recharge from precipitation in 

the recharge zone by the Recharge Package (RCH), flux from 

the underlying Trinity aquifer by the Well Package (WEL), 

discharge from groundwater withdrawals by the WEL and 

discharge to springs by the Drain Package (DRN) (Harbaugh, 

2005). For SEAWAT, dissolved-solids concentrations for water 

are required to be specified for cells in which water is entering the aquifer (RCH or WEL) and are assumed to be equal to 

the concentration in the cell for those cells where groundwater 

leaves the aquifer.

Lateral Model Boundaries 

The western, northern, and eastern lateral extents of the 

active part of the model are the same as for the diffuse-flow 

model (Lindgren, 2006), whereas the southern extent has been 
14 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

30°

29°

100° 99° 98°

30°

29°

100° 99° 98°

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

B. Thickness of Edwards aquifer active model area

A. Faults and altitude of top of Edwards aquifer active model area 

900

850

600

850

950

800

900

750

650

550

1,000

750

1,100

800

500

450

700

–500

500

–2,500 –3,000

–1,000 –3,500

1,000

–4,000

0

500

–500

–1,500

0

–2,000

500

500

Thickness of Edwards aquifer

 active model area, in feet

 24

 700

 1,228

Model extent

Recharge zone boundary

 (modified from Puente, 1978)

Line of equal thickness—

 Intervals 50 and 100 feet. 

 Modified from Hovorka and

 others, 1996

Active model-area boundary

Spring—Tail points in direction

 of flow

EXPLANATION

q

 –5,507

 –2,000

 1,408

Model extent

Recharge zone boundary

 (modified from Puente, 1978)

Fault (Lindgren and others,

 2004)

Structure contour—Shows 

 altitude of top of Edwards 

 aquifer. Interval 500 feet. 

 Datum is National Geodetic

 Vertical Datum of 1929.

 Modified from Collins

 and Hovorka, 1997

 1,000

 500

 0

 –500

 –1,000

 –1,500

 –2,000

 –2,500

 –3,000

 –3,500

 –4,000

Active model-area boundary

Spring—Tail points in direction

 of flow

EXPLANATION

Altitude of top of Edwards

 aquifer active model area,

 in feet—Datum is National

 Geodetic Vertical Datum of

 1929

q

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras

Springs

Leona Springs

Barton Springs 

San Marcos

Springs

q

q

q

q

q

q

q

q

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras

Springs

Leona Springs

Barton Springs 

San Marcos

Springs

q

q

q

q

q

q

q

q

Figure 6. A, Faults and altitude of top of Edwards aquifer active model area, and B, thickness of Edwards aquifer active model 

area, San Antonio region, Texas.
Updated Conceptualization and Numerical Model of the Edwards Aquifer 15

Figure 7. A, Altitude of top, and B, thickness of the middle part of the Edwards aquifer (model layer 4), San Antonio region, 

Texas.

30°

29°

100° 99° 98°

30°

29°

100° 99° 98°

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

A. Altitude of top of middle part of Edwards aquifer

B. Thickness of middle part of Edwards aquifer

Altitude of top of middle

 part (model layer 4) of 

 Edwards aquifer, in feet—

 Datum is National 

 Geodetic Vertical 

 Datum of 1929

 –5,160.1 to –3,612.6

 –3,612.5 to –2,615.7

 –2,615.6 to –1,778.4 

 –1,778.3 to –1,024.9 

 –1,024.8 to –321.6

 –321.5 to 311.0

 311.1 to 1,233.0

Model extent

Active model-area boundary

EXPLANATION

Thickness of middle

 part (model layer 4)

 of Edwards aquifer,

 in feet

 12.9 to 21.0

 21.1 to 26.0 

 26.1 to 31.0

 31.1 to 44.6

Model extent

Active model-area boundary

EXPLANATION
16 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement Figure 8. Number of active model layers per cell, San Antonio region, Texas. 30° 29° 98° 99° 100° Base from U.S. Geological Survey digital data Texas Centric Mapping System Albers equal area projection North American Datum of 1983

0 10 20 KILOMETERS

0 10 20 MILES

EXPLANATION

Inactive model areaModel boundary

 7 8

 5 6

 3 4

Active model area, by number of active 

 1 2

model layers per

 cell—For each location, 1–8 model layers are active

KIMBLE COUNTY SUTTON COUNTY

TRAVIS COUNTY GILLESPIE COUNTY

BLANCO COUNTY

HAYS COUNTY

KERR COUNTY

KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL COUNTY

BANDERA COUNTY

GUADALUPE COUNTY

GONZALES COUNTY

BEXAR COUNTY

MEDINA COUNTY

UVALDE COUNTY

KINNEY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA COUNTY

KARNES COUNTY FRIO COUNTY ZAVALA COUNTY MAVERICK COUNTY

EDWARDS COUNTY
Updated Conceptualization and Numerical Model of the Edwards Aquifer 17

extended downdip. The northern extent of the model corresponds to the northern extent of the recharge zone (fig. 5). 

The northern part of the eastern extent of the active part of 

the model corresponds with the Colorado River (fig. 5). The 

western and southern extents of the active part of the model 

are less well-defined hydrogeologic boundaries. The western 

extent of the active part of the model coincides with a poorly 

defined groundwater divide (fig. 1) near Brackettville in Kinney County (Garza, 1966; LBG-Guyton Associates, 1995)—

minimal flow across this boundary is generally assumed. Thus, 

the southern part of the western model boundary was simulated as a no-flow boundary similar to the diffuse-flow model 

(Lindgren, 2006). Last, the northern part of the western extent 

of the active part of the model was simulated as a specifiedflux boundary using WEL which better simulates heads near 

the western boundary. This boundary condition was also 

similar to the boundary condition for the diffuse-flow model 

(Lindgren, 2006). 

The 10,000-mg/L dissolved solids concentration contour 

(A.L. Schultz, private consultant, written commun., 2000) 

has been used several times in previous MODFLOW models 

as a no-flow boundary (Lindgren and others, 2004; Lindgren, 

2006; Lindgren and others, 2011b). Because the objective 

of this modeling effort was to simulate the transition zone, 

the southern boundary (no-flow, constant dissolved-solids 

concentration) needed to be moved to minimize boundary 

effects. The initial goal was to set the southern boundary at the 

estimated 35,000-mg/L dissolved-solids concentration contour, 

but insufficient data prevented a full characterization of this 

concentration contour on a regional scale. Therefore, the 

boundary was based on an approximate location of the 

35,000-mg/L contour (fig. 9) where data were available 

(described in the “Initial Dissolved-Solids Concentrations” 

section of the “History-Matching Phase” section), and in areas 

where data were scarce, the boundary was set sufficiently far 

downdip from the 10,000-mg/L contour. This boundary was 

based on the assumption that any groundwater flow moves 

parallel to this boundary, and flow across this boundary is 

assumed to be minimal. 

The eastern extent of the active part of the model is 

consistent with the extent of Lindgren (2006). However, the 

northern part of the eastern boundary located at the Colorado 

River is simulated as a no-flow boundary condition in this 

assessment (fig. 5). The Colorado River serves as a regional 

sink (groundwater discharge), but is a small component of the 

groundwater budget far away from the predictive area of interest. Therefore, any effects of this boundary on results in the 

predictive area of interest should be minimal. Spring discharge 

is simulated within the model at Barton Springs near the eastern boundary; Barton Spring discharges to the Colorado River 

(fig. 5).

Groundwater Recharge

Groundwater recharge is specified where water enters the 

Edwards aquifer. The two sources of recharge to the Edwards 

aquifer are stream losses and rainfall infiltration in the 

recharge zone and, to a lesser extent, flux from the underlying Trinity aquifer. Concentrations of dissolved solids for flux 

from the RCH and for the underlying Trinity aquifer were set 

to 325 mg/L, which is assumed to be background dissolvedsolids concentration for this assessment and is near the median 

of 313 mg/L dissolved-solids concentration determined by 

Musgrove and others (2010) from groundwater sampled in the 

San Antonio segment of the Edwards aquifer.

Recharge to the Edwards aquifer occurs primarily by 

seepage from streams to the aquifer as these streams cross the 

recharge zone (fig. 5) (Lindgren and others, 2004; Maclay, 

1995). Additional recharge is from rainfall infiltration in the 

interstream areas of the recharge zone. A Neumann boundary 

condition (Franke and others, 1987), which is also referred 

to as a specified-flow (Anderson and Woessner, 1992) or 

specified-flux (Winston, 2015) boundary condition, was 

implemented using the RCH (Harbaugh, 2005) to simulate 

groundwater recharge. Recharge was assigned in the same 

manner as was done in Lindgren and others (2004). Monthly 

recharge rates were applied for the transient simulation period 

(1999–2009). The recharge zone corresponds to cells in model 

layers 5 through 8 in the outcrop area (figs. 3 and 5), and, 

therefore, recharge was applied into the topmost active model 

layer for each cell in the RCH. 

A Neumann boundary condition was implemented for 

the northern model boundary to account for inflow from the 

underlying Trinity aquifer (figs. 2 and 5). The boundary 

condition was implemented using the WEL (Harbaugh, 2005) 

similar to what was used in previous models (Lindgren and 

others, 2004; Lindgren, 2006; Lindgren and others, 2011b). 

The assumption was made that water-level fluctuations in 

the Trinity aquifer near the northern boundary are similar 

to water-level fluctuations in the Edwards aquifer, resulting 

in a constant hydraulic gradient and rate of flow across the 

boundary, which is consistent with assumptions in previous 

models (Lindgren and others, 2004; Lindgren, 2006, Lindgren 

and others, 2011b). For this report, the constant fluxes at this 

boundary from the one-layer diffuse-flow model (Lindgren, 

2006) were applied to the corresponding model cell locations 

(row, column) for model layer 8 (bottom model layer) for all 

stress periods.

Groundwater Discharge

Groundwater discharge is specified for the model or calculated by the model where water leaves the Edwards aquifer. 

Groundwater primarily discharges from the Edwards aquifer 

through withdrawals at wells and as discharge at springs. A 

small amount of groundwater also discharges to the Colorado 

River. In SEAWAT, dissolved-solids concentrations at groundwater-discharge locations are equal to the dissolved-solids 

concentration of the cell in which each particular boundary 

condition lies. Therefore, if the dissolved-solids concentration 

at a given cell is 325 mg/L, water withdrawn from a well or 
18 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

Figure 9. Lines of equal dissolved-solids concentrations (A.L. Schultz, private consultant, written commun., 2000), dissolvedsolids concentration locations of measured or estimated values based on data collected from 1930 through 2009, and distribution 

of interpolated dissolved-solids concentrations, San Antonio region, Texas.

0 10 20 30 KILOMETERS

0 10 20 30 MILES

29°

30°

31°

98°

99°

100°

Base from U.S. Geological Survey digital data

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

Interpolated dissolved-solids concentration,

 in milligrams per liter (mg/L)

 Less than or equal to 1,000

 1,001 to 3,000

 3,001 to 10,000

 10,001 to 35,000

 Greater than 35,000

Recharge zone boundary (modified from Puente, 1978)

Line of equal dissolved-solids concentration,

 in mg/L (A.L. Schultz, written commun., 2000)

Model boundary

Active model-area boundary

Spring—Tail points in direction of flow

Dissolved-solids concentration, in mg/L

 Less than or equal to 1,000

 1,001 to 3,000

 3,001 to 10,000

 10,001 to 35,000

 Greater than 35,000

EXPLANATION

q

10,000 mg/L

3,000 mg/L

1,000 mg/L

UNITED STATES

MEXICO

COAHUILA

TEXAS

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras Springs

Leona Springs

San Marcos Springs

Barton Springs 

q

q

q

q

q

q

q

q

LLANO COUNTY

BURNET

COUNTY

BLANCO

COUNTY

WILLIAMSON COUNTY

BASTROP

COUNTY

TRAVIS

COUNTY

HAYS

COUNTY

COMAL

COUNTY

KENDALL

COUNTY

GILLESPIE COUNTY

KERR

COUNTY 

EDWARDS COUNTY

REAL

COUNTY 

BANDERA

COUNTY

BEXAR

COUNTY

GUADALUPE

COUNTY

CALDWELL

COUNTY

GONZALES

COUNTY

WILSON

COUNTY

MEDINA

COUNTY

KINNEY UVALDE COUNTY

COUNTY 

MAVERICK

COUNTY 

ZAVALA COUNTY 

FRIO

COUNTY

ATASCOSA

COUNTY

LA SALLE COUNTY

DIMMIT COUNTY 
Updated Conceptualization and Numerical Model of the Edwards Aquifer 19

exiting the model as spring discharge in that cell is assigned a 

dissolved-solids concentration of 325 mg/L. 

Well Discharge

Groundwater withdrawals by wells were simulated as 

a Neumann boundary condition using the WEL (Harbaugh, 

2005) with monthly time steps. Groundwater-withdrawal 

data were compiled and distributed within the model grid 

to simulate the period from 1999 through 2009. Withdrawal 

data from 2000 through 2009 were provided by water-use 

category by the Edwards Aquifer Authority (EAA), SAWS, 

Bexar Metropolitan Water System (Bexar Met), and Fort Clark 

Municipal Utility District (MUD). The 1999 withdrawal data 

were obtained from the conduit-flow MODFLOW model of 

the Edwards aquifer (Lindgren and others, 2004).

Site-specific withdrawals were provided by several agencies and by the water-use categories that included municipal, 

industrial (includes manufacturing, mining, and power generation), irrigation, and livestock (Edwards Aquifer Authority, 

written commun., 2004, 2012; Bexar Metropolitan Water 

System, written commun., 2004, 2012; Fort Clark Municipal 

Utility District, written commun., 2004, 2012). Site-specific 

withdrawal data provided as annual totals were parsed to 

the appropriate stress period (months) by the same methods 

discussed in Lindgren and others (2011b). Municipal, industrial, irrigation, and livestock withdrawals represent monthly 

withdrawals except for some withdrawals from EAA. Withdrawal data from SAWS were provided by well field and distributed evenly per well based on the number of active wells 

(San Antonio Water System, 2012). Irrigation withdrawal data 

were also provided by the Texas Water Development Board 

(TWDB) by county (Texas Water Development Board, 2004) 

and were spatially distributed to a site-specific location for the 

land-use categories of row crops, small grains, and orchards/

vineyards.

The vertical assignment of withdrawals to a model layer 

was done on the basis of the percentage of the well open 

interval in each of the model layers and the transmissivity of 

the model layer. If a well was open in more than 50 percent 

of a layer, then this layer was included for that well, except 

for layer 4 (RDM) in which no water was withdrawn. Where 

well construction information was not available, the well was 

assumed to be open to the entire thickness of the Edwards 

aquifer. Well locations and average annual groundwater 

withdrawals during 1999–2009 are shown in figure 10, and 

ranged from less than 1 to more than 22,000 acre-feet per year 

(acre-ft/yr).

Spring Discharge

A Cauchy boundary condition (Anderson and Woessner, 

1992; Franke and others, 1987), also known as a head-dependent flux boundary condition, was implemented to simulate 

groundwater discharge from major springs (Las Moras, Leona, 

San Pedro, San Antonio, Comal, San Marcos, and Barton 

Springs; fig. 1) using the DRN (Harbaugh, 2005). The altitude 

of the spring orifice and the spring conductance in the DRN 

control the spring discharge from the model (Harbaugh, 

2005). Because Barton Springs and Las Moras Springs are 

near model boundaries, these springs were simulated to allow 

groundwater to discharge from the model, but spring discharge 

observation data were not used for the purpose of history 

matching. Data from Barton Springs were not used for history 

matching because the Barton Springs segment of the Edwards 

aquifer is not the focus of this model or previous models of the 

study area (Lindgren and others, 2004; Lindgren, 2006; Lindgren and others, 2011b). Data from Las Moras Springs were 

not used for history matching because much of its discharge 

might be derived from west of the western model boundary 

(Lindgren and others, 2004) (fig. 1). Hueco Springs (fig. 1), 

in Comal County north of Comal Springs, was not simulated 

because of the complexity of the flow relations between 

the Edwards and Trinity aquifers in the area of the springs 

and because much of the Hueco Springs discharge might be 

derived from the Trinity aquifer (Lindgren and others, 2004).

Initial Conditions

Transient models require initial conditions, also referred 

to as starting conditions, for the beginning of the model 

simulation. Previous MODFLOW models (Lindgren and 

others, 2004; Lindgren, 2006; Lindgren and others, 2011b) 

required initial heads, but because the history-matched model 

(1999–2009) documented in this report is a transient SEAWAT 

model, both initial heads and initial dissolved-solids concentrations at the start of 1999 are required. 

Initial Heads

A constant-density version of the model was initially 

run under steady-state conditions by using MODFLOW-2005 

(Harbaugh, 2005) to determine initial-head conditions for 

each of the eight layers in the model. This steady-state model 

was based on long-term (1978–1999) average conditions of 

groundwater recharge and discharge. Heads from the steadystate MODFLOW model were used as initial heads in the 

transient SEAWAT model. It is expected that matches to observation data during the first year of the history-matching period 

may be less desirable than matches to the following years as 

any errors in estimation of initial heads should be dampened 

out by the second year of the model simulation. 

Initial Dissolved-Solids Concentrations

The initial dissolved-solids concentrations (fig. 9) 

used in the model were spatially interpolated from existing 

concentration (before 1999) data including water-quality 

samples and geophysical logs from Schultz (1992, 1993 and 

1994). Because data were sparse in parts of the study area, 

two additional datasets were added: (1) groundwater withdrawal locations from the model (assuming a dissolved-solids 

background concentration of 325 mg/L) and (2) transect well 
20 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement Figure 10. Average annual groundwater withdrawals during 1999–2009, and the production wells of predictive interest near the brackish-water transition zone, San Antonio region, Texas. Base from U.S. Geological Survey digital data Albers equal area projection, Texas Centric Mapping System North American Datum of 1983

0 10 20 MILES0 10 20 KILOMETERS

EXPLANATION

Active model areaModel extentRecharge zone boundary (modified from Puente, 1978)Approximate line of 1,000 milligrams per liter 

dissolved-solids concentration (A.L. Schultz, 

 written commun., 2000)Production well of predictive interest* and well number

Average annual groundwater withdrawals during 1999–2009 by well, in acre-feet

22,812 to 13,007

 13,006 to 6,777 6,776 to 4,074 4,073 to 1,696

 1,695 to 611 610 to 139 138 to 0.000001

5

* Prediction of dissolved-solids concentration change during simulated 

7-year (1950–56) drought of record conditions made at these production wells

UVALDE COUNTY

EDWARDS COUNTY

BEXAR COUNTY

KERR COUNTY

MEDINA COUNTY

FRIO COUNTY

KINNEY

COUNTY

ZAVALA COUNTY

REAL COUNTY

HAYS COUNTY

WILSON COUNTY

ATASCOSA COUNTY

BANDERA COUNTY

KENDALL

COUNTY

KARNES COUNTY

COMAL COUNTY

GUADALUPE COUNTY

GILLESPIE COUNTY

TRAVIS COUNTY BLANCO COUNTY

MAVERICK

COUNTY

GONZALES

COUNTY

CALDWELL

COUNTY

98° 99° 100°

30°

29°

Area enlarged

29°30'29°20'

98°20' 98°30'

24

22

21

9

11

13

23

6

19

17

16

14

15

18

12

20

25

5 3

41278

10
Updated Conceptualization and Numerical Model of the Edwards Aquifer 21

dissolved-solids concentration estimates obtained from Lambert and others (2010).

Dissolved-solids concentration data from 1930 through 

2009 were obtained from the USGS National Water Information System (NWIS) (U.S. Geological Survey, 2013), the 

TWDB Groundwater Database (Texas Water Development 

Board, 2013), and the USGS National Produced Waters 

Geochemical Database (Blondes and others, 2013). Dissolvedsolids concentration data points for data collected from 1999 

through 2009 were used for the observation dataset (discussed 

in detail in the “Calibration through History Matching” section). Dissolved-solids concentration data before 1999 from 

these databases (including the work of Schultz) were used to 

form the initial concentrations dataset (fig. 9). If construction 

information was not available for the well or a sampling depth 

was not given, then the dissolved-solids concentration was not 

used.

The work of Schultz (1992, 1993, 1994; A.L. Schultz, 

private consultant, written commun., 2000) to characterize 

dissolved-solids concentrations was used for estimating the 

distribution of the transition zone within the Edwards aquifer. 

Schultz’s publications (1992, 1993, and 1994) provide maps 

of well locations, depths, and dissolved-solids concentrations. 

These data were used to create two-dimensional contours 

of equal dissolved-solids concentration regionally for the 

Edwards aquifer. The dissolved-solids data used to create the 

concentration contours in Schultz’s publications (1992, 1993, 

and 1994) were georeferenced by using Esri ArcMap 10.0 

(Esri, 2015a) and digitized to supplement the other dissolvedsolids data compiled. 

Few data were available to characterize dissolvedsolids concentrations downdip of the historical 10,000-mg/L 

dissolved-solids contour by Schultz (A.L. Schultz, private 

consultant, written commun., 2000). Many boreholes have 

been drilled in south-central Texas for oilfield exploration 

and are present downdip of the transition zone. The Railroad Commission of Texas maintains borehole geophysical 

logs collected from wells that intersect the Edwards aquifer 

that can provide valuable information about the water quality of the aquifer. Available geophysical logs were compiled, 

reviewed, and assessed to determine if they penetrated the 

Edwards aquifer and contained pertinent information to calculate dissolved solids. Within the study area, 120 geophysical 

logs contain the necessary information. Data from geophysical logs obtained from the Railroad Commission of Texas Oil 

and Gas Well Logs database (Railroad Commission of Texas, 

2012) were converted to estimates at depth of dissolved-solids 

concentrations.

Multiple techniques using borehole geophysical logs 

have been developed to calculate resistivity of the formation 

water (Rw) from saturated, clean (shale-free) water zones 

(Crain, 2013). One method used for this report and for previous studies of the transition zone, the Archie method, uses 

deep formation resistivity, porosity, and a cementation factor. 

Deep formation resistivity is the resistivity beyond the invaded 

zone, which is the volume close to the borehole wall where 

some or all of the moveable fluids have been displaced by mud 

filtrate (Schlumberger, 2015). A cementation factor (m) of 2 is 

a common value for carbonates and has been used successfully 

in the transition zone (Schultz, 1992) and was used for this 

report. The resistivity of the formation water using this method 

was calculated by the following equation (Evenick, 2008):

Rw = Rt øm (1)

where 

Rw is the resistivity of the formation water, in 

ohm-meters; 

Rt is the deep resistivity of the formation beyond 

the invaded zone, in ohm-meters;

ø is the porosity of the water saturated zone (the 

percent of the zone that is pores or voids), 

dimensionless, in decimal format; and

m is the cementation factor, dimensionless (a 

cementation factor of 2 was used for this 

report).

Another method to calculate Rw is based on the understanding of spontaneous potential (SP) log responses. When 

the drilling mud is fresher than the formation fluid, the shale 

baseline is typically the most positive SP reading and has negative SP deflections in clean (shale-free) zones. The static selfpotential (SSP) is the magnitude of deflections (in millivolts) 

from the shale baseline, which is a response of the difference 

between drilling mud resistivity and Rw. The SP method uses 

this relation to calculate Rw. For the SP method to work, the 

zone assessed must be clean (shale-free) with at least some 

permeability, and the drilling mud and formation fluid must 

have different resistivities. The SSP, log header information, 

and formation constants are used to solve multiple equations 

to calculate Rw for a specific zone. The SP method, which uses 

SP log responses, borehole temperature, and mud filtrate resistivity values to calculate Rw, is explained in Crain (2013). 

Once Rw was determined by using the Archie or SP 

method, the fluid conductivity was calculated using the following equation:

C = 10,000 / Rw (2)

where

C is the conductivity of the formation water, in 

microsiemens per meter; and 

Rw is the resistivity of the formation water, in 

ohm-meters.

The fluid conductivity was then temperature corrected to 

estimate the specific conductance at 25 degrees Celsius. During previous studies of the transition zone in the San Antonio 

segment of the Edwards aquifer, observations of specific conductance and dissolved-solids concentration were related by 

regression analysis (appendix 2 of Thomas and others, 2012). 

Based on the established correlation between dissolved-solids 

concentration and specific conductance in the transition zone 

(Lambert and others, 2009), specific conductance values were 
22 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

converted to dissolved-solids concentration estimates using 

the following equation:

EDS = 0.6522SC + 25.77 (3)

where

EDS is the estimated dissolved-solids 

concentration, in milligrams per liter; and 

SC is the specific conductance, in microsiemens 

per centimeter at 25 degrees Celsius.

Several factors can introduce errors into estimates of 

dissolved-solids concentrations calculated from geophysical 

logs and quality-control measures were used to assess possible 

errors. Where possible, the Archie and SP methods were used 

to estimate dissolved-solids concentrations, and an average 

of the two dissolved-solids concentrations was computed. 

Three-dimensional grids were created from dissolved-solids 

estimates and reviewed for outliers, which were removed from 

the dataset.

Dissolved-solids concentration data from the databases, 

reports (Schultz, 1992, 1993, 1994), and geophysical logs 

were combined into one dataset. Three-dimensional spatial 

interpolation was attempted to the 8-layer grid using threedimensional kriging, but lack of information on directional 

variograms (direction-dependent functions that quantify 

the spatial continuity or correlation of datasets) (Isaaks and 

Srivastava, 1989; Oliver and others, 2008) and data gaps made 

interpolation difficult and introduced many anomalies. Therefore, the dataset was interpolated in two-dimensions instead 

(the depth parameter was excluded) through kriging with 

PEST’s “ppk2fac” and “fac2real” utilities (Doherty, 2011). 

 A near singular matrix during the interpolation process 

required declustering of the dataset to remove similar values 

near one another. The interpolated dissolved-solids concentration distribution (fig. 9) was used as the initial concentrations 

for all model layers. 

Requirements for SEAWAT

SEAWAT version 4 is a well-known and widely tested 

three-dimensional variable-density groundwater-flow computer program developed by the USGS that couples the 

groundwater-flow code MODFLOW-2000 (Harbaugh and 

others, 2000) with the subsurface contaminant transport code 

MT3DMS (Zheng and Wang, 1999). SEAWAT version 4 can 

explicitly simulate variable-density and variable-viscosity 

effects on groundwater flow based on changes in both 

dissolved-solids concentration and temperature; however, variable-viscosity effects and heat transport were not considered 

in this assessment. The effect of changing density on groundwater flow is simulated in SEAWAT by coupling the groundwater-flow equation and the solute-transport advective-dispersive equation through a linear equation of state that relates 

water density in a cell to the dissolved-solids concentration 

of the water in that cell. Although SEAWAT was originally 

applied on a regional scale to simulate saltwater intrusion 

and submarine groundwater discharge in a coastal aquifer in 

Florida (Langevin, 2000), it has been tested and used in other 

settings from laboratory to regional scales and is applicable to 

noncoastal aquifer systems with areas of increased salinities, 

such as the Edwards aquifer. More information on SEAWAT 

can be found at the U.S. Geological Survey’s SEAWAT Web 

page (U.S. Geological Survey, 2014) and in the user manuals 

(Guo and Langevin, 2002; Langevin and others, 2003; Langevin and others, 2008).

SEAWAT version 4 requires additional information not 

required by MODFLOW-2000 in order to simulate solute 

transport of dissolved solids and density-dependent groundwater flow: advection and dispersion information, the linear 

equation of state that relates water density to dissolved-solids 

concentration, dissolved-solids concentrations at boundary conditions, effective porosity, and initial concentrations. 

Advection was solved using the standard finite-difference 

method with upstream-weighting. Numerical dispersion occurs 

when using the finite-difference upstream-weighting method 

to solve for advection; therefore, longitudinal, and transverse 

dispersivity parameters, as well as molecular diffusion, were 

set to 0.0. The reference freshwater density for SEAWAT version 4 was set to 62.4 pounds per cubic foot (lbs/ft3

) and the 

slope of the linear equation of state was set to 0.714. 

The Source/Sink Mixing (SSM) Package of SEAWAT 

version 4 contains information related to concentrations 

and locations of model sources and sinks (boundary conditions). At a model sink (where the water leaves a cell through 

a spring, simulated by using the DRN or production well, 

simulated by using the WEL), the dissolved-solids concentration of the sink within a cell will equal the concentration of 

water in the cell. At model sources (where water enters a cell 

through net recharge, simulated by using the RCH or from 

the underlying Trinity aquifer, simulated by using the WEL), 

the dissolved-solids concentration is specified in the SSM. 

The constant dissolved-solids concentration boundary condition for the southern model boundary is also specified in the 

SSM and was adjusted during history matching (discussed in 

the “Calibration Though History Matching” section) and the 

concentration of each cell was based on the history-matched 

initial dissolved-solids concentration at that cell. The Basic 

Transport Package (BTN) (Zheng and Wang, 1999) contains 

information about the effective porosity and initial dissolvedsolids concentration arrays in each layer of the model. These 

arrays were created as external arrays in SEAWAT so that 

parameter changes could be made to individual arrays during 

history matching. 

Calibration through History Matching

Models are calibrated through the process of history 

matching. History matching involves adjusting model input 

parameters to improve the agreement between historical 

observations (heads, spring discharges, and dissolved-solids 

concentrations representative of a defined time) and modelsimulated equivalents with the presumption that a model 
Updated Conceptualization and Numerical Model of the Edwards Aquifer 23

that can reproduce past conditions will provide more reliable 

predictions. 

Unfortunately, obtaining a sufficient match does not 

ensure that the history-matched parameter values (1) are 

representative of real-world properties, (2) are unique, because 

the data used for history matching is rarely sufficient to fully 

inform all aspects of the adjustable parameters, or (3) will 

provide the most conservative estimate of the prediction of 

interest (White and others, 2014). If the predictions of interest 

do not resemble the observation dataset, then history matching 

may do little to improve the reliability of predictions (Moore 

and Doherty, 2005). Therefore, calibration through history 

matching is generally a necessary step in the development of 

groundwater models, but it should not be considered the final 

product of the development of the model. Uncertainty quantification of the input parameters and model predictions is needed 

to supplement the history-matched best-fit parameter values. 

The best-fit parameter values are herein defined as those 

values that minimize the sum of squared differences (residuals) between the observation dataset and the model-simulated 

equivalents.

History matching, parameter uncertainty analysis, and 

predictive uncertainty analysis were completed with PEST 

(Doherty, 2005) and PEST++ (Welter and others, 2012), which 

are suites of model-independent, open-source, public-domain 

software that implement several parameter estimation and 

predictive uncertainty analysis algorithms. Use of this software, coupled with explicit mathematical regularization—the 

mathematical processes that help achieve numerical stability in the inverse problem and govern parameter plausibility and acceptability (Doherty and Hunt, 2010)—facilitates 

flexibility of the inverse problem to match the observation 

dataset while accommodating a large number of parameters. 

More importantly, a highly parameterized (Doherty and Hunt, 

2010) approach to history matching facilitates robust uncertainty quantification because more uncertain model inputs 

are explicitly recognized and included in the analyses than 

would have been with few adjustable parameters (Voss, 2011a, 

2011b). Before history matching, preprocessing for PEST was 

needed: (1) to determine parameterization in accordance with 

the model objectives and predictions of interest, (2) to process 

the observation dataset, and (3) to determine realistic expected 

values and ranges of parameters. The initial value and specified acceptable ranges of the adjustable parameters are considered part of the “prior,” which is specification of expert knowledge in the parameter estimation process. Parameter ranges 

represent the 95-percent credible interval for each parameter.

Parameterization 

Parameterization describes how uncertain model inputs 

are discretized in space and time in the history-matching process. The parameterization selected for this modeling analysis 

used a combination of pilot points (Doherty, 2003; Doherty 

and Hunt, 2010; Doherty and others, 2010a) and piece-wise 

zones of uniform property values (Hill and Tiedeman, 2007). 

Pilot points are a method of spatial parameter definition where 

hydraulic properties are assigned to a set of points distributed 

throughout the model domain rather than directly to the grid 

or mesh elements of a numerical model; property values are 

then assigned to model elements through spatial interpolation 

(Doherty, 2003). The selected parameterization recognized 

many more model inputs as uncertain when compared to 

previous models of the Edwards aquifer (Lindgren and others, 

2004; Lindgren, 2006; Lindgren and others, 2011b). Pilot 

points were evenly distributed throughout the model domain, 

with a spacing of approximately 3.5 miles (fig. 11). The quantity of pilot points increased with depth in the model as the 

active extent of layers increase from layer 1 through 

layer 8 (fig. 8). The upper part of the Edwards aquifer 

(model layers 1–3) had 404 pilot points, the middle part of 

the Edwards aquifer (model layer 4) had 433 pilot points, and 

the lower part of the Edwards aquifer (model layers 5–8) had 

454 pilot points. An anisotropic variogram (Isaaks and Srivastava 1989; Oliver and others, 2008) was used to distribute 

hydraulic conductivity from pilot points to the model grid and 

to incorporate apparent anisotropy in the calibrated spatial 

distribution of hydraulic conductivity obtained from Lindgren 

(2006) and the geostatistical analysis of Painter and others 

(2002). Use of an anisotropic variogram with pilot points is 

a more flexible application of expert knowledge compared 

to previous modeling efforts (Lindgren and others, 2004; 

Lindgren, 2006), and also respects the highly stochastic nature 

of the spatial distribution of hydraulic conductivity within the 

Edwards aquifer. The primary axis of anisotropy was specified as 50 degrees east of north, with an anisotropy ratio of 5.0 

(Isaaks and Srivastava, 1989).

Hydraulic and Transport Parameters

Horizontal hydraulic conductivity for the upper and lower 

parts of the Edwards aquifer, vertical hydraulic conductivity for the middle part of the Edwards aquifer, and specific 

storage for the upper and lower parts of the Edwards aquifer (fig. 3) were adjusted during history matching by using 

multiplier pilot points (Dausman and others, 2010a). Instead 

of multipliers, the hydraulic characteristic (representative of 

fault hydraulic conductivity divided by fault width) of each 

horizontal-flow barrier within the Horizontal-Flow Barrier 

Package (HFB) (Harbaugh and others, 2000) was adjusted. 

The only solute transport-related parameters adjusted in the 

history-matching effort by using the BTN were multiplier 

pilot points for effective porosity for the upper and lower parts 

of the Edwards aquifer. During each forward model run, the 

multiplier pilot points are spatially interpolated to create arrays 

with the same dimension as the model input arrays. These 

interpolated arrays of multiplier values are then multiplied by 

the appropriate initial uniform property array, which resulted 

in updated arrays that were then either saved to file for direct 

use in the model or, for select properties, were further adjusted 

through the use of scaling parameters.

The selected parameterization was based on many trial 

iterations of parameterization and history matching. This 

iterative process indicated the need for dissolved-solids 
24 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement KIMBLE COUNTY SUTTON COUNTY

TRAVIS COUNTY GILLESPIE COUNTY

BLANCO COUNTY

HAYS COUNTY

KERR COUNTY

KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL COUNTY

BANDERA COUNTY

GUADALUPE COUNTY

GONZALES COUNTY

BEXAR COUNTY

MEDINA COUNTY

UVALDE COUNTY

KINNEY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA COUNTY

KARNES COUNTY FRIO COUNTY ZAVALA COUNTY MAVERICK COUNTY

EDWARDS COUNTY

30°29°

98° 99° 100°

Base from U.S. Geological Survey digital data Texas Centric Mapping SystemAlbers equal area projectionNorth American Datum of 1983

0 10 20 KILOMETERS

0 10 20 MILES

Recharge zone boundary (modified from Puente, 1978)Model boundaryPilot point

Active model areaInactive model area

EXPLANATION

Figure 11. Locations of pilot points where hydraulic properties are adjusted during history matching, San Antonio region, Texas.
Updated Conceptualization and Numerical Model of the Edwards Aquifer 25

concentration threshold parameters and scaling parameters 

to best represent the hydraulic property variation associated 

with the different water-quality zones (fig. 4) of the Edwards 

aquifer. The updip extent of the transition zone is where 

mixing between fresh and saline water has caused formation 

of high-permeability, high-porosity preferential pathways of 

groundwater flow (Hovorka and others, 1998). The saline zone 

of the Edwards aquifer (fig. 4) is relatively stagnant compared 

to the freshwater zone, and thought to have relatively low 

permeability and porosity associated with largely unaltered 

carbonates (Deike, 1990; Groschen, 1994; Hovorka and others, 1998). To represent these relations in the history-matching 

effort, two threshold parameters were estimated, with each 

threshold parameter having three scaling parameters: hydraulic conductivity, specific storage, and effective porosity. The 

cumulative effect of this parameterization process is a flexible 

and automated property adjustment procedure that approximates the location of preferential flow pathways within the 

Edwards aquifer.

The first parameter (high_thresh) represents the dissolved-solids threshold concentration above which aquifer 

properties need to be scaled to be more representative of the 

area near the transition-zone interface where karstic alteration 

has developed high-permeability, high-porosity preferential 

flow areas. For dissolved-solids concentrations greater than

high_thresh, the hydraulic conductivity, specific storage, and 

effective porosity need to be scaled up, that is, increased by 

multiplying by a scaling factor greater than 1.0. Because little 

knowledge exists about the exact nature of this complex process, the concentration where this threshold occurs was treated 

as an adjustable parameter in history matching so that the bestfit value could be informed by the observation dataset. The 

best-fit value of high_thresh was 329 mg/L and the associated 

scaling parameters for hydraulic conductivity, specific storage, 

and effective porosity were 1.5, 1.7, and 2.0, respectively. 

These best-fit scaling parameters indicate that more groundwater needs to be transmitted through the transition zone 

(higher hydraulic conductivity) at a lower advective velocity 

(higher effective porosity). The second threshold parameter 

(low_thresh) represents the dissolved-solids threshold concentration above which hydraulic conductivity, specific storage, 

and effective porosity need to be scaled down to be more representative of the unaltered, lower permeability, lower porosity 

saline zone. The best-fit value of low_thresh was 2,074 mg/L 

and the associated scaling parameters for hydraulic conductivity, specific storage, and effective porosity were 0.0011, 0.001, 

and 0.75. The saline-zone scaling parameters indicate that less 

groundwater needs to be transmitted in the unaltered saline 

part of the Edwards aquifer. Implementation of this concentration-based parameterization allows for the flexibility of inclusion of additional expert knowledge of the aquifer in history 

matching, which effectively supplements the sparse transition 

and saline-zone observation dataset. Furthermore, the best-fit 

estimates of these threshold parameters and scaling parameters 

are consistent with the hydrogeologic and geochemical understanding of the existence of preferential-flow areas near the 

“bad-water line” of 1,000 mg/L dissolved-solids (which falls 

within the range of high_thresh to low_thresh estimates of 

329 to 2,074 mg/L). 

The initial spatial distribution of hydraulic (horizontal 

hydraulic conductivity, vertical hydraulic conductivity, and 

specific storage) and transport (effective porosity) properties 

in the model domain were specified as uniform and, where 

applicable, were assigned an initial value that is the mean of 

the associated calibrated spatial property distribution from 

previous models (Lindgren and others, 2004; Lindgren, 2006). 

Lower and upper bounds of the 95-percent credible interval of 

each parameter were used to set the range of potential values, 

which allows the parameter values maximum flexibility within 

a range of values based on expert knowledge.

The horizontal hydraulic conductivity pilot points 

(Doherty, 2003; Doherty and Hunt, 2010; Doherty and others, 

2010a) for the upper and lower parts of the Edwards aquifer 

and the vertical hydraulic conductivity pilot points for the 

middle part of the Edwards aquifer were adjusted during the 

history-matching effort to better reproduce the observation 

dataset. The initial values of these multiplier pilot points were 

1, and the prior 95-percent credible interval of these multiplier 

pilot points ranged from 1.0 × 10-5 to 1.0 × 105. The initial 

spatial distribution of horizontal hydraulic conductivity for the 

upper and lower parts of the Edwards aquifer was a uniform 

value of 1,500 ft/day based on the average hydraulic conductivity from the calibrated model of Lindgren (2006). The 

initial spatial distribution of vertical hydraulic conductivity for 

the middle part of the Edwards aquifer was a uniform value of 

3.0 ft/day based on the lower bound of the range of calibrated 

hydraulic conductivity values in Lindgren (2006). The best-fit 

horizontal hydraulic conductivity (fig. 12A and 12 C) distributions for the upper and lower parts of the Edwards aquifer 

(including the concentration threshold scaling) and the best-fit 

vertical hydraulic conductivity (fig. 12B) distribution of the 

middle part of the Edwards aquifer indicate considerable 

spatial variability in best-fit values of conductivities needed 

for history matching. In general, hydraulic conductivities 

are smaller farther downdip in the model as expected. Also, 

although vertical hydraulic conductivities of the middle part 

of the Edwards aquifer are smaller than horizontal hydraulic 

conductivities of the upper and lower parts of the Edwards 

aquifer, spatial variability still exists.

As the model is a fully confined approximation of the 

Edwards aquifer, the specific storage parameters control storage changes in the aquifer. Specific storage multiplier pilot 

points were adjusted for the upper and lower parts of the 

Edwards aquifer, and the initial value of the specific storage 

distribution was 1.1 × 10-6 1/ft and was uniform. The initial 

value of all specific storage multiplier pilot points was 1, and 

the prior 95-percent credible interval of these multiplier pilot 

points allowed specific storage to range from 1.0 × 10-9 1/ft to 

1.0 × 10-4 1/ft. The best-fit specific storage distributions for the 

upper and lower parts of the Edwards aquifer (fig. 13), which 

include the concentration threshold scaling, indicate considerable spatial variability in values of specific storage, with 
26 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement Figure 12. The history-matched best-fit distributions of hydraulic conductivities of the Edwards aquifer, San Antonio region, Texas. 

A, Horizontal hydraulic conductivity of 

the upper part of the Edwards aquifer. 

B, Vertical hydraulic conductivity of the middle part of the Edwards aquifer. 

C, Horizontal hydraulic conductivity of the lower part of the 

Edwards aquifer.

Base from U.S. Geological Survey digital data Texas Centric Mapping SystemAlbers equal area projectionNorth American Datum of 1983

Base from U.S. Geological Survey digital data Texas Centric Mapping SystemAlbers equal area projectionNorth American Datum of 1983

Base from U.S. Geological Survey digital data Texas Centric Mapping SystemAlbers equal area projectionNorth American Datum of 1983

KIMBLE COUNTY SUTTON COUNTY

TRAVIS

COUNTY GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY KERR COUNTY KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY

GUADALUPE

COUNTY

GONZALES

COUNTY

MEDINA

COUNTY

UVALDE COUNTY KINNEY

COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

FRIO COUNTY ZAVALA COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

KIMBLE COUNTY SUTTON COUNTY

TRAVIS

COUNTY GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY KERR COUNTY KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY

GUADALUPE

COUNTY

GONZALES

COUNTY

MEDINA

COUNTY

UVALDE COUNTY KINNEY

COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

FRIO COUNTY ZAVALA COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

KIMBLE COUNTY SUTTON COUNTY

TRAVIS

COUNTY GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY KERR COUNTY KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY

GUADALUPE

COUNTY

GONZALES

COUNTY

MEDINA

COUNTY

UVALDE COUNTY KINNEY

COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

FRIO COUNTY ZAVALA COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

30°29°

98° 99° 100°

30°29°

98° 99° 100°

30°29°

98° 99° 100°

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

A. Horizontal hydraulic conductivity of upper part of Edwards aquifer

B. Vertical hydraulic conductivity of middle part of Edwards aquifer

C. Horizontal hydraulic conductivity of lower part of Edwards aquifer

Hydraulic conductivity, in feet per day 0.001 1,000 5,000 10,000 50,000 100,000 or largerModel extentActive model-area boundary

EXPLAN

ATION

Edwards aquifer Upper part (model layers 1–3) Middle part (model layer 4) Lower part (model layers 5–8)
Updated Conceptualization and Numerical Model of the Edwards Aquifer 27

the highest values occurring in the recharge zone (fig. 13B) 

in Kinney and Uvalde counties. The specific storage for the 

middle part of the Edwards aquifer was not used as an adjustable parameter and was specified as a uniform value of 1.1 × 

10-6 1/ft.

The hydraulic effect of faults within the Edwards aquifer 

is pronounced but uncertain; some faults are believed to be 

barriers to flow, whereas others are believed to concentrate 

flow preferentially (Groschen, 1994; Hovorka and others, 

1998). In recognition of this uncertainty, the hydraulic characteristic of each of the HFB-simulated faults (same as those 

included in the diffuse-flow model [Lindgren, 2006]) was 

adjusted during the history-matching process. The historymatched HFB hydraulic characteristic value of Lindgren 

(2006) for each fault was used as the initial value in this analysis. The prior 95-percent credible interval of possible hydraulic 

characteristic values was specified to reflect unknown specifics 

of any fault; therefore, values were allowed to range from 0 to 

1.0 × 1010 per day. The best-fit HFB hydraulic characteristic 

values ranged from 0 to 10,430 per day.

Effective porosity multiplier pilot points were adjusted 

for the upper and lower parts of the Edwards aquifer. The 

initial value of the effective porosity distribution was 0.25 and 

was uniform. The initial value of all effective porosity multiplier pilot points was 1, and the prior 95-percent credible interval of these multiplier pilot points allowed effective porosity 

to range from 0.025 to 0.5 (2.5 percent to 50 percent). The 

best-fit effective porosity (fig. 14) distributions for the upper 

and lower parts of the Edwards aquifer, which include the 

effects of the dissolved-solids concentration threshold scaling 

parameters, indicate considerable spatial variability within the 

aquifer. The highest values of effective porosity are within the 

updip area of the transition zone (fig. 4) where the application 

of the high_thresh scaling parameter has increased the 

values (fig. 14). The effective porosity for the middle part of 

the Edwards aquifer was specified a uniform value of 0.25 as 

the effective porosity of this layer was assumed to minimally 

affect solute transport (0.25 is representative of the mean of 

the prior).

Model Sources and Sinks

The terms “model sources” and “model sinks” respectively refer to the model boundary conditions that control 

water entering and leaving the aquifer. Including these 

boundary condition elements in the history matching process is important because these elements strongly affect the 

model simulated results and are largely uncertain. Groundwater recharge is a primary model source. The selected 

recharge parameterization attempts to characterize several 

sources of recharge uncertainty. The USGS has published 

monthly recharge estimates for the Edwards aquifer that are 

distributed across surface-water basins (Lindgren and others, 

2004; Lindgren, 2006; Slattery, 2004). These estimates were 

prepared using the best available information; nevertheless, 

the estimates are uncertain, largely as a result of the complex 

and transient nature of recharge within the highly karstified 

Edwards aquifer. To represent the uncertainty of recharge 

within the history-matching effort, the recharge parameterization included two distinct sets of parameters that were used in 

sequence. The result of this parameterization was used to build 

a RCH for input into the model.

 The first set of recharge parameters were scaling 

parameters that were applied to each surface-water basin (fig. 

15) that overlies the outcrop of the Edwards aquifer. These 

parameters corresponded to individual basins within the 

recharge zone from Lindgren (2006). These scaling parameters 

were assigned to either a wetter month when rainfall is higher 

than average (April–June, September–November) or drier 

month when rainfall is less than average (January–March, 

July, August, December) (Larkin and Bomar, 1983). The initial 

value of all month-dependent recharge multiplier parameters was 1, and the prior 95-percent credible bounds of the 

multipliers ranged from 0.25 to 1.75. The best-fit wetter and 

drier month surface-water basin recharge multipliers (fig. 15) 

indicate that multipliers for individual surface-water basins are 

quite different between the wetter and drier months. Decreases 

in recharge are seen in the red color range and increases in 

recharge are seen in the blue color range. 

The second set of recharge parameters attempted to 

capture the potential for annual variability in recharge estimates by using a global scaling parameter for each calendar 

year from 1999 through 2009. That is, for each year in the 

model simulation period, a single multiplier parameter was 

applied to the input recharge arrays for each simulated month 

in that year. The initial value of all annual recharge multiplier 

parameters was 1, and the prior 95-percent credible bounds of 

the multipliers ranged from 0.25 to 1.75. The best-fit annual 

recharge scaling parameters indicate considerable variability 

from year to year with no discernable pattern (fig. 16).

In addition to recharge originating from precipitation, 

groundwater flow (flux) into the Edwards aquifer from the 

underlying Trinity aquifer (figs. 2–3) may also be a substantial 

model source, especially during drought conditions; however, 

the magnitude and spatial distribution of this flux is highly 

uncertain (Wong and others, 2013). The calibrated spatial 

distribution of the Trinity aquifer flux from Lindgren (2006) 

was used as the initial flux for this modeling effort. To capture 

this source of uncertainty in the history-matching process, 

the Trinity aquifer flux was scaled by a group of multiplier 

parameters. These multiplier parameters consisted of cells of 

the Trinity aquifer flux grouped into 10 zones. This parameterization provided the ability to spatially adjust the Trinity 

aquifer flux during history matching. The Trinity aquifer flux 

from Lindgren (2006) was scaled by the multiplier parameters 

before being used by the model in each forward run. The 

initial value of all Trinity aquifer flux multiplier parameters 

was 1, and the prior 95-percent credible interval bounds of 

the multipliers ranged from 0.01 to 100. Flux from the Trinity 

aquifer into the Edwards aquifer is represented in model 

layer 8 (fig. 5) at the northern boundary of the active model 

area. The best-fit spatial distribution of Trinity aquifer flux 
28 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

30°

29°

100° 99° 98°

30°

29°

100° 99° 98°

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

A. Specific storage of upper part of Edwards aquifer

B. Specific storage of lower part of Edwards aquifer

Specific storage of upper

 part (model layers 1–3)

 of Edwards aquifer, in

 per foot

 Less than 0.0000007

 0.000005

 0.000001

 0.00005

 0.0001

Model extent

Active model-area

 boundary

EXPLANATION

Specific storage of lower

 part (model layers 5–8)

 of Edwards aquifer, in

 per foot

 Less than 0.0000007

 0.000005

 0.000001

 0.00005

 0.0001

Model extent

Active model-area

 boundary

EXPLANATION

Figure 13. The history-matched best-fit distributions of specific storage of the Edwards aquifer, San Antonio region, Texas for 

A, the upper part of the Edwards aquifer, and B, the lower part of the Edwards aquifer.
Updated Conceptualization and Numerical Model of the Edwards Aquifer 29

Figure 14. The history-matched best-fit distributions of effective porosity of the Edwards aquifer, San Antonio region, Texas for 

A, the upper part of the Edwards aquifer, and B, the lower part of the Edwards aquifer.

30°

29°

100° 99° 98°

30°

29°

100° 99° 98°

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

A. Effective porosity of upper part of Edwards aquifer

B. Effective porosity of lower part of Edwards aquifer

Effective porosity of 

 lower part (model 

 layers 5–8) of

 Edwards aquifer,

 unitless

 0.16

 0.21

 0.25

 0.31

 0.50

Model extent

Active model-area

 boundary

Spring—Tail points in

 direction of flow

EXPLANATION

q

Effective porosity of 

 upper part (model 

 layers 1–3) of

 Edwards aquifer,

 unitless

 0.13

 0.20

 0.25

 0.42

 0.50

Model extent

Active model-area

 boundary

Spring—Tail points in

 direction of flow

EXPLANATION

q

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras

Springs

Leona Springs

Barton Springs 

San Marcos

Springs

q

q

q

q

q

q

q

q

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras

Springs

Leona Springs

Barton Springs 

San Marcos

Springs

q

q

q

q

q

q

q

q
30 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement Figure 15. The history-matched best-fit distributions of recharge multipliers for the surface-water basins overlying the Edwards aquifer, San Antonio region, Texas. 

A, Recharge multipliers for the wetter months of April–June, September–November. 

B, Recharge multipliers for the drier months of January–March, July, August, and 

December.

30°29°

98° 99° 100°

Base from U.S. Geological Survey digital data, 2012 Texas Centric Mapping SystemAlbers equal area projectionNorth American Datum of 1983

Recharge zone (outcrop) of the Edwards

 aquifer modified from Puente, 1978

KIMBLE COUNTY SUTTON COUNTY

TRAVIS COUNTY GILLESPIE COUNTY

BLANCO COUNTY

HAYS COUNTY

KERR COUNTY

KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL COUNTY

BANDERA COUNTY

GUADALUPE COUNTY GONZALES

COUNTY

BEXAR COUNTY

MEDINA COUNTY

UVALDE COUNTY

KINNEY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA COUNTY

KARNES COUNTY FRIO COUNTY ZAVALA COUNTY MAVERICK COUNTY

EDWARDS COUNTY

Medina

Lake

Canyon Lake

South Llano River

Dry Draw

P ni ot Creek

Dry Frio River

Long Kent Creek

Johnson Creek

Guadalupe River

Pedernales River

Miller Creek Blanco River

S al ughter Creek

Colorado Rive r

Onion Creek

Cedar Creek

York Creek

San Marcos Rive r

Plum Creek

Picosa Creek

l E C m er ke

Chaparrosa Creek

Turkey Creek

Palo Blanco Creek

Elm Creek

Nueces River

Blanco Creek

Sabinal River

Seco Creek

Leona River

Frio River

dno H Co keer

San Miguel Cree k Black Cre ek

Medina Riv er

Atascosa River

Lucas Creek

Hondo Creek

San Antonio Creek

Cibolo Creek

Ecleto Creek

Sandies Cree k

E ml Creek

Elm Creek

1 2

3 4 5 6

7

8

9

10

12

13

11

16

14

15

0 10 20 KILOMETERS

0 10 20 MILES

A. Recharge multipliers for wetter months (April–June, September–November)

Recharge multiplier, unitless 0.25

1.00

 1.75Active model areaModel extentContributing area—Basin extent shown  terminates at downdip extent of  Edwards aquifer recharge zone Major river basin boundary Stream basin or contributing area boundaryRecharge zone boundary

EXPLANATION

Edwards aquifer recharge subzones San Antonio segment Nueces-West Nueces River Frio-Dry Frio River Sabinal River Area between Sabinal and Medina Rivers Medina River

 Area between Medina River and

 Cibolo-Dry Comal Creek

 Cibolo Creek and Dry Comal Creek 

 Guadalupe River Blanco River Barton Springs segment Onion Creek Little Bear Creek Bear Creek Slaughter Creek Williamson Creek Barton Creek Colorado River

1

23

4

567

89

10

11

12131415

16
Updated Conceptualization and Numerical Model of the Edwards Aquifer 31 Figure 15. The history-matched best-fit distributions of recharge multipliers for the surface-water basins overlying the Edwards aquifer, San Antonio region, Texas. 

A, Recharge multipliers for the wetter months of April–June, September–November. 

B, Recharge multipliers for the drier months of January–March, July, August, and 

December.—Continued

30°29°

98° 99° 100°

Base from U.S. Geological Survey digital data, 2012 Texas Centric Mapping SystemAlbers equal area projectionNorth American Datum of 1983

Recharge zone (outcrop) of the Edwards

 aquifer modified from Puente, 1978

KIMBLE COUNTY SUTTON COUNTY

TRAVIS COUNTY GILLESPIE COUNTY

BLANCO COUNTY

HAYS COUNTY

KERR COUNTY

KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL COUNTY

BANDERA COUNTY

GUADALUPE COUNTY GONZALES COUNTY

BEXAR COUNTY

MEDINA COUNTY

UVALDE COUNTY

KINNEY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA COUNTY

KARNES COUNTY FRIO COUNTY ZAVALA COUNTY MAVERICK

EDWARDS COUNTY

COUNTY

Medina

Lake

Canyon Lake

South Llano River

Dry Draw

P ni ot Creek

Dry Frio River

Long Kent Creek

Johnson Creek

Guadalupe River

Pedernales River

Miller Creek Blanco River

S al ughter Creek

Colorado Rive r

Onion Creek

Cedar Creek

York Creek

San Marcos Rive r

Plum Creek

Picosa Creek

l E C m er ke

Chaparrosa Creek

Turkey Creek

Palo Blanco Creek

Elm Creek

Nueces River

Blanco Creek

Sabinal River

Seco Creek

Leona River

Frio River

dno H Co keer

San Miguel Cree k Black Cre ek

Medina Riv er

Atascosa River

Lucas Creek

Hondo Creek

San Antonio Creek

Cibolo Creek

Ecleto Creek

Sandies Cree k

E ml Creek

Elm Creek

1 2

3 4 5 6

7

8

9

10

12

13

11

16

14

15

0 10 20 KILOMETERS

0 10 20 MILES

B. Recharge multipliers for drier months (January–March, July, August, December)

Recharge multiplier, unitless 0.25

1.00

 1.75Active model areaModel extentContributing area—Basin extent shown  terminates at downdip extent of  Edwards aquifer recharge zone Major river basin boundary Stream basin or contributing area boundaryRecharge zone boundary

EXPLANATION

Edwards aquifer recharge subzones San Antonio segment Nueces-West Nueces River Frio-Dry Frio River Sabinal River Area between Sabinal and Medina Rivers Medina River

 Area between Medina River and

 Cibolo-Dry Comal Creek

 Cibolo Creek and Dry Comal Creek 

 Guadalupe River Blanco River Barton Springs segment Onion Creek Little Bear Creek Bear Creek Slaughter Creek Williamson Creek Barton Creek Colorado River

1

23

4

567

89

10

11

12131415

16
32 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

(fig. 17B) displays more variability when compared to the 

initial spatial distribution of Trinity aquifer flux (which is the 

same as the flux from Lindgren, 2006) (fig. 17A). 

Springs are major model sinks where water flows out 

of the Edwards aquifer (table 1). For modeling purposes, 

springs were represented as outflow (drain) features, which is 

a gross simplification of these complex hydrologic features. 

As a result of this simplification, the model inputs associated with these boundary conditions are uncertain. Thus, the 

drain conductance and altitude for each spring in DRN of the 

model were adjusted during history matching. The calibrated 

values for drain conductance and altitude (table 11 of Lindgren 

and others, 2004, p. 131) were used as initial values in this 

modeling analysis. The prior 95-percent credible interval for 

the drain conductance parameters ranged from 0 to 1.0 × 1010

feet squared per day (ft2/day). The lower bound of the prior

95-percent credible interval for the drain altitude parameters 

was set at the altitude of the bottom of the model cell that 

contained the drain; the upper bound was 1.0 × 1010 ft above 

National Geodetic Vertical Datum of 1929 (NGVD 29). The 

best-fit drain conductance and drain altitude values for Leona, 

San Pedro, San Antonio, Comal, San Marcos, and Barton 

Springs (table 1) are similar to calibrated values for these 

springs from previous models (Lindgren and others, 2004; 

Lindgren, 2006). The best-fit drain conductance for this model 

for Las Moras Springs was two orders of magnitude higher 

than for the model of Lindgren and others (2004), but this 

spring is near the active model boundary in the western part 

of the model area and data from this spring were not used for 

history matching.

Initial Dissolved-Solids Concentrations

As explained in the “History-Matching Phase–Initial 

Conditions–Initial Dissolved-Solids Concentrations” section of this report, initial dissolved-solids concentrations 

were based on water-quality and geophysical log data, and a 

relatively large amount of dissolved-solids concentration data 

exist for the freshwater zone when compared to the sparse 

amount of dissolved-solids concentration data for the transition and saline zones (fig. 9). Because of this disparity and the 

associated uncertainty, the initial dissolved-solids concentration distribution was parameterized and adjusted through use 

of pilot points of initial concentration multipliers for the upper, 

middle, and lower parts of the Edwards aquifer. As the initial 

concentrations were adjusted, the concentration assigned to 

the southern specified-concentration boundary condition was 

also automatically adjusted to account for the uncertainty in 

the specification of this boundary condition (fig. 18). The bestfit initial dissolved-solids concentration distributions for the 

upper, middle, and lower parts of the Edwards aquifer indicate 

some changes with depth (such as in Frio County in the lower 

part of the Edwards aquifer) where the observation dataset has 

informed the initial concentrations. However, in general, the 

initial concentrations were similar between the upper, middle, 

and lower parts of the Edwards aquifer because there was not 

enough information in the observation dataset to uniquely 

estimate the initial concentrations in the different parts of the 

aquifer (fig. 18).

Historical Observation Dataset

More than 8,000 observations were used for history 

matching; data consisted of groundwater heads, spring discharge, dissolved-solids concentrations, and dissolved-solids 

Table 1. The history-matched best-fit Drain (DRN) package 

altitude and conductance parameters that control spring 

discharge from the model, San Antonio region, Texas.

[DRN, MODFLOW Drain Package; ft, feet; NGVD, National Geodetic 

Vertical Datum of 1929; ft2/day, feet squared per day]

Springs in the San 

Antonio region (fig. 1)

DRN Altitude, ft above 

NGVD 29

DRN Conductance, ft2

/day

Las Moras Springs 1,009 104,626,550 

Leona Springs 799 142,869 

San Pedro Springs 684 216,845 

San Antonio Springs 664 726,731 

Comal Springs 558 5,139,129 

San Marcos Springs 551 4,701,207 

Barton Springs 403 2,160,535 

Figure 16. The history-matched best-fit annual recharge 

multipliers, 1999–2009.

1999 2000 2001 2002 2003 2004

Year

2005 2006 2007 2008 2009

Annual recharge multiplier

1.8

1.6

1.4

1.2

1.0

0.8

0.6

0.4

0.2

0
Updated Conceptualization and Numerical Model of the Edwards Aquifer 33

Figure 17. Inflow of groundwater from the Trinity aquifer (flux) into the overlying Edwards aquifer, San Antonio region, Texas. A, Initial 

flux. B, History-matched best-fit flux.

30°

29°

100° 99° 98°

30°

29°

100° 99° 98°

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

Base from U.S. Geological Survey digital data 

Texas Centric Mapping System

Albers equal area projection

North American Datum of 1983

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

SUTTON COUNTY KIMBLE COUNTY

TRAVIS

COUNTY

GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY

KERR COUNTY

KENDALL

COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY GUADALUPE

COUNTY

GONZALES

COUNTY

KINNEY UVALDE COUNTY MEDINA COUNTY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

ZAVALA COUNTY FRIO COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

A. Initial Trinity aquifer flux into the overlying Edwards aquifer

B. History-matched best-fit Trinity aquifer flux into the overlying Edwards aquifer

Inflow of groundwater

 from the Trinity aquifer

 (flux) into the overlying 

 Edwards aquifer, in cubic

 feet per day

 Less than 2,000

 6,000

 12,000

Active model area

Model extent

Recharge zone boundary 

 of the Edwards aquifer

 (modified from Puente,

 1978) 

EXPLANATION
34 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement Base from U.S. Geological Survey digital data Texas Centric Mapping System Albers equal area projection North American Datum of 1983 Base from U.S. Geological Survey digital data Texas Centric Mapping System Albers equal area projection North American Datum of 1983 Base from U.S. Geological Survey digital data Texas Centric Mapping System Albers equal area projection North American Datum of 1983 KIMBLE COUNTY SUTTON COUNTY

TRAVIS

COUNTY GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY KERR COUNTY KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY

GUADALUPE

COUNTY

GONZALES

COUNTY

MEDINA

COUNTY

UVALDE COUNTY KINNEY

COUNTY

WILSON COUNTY

DEWITT

COUNTY

DEWITT

COUNTY

DEWITT

COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

FRIO COUNTY ZAVALA COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

KIMBLE COUNTY SUTTON COUNTY

TRAVIS

COUNTY GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY KERR COUNTY KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY

GUADALUPE

COUNTY

GONZALES

COUNTY

MEDINA

COUNTY

UVALDE COUNTY KINNEY

COUNTY

WILSON COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

FRIO COUNTY ZAVALA COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

KIMBLE COUNTY SUTTON COUNTY

TRAVIS

COUNTY GILLESPIE COUNTY

BLANCO

COUNTY

HAYS

COUNTY KERR COUNTY KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL

COUNTY

BANDERA COUNTY

GUADALUPE

COUNTY

GONZALES

COUNTY

MEDINA

COUNTY

UVALDE COUNTY KINNEY

COUNTY

WILSON COUNTY

ATASCOSA

COUNTY

KARNES COUNTY

FRIO COUNTY ZAVALA COUNTY

MAVERICK

COUNTY

EDWARDS

COUNTY

BEXAR

COUNTY

30°29°

98° 99° 100°

30°29°

98° 99° 100°

30°29°

98° 99° 100°

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

0 10 20 KILOMETERS

0 10 20 MILES

A. Initial dissolved-solids concentrations in upper part of Edwards aquifer

B. Initial dissolved-solids concentrations in middle part of Edwards aquifer

C. Initial dissolved-solids concentrations in lower part of Edwards aquifer

Dissolved-solids concentrations, in milligrams per liter

67

 15,000 180,941Model extentActive model-area boundarySpringTransect well

EXPLAN

ATION

q

Edwards aquifer Upper part (model layers 1–3) Middle part (model layer 4) Lower part (model layers 5–8)

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras Springs

Leona Springs

San Marcos Springs

Barton Springs 

q

q

qq

q

q

q

q

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras Springs

Leona Springs

San Marcos Springs

Barton Springs 

q

q

qq

q

q

q

q

San Pedro Springs

San Antonio Springs

Hueco Springs

Comal Springs

Las Moras Springs

Leona Springs

San Marcos Springs

Barton Springs 

q

q

qq

q

q

q

q

Figure 18. The history-matched best-fit distributions of initial dissolved-solids concentrations of the Edwards aquifer, San Antonio region, Texas for 

A, the upper part of the 

Edwards aquifer, 

B, the middle part of the Edwards aquifer, and 

C, the lower part of the Edwards aquifer.
Updated Conceptualization and Numerical Model of the Edwards Aquifer 35

mass constraints at withdrawal (production) wells (where mass 

refers to the solute mass) (fig. 19). Annual total discharge volumes of the five history-matched springs in the model (Leona, 

San Pedro, San Antonio, Comal, and San Marcos) were also 

used for history matching. 

The dissolved-solids mass constraints represent the 

observation that production wells, or withdrawal wells, in the 

Edwards aquifer do not produce poor-quality (high dissolvedsolids concentration) water. These constraints were formulated 

so that the total mass of dissolved solids for a given well 

was equal to the specified freshwater dissolved-solids concentration of the Edwards aquifer times the volume of water 

extracted from the well during the model simulation period. 

Observation weights were subjectively specified to guide 

the history-matching process to reproduce the observations 

that most resemble the model predictions. As explained by 

Doherty and Welter (2010, p. 2) “subjective decisions must 

then be made on weighting schemes to apply to the calibration 

data set when fitting model outputs to that data set in order to 

estimate parameters. As estimated parameter values depend on 

these decisions, these estimates, as well as assessments made 

of the uncertainty associated with these estimates, inherit a 

subjective component.” The observation weights were specified so that the composite, weighted objective function—a 

function of residuals (differences between observations and 

model-simulated equivalents)—at the start of the historymatching process was composed of an equal contribution from 

each of the following types of observations: 

1. Groundwater-head time series for well J-17; 

2. Remaining groundwater-head time series at another 

192 wells; 

3. Spring discharge time series for Comal and San Marcos 

Springs; 

4. Remaining spring discharge time series for Leona, San 

Pedro, and San Antonio Springs; 

5. Total annual discharge volumes for Comal and San 

Marcos Springs; 

6. Remaining total annual discharge volumes for Leona, 

San Pedro, and San Antonio Springs; 

7. Total mass of dissolved-solids (mass constraints) for 

production wells, and 

8. Dissolved-solids concentration data.

Groundwater Head Observations

Groundwater heads for wells within the active model 

domain (4,232 head observations, fig. 19, 20A) for the period 

of 1999–2009 were used for history matching. Well construction information was used, when available, to assign heads to 

distinct layers instead of using the well depth. Well J-17 was 

more heavily weighted because it is the focus of regulatory 

scrutiny and heads measured at this well are known to be 

statistically correlated with discharge at Comal Springs (Miller 

and Long, 2006). 

Spring Discharge Observations

Five major springs were used for history matching: 

Leona, San Pedro, San Antonio, Comal, and San Marcos 

Springs (fig. 19, 20B). Comal and San Marcos Springs were 

more heavily weighted because predictions will be made at 

these springs and because the discharge at these two springs 

is important for ecological and recreational purposes. Total 

annual discharge volumes for each of the five springs were 

calculated from the spring discharge series and were used as 

an additional data type in the observation dataset during history matching. A total of 636 spring discharge observations 

and 55 total annual volume observations were included in the 

observation dataset.

Dissolved-Solids Concentration

Each dissolved-solids concentration observation (693 

observations) (fig. 19, 20C) was included in history matching; each observation was used twice, with different weights 

applied to each observation as part of a standard and inverseweighting process. One standard weight was assigned proportional to the magnitude of the concentration, whereas the other 

was assigned inversely proportional to the magnitude of the 

concentration. Including the inverse-weighted concentration 

observations ensures that the history-matching effort appropriately takes into account low dissolved-solids concentrations 

that would otherwise be masked by large dissolved-solids 

concentrations. A total of 693 standard-weighted and 693 

inverse-weighted dissolved-solids concentration observations 

were included in the observation dataset.

Constraints on Total Mass of Dissolved Solids at Production 

Wells

Constraints on the total mass of dissolved solids removed 

from water from withdrawal (production) wells (fig. 19) were 

specified to minimize the extraction of water with concentrations above 325 mg/L (the background freshwater concentration), because simulated wells should not be withdrawing 

brackish or saline water. For each simulated stress period, the 

volume of water withdrawn at a production well was multiplied by the simulated dissolved-solids concentration in the 

same model cell. These values are then summed for all stress 

periods. This calculation yields a total mass of dissolved solids 

removed by each production well. The “observed” value 

for each production well is the product of the background 

concentration times the total volume extracted from each well, 

effectively enforcing a condition that only freshwater exists at 

production wells. Formulation of these types of constraints is 

valid because the only source of dissolved solids in the model 

originates from the saline zone. A total of 1,973 mass constraints were included in the observation dataset. 
36 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement 30° 29° 98° 99° 100° Base from U.S. Geological Survey digital data Texas Centric Mapping System Albers equal area projection North American Datum of 1983

0 10 20 KILOMETERS

0 10 20 MILES

Groundwater withdrawal well locationActive model areaModel extent

EXPLANATION

q

Recharge zone boundary (modified from Puente, 1978)

Approximate line of 1,000 milligrams per liter

dissolved-solids concentration (A.L. Schultz, 

 written commun., 2000)

Spring—Tail points in direction of flow

Observation well, by type Transect Dissolved-solids concentration Groundwater head

KIMBLE COUNTY SUTTON COUNTY

TRAVIS COUNTY GILLESPIE COUNTY

BLANCO COUNTY

HAYS COUNTY

KERR COUNTY

KENDALL COUNTY

REAL COUNTY

CALDWELL

COUNTY

COMAL COUNTY

BANDERA COUNTY

GUADALUPE COUNTY

GONZALES COUNTY

BEXAR COUNTY

MEDINA COUNTY

UVALDE COUNTY

KINNEY COUNTY

WILSON COUNTY

DEWITT

COUNTY

ATASCOSA COUNTY

KARNES COUNTY FRIO COUNTY ZAVALA COUNTY MAVERICK COUNTY

EDWARDS COUNTY

San 

Pedro Springs

San Antonio Springs

Comal Springs

Leona Springs

San Marcos Springs

q

q

q

q

q

Bexar County index well (J–17)

Figure 19. Locations of observation data used for history matching (groundwater withdrawals, springs, dissolved-solids concentrations, and groundwater heads), San Antonio region, Texas, 1999–2009. 
Updated Conceptualization and Numerical Model of the Edwards Aquifer 37

Simulated spring discharge,

in cubic feet per second

–100 0 100 200 300 400 500 600

Simulated log10 dissolved-solids concentration,

in milligrams per liter

1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5

Simulated groundwater head, in feet

200 400 600 800 1,000 1,200 1,400

–100

0

100

200

300

400

500

600

–150

–100

–50

0

50

100

150

Observed spring discharge,

in cubic feet per second

Residual values (observed spring discharge

minus simulated spring discharge), 

in cubic feet per second

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

5.5

–1.0

–0.5

0

0.5

1.0

Observed log10 dissolved-solids concentration,

in milligrams per liter

Residual log10 values (observed dissolved-solids

concentration minus simulated dissolved-solids

concentration), in milligrams per liter

200

400

600

800

1,000

1,200

1,400

–200

–100

0

100

200

Observed groundwater head, in feet

Residual values (observed groundwater head

minus simulated groundwater head), in feet

B. Spring discharge (636 observations)

A. Groundwater head (4,232 observations)

C. Dissolved-solids concentration (693 observations)

EXPLANATION

Zero-residual line

Observed to simulated data point

Residual data point

Correlation coefficient

Logarithm to base 10

R=0.999

log10

1:1 line

1:1 line

1:1 line

R=0.988

R=0.999

R=0.832

Figure 20. Observed values compared to simulated values and residuals for A, groundwater heads, B, spring discharges, and 

C, dissolved-solids concentrations, San Antonio region, Texas.
38 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

Goodness of Fit

In general, the best-fit parameters yield an acceptable fit 

with the observation dataset for the history-matched model 

(1999–2009). As an indication of goodness of fit, scatter plots 

of observed values compared to model-simulated equivalents 

were prepared with equal-value (1:1) lines, also known as 

line-of-equality plots. These scatterplots (fig. 20) indicate a 

good agreement between observations and model-simulated 

equivalents, and the associated residuals indicate little or no 

bias over the range of model-simulated values. The correlation 

coefficient [equation 5.3.2 in Doherty (2015)] was calculated 

for each scatterplot and is shown on each plot (fig. 20). The 

calculated correlation coefficients were 0.999, 0.988, and 

0.832 for groundwater heads, spring discharges, and dissolved-solids concentrations, respectively. The fit is generally 

accepted to be better the closer the correlation coefficient is to 

1.0. Additionally, comparison of observed and model-simulated time series for key prediction locations (Comal Springs, 

San Marcos Springs, and well J-17) indicates that the model 

is generally capable of reproducing important model outputs 

over the range of hydrologic conditions observed during the 

model simulation period (fig. 21). 

The final, best-fit parameter values are within the acceptable (prior) 95-percent credible interval and, where applicable, 

are in general agreement with aquifer-specific expert knowledge of the Edwards aquifer system and results of previous 

history-matched models (Lindgren and others, 2004; Lindgren, 

2006). The best-fit spatial distributions of horizontal hydraulic conductivity (fig. 12), in general, agree with the geostatistically derived horizontal hydraulic conductivity spatial 

distribution of Painter and others (2002), which was used in 

the model of Lindgren and others (2004, fig. 8) and with the 

spatial distribution of history-matched hydraulic conductivity of Lindgren (2006, fig. 4). In general, regions of relatively 

high and low hydraulic conductivities in this model (fig. 12) 

correspond to relatively high and low hydraulic-conductivity 

regions in the previous models (except for the high hydraulic 

conductivity area in the recharge zone at the border between 

Bexar and Comal Counties). This was also the case for the 

best-fit specific storage results (fig. 13) when compared to 

the history-matched storage zone results of previous models 

(fig. 24 of Lindgren and others, 2004). Previous models 

focused on the best-fit spatial distributions of these hydraulic 

parameters because these were the only model inputs used 

to understand the water budget or regional flow; whereas 

in this model, the best-fit spatial property distributions of 

these parameters were supplemented with linear uncertainty 

analysis, which quantified both parameter and predictive 

uncertainty. 

Whereas traditional measures used to evaluate calibration results indicate the model is acceptably history matched, 

many of the uncertain parameters that were explicitly recognized during the history-matching process are likely to 

remain uncertain even after history matching. Unfortunately, 

if the predictions of interest are sensitive to these uncertain 

parameters, then considerable prediction uncertainty may also 

remain following history matching. 

Parameter Uncertainty—Schur’s Complement 

for Linear-Based Conditional Uncertainty 

Propagation 

The prior—the known or expected stochastic character 

of adjustable parameters determined from expert knowledge, literature values, field tests, previous models, and other 

sources of information—is depicted in figure 22 as the dashed 

(log-transformed) Gaussian distributions. In addition to the 

best-fit parameter values derived from history matching (figs. 

12–18, table 1), the posterior (after history-matching) uncertainty of model parameters was also estimated using linear 

uncertainty analysis. Linear uncertainty analysis is predicated 

on the linearity assumption between model parameters and 

observations, which also implies an assumption of Gaussian 

distributions of parameters, which means the primary quantitative metric for uncertainty in a linear framework is variance 

(Doherty and others 2010b). Specifically, Schur’s complement 

(Meyer, 2000) for conditional uncertainty propagation (White 

and others, 2015) is:

∑θ = ∑θ – ∑θ Ј T(Ј∑θ Ј T + ∑εε)-1Ј∑θ

(4)

where 

∑θ is the prior parameter covariance matrix, 

Ј is the Jacobian matrix, and 

 ∑εε is the covariance matrix of measurement 

noise. 

The matrix ∑θ of equation 4 is the posterior parameter 

covariance matrix, which can be seen as the prior parameter 

covariance matrix minus the conditioning provided by the 

observations; conditioning is embodied in the linear mapping 

provided by the Jacobian matrix. As such, use of equation 4 

assumes a linear relation between adjustable parameters and 

model-simulated observation equivalents and that parameter 

and measurement noise uncertainty can be described by a multivariate Gaussian (or log-Gaussian) distribution. Fienen and 

others (2010) discusses a complete derivation of equation 4 

from Bayes’ equation. The diagonal elements of of equation 4 

contain the post-history matching (posterior) parameter uncertainty, expressed as variances (squares of standard deviations) 

(fig. 22, shaded distributions). 

Comparison of the prior (before history matching) and 

posterior (after history matching) parameter uncertainty plots 

(fig. 22) (that is, the variances associated with different model 

parameters) depicts how the information within the observation dataset used for history matching informs many parameters. The information transfer from the observation dataset to 

the parameters is seen in two ways: (1) through a reduction in 

posterior variances (reduction of width and increase in height 

of distribution), and (2) through shifts of the posterior distribution along the x-axis when compared to the prior distribution. 
Updated Conceptualization and Numerical Model of the Edwards Aquifer 39

550

500

450

400

350

300

250

200

150

710

700

690

680

670

660

650

640

630

450

400

350

300

250

200

150

100

50

0

Groundwater head, in feet Spring discharge, in cubic feet per second Spring discharge, in cubic feet per second

2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010

Year

A. Comal Springs

B. San Marcos Springs

C. Bexar County index well (J–17)

EXPLANATION

Simulated

Observed

Figure 21. Simulated and observed time-series of A, Comal Springs discharge, B, San Marcos Springs discharge, and C, Bexar 

County index well J-17 groundwater head, San Antonio region, Texas.
40 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

Figure 22. Distributions of prior (before history matching) and posterior (after history matching) parameter uncertainty.

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing density Increasing density Increasing density

Increasing density

FLUX MULTIPLIER PARAMETERS

PILOT POINT MULTIPLIER PARAMETERS

Posterior parameter uncertainty distribution for each 

 parameter within the group—The intensity of the 

 shaded area reflects the overlapping of distributions

 Aquifer

 Upper part of 

 the Edwards

 Middle part of 

 the Edwards

EXPLANATION

Prior parameter uncertainty 

 distribution—Uniform for 

 all parameters in this 

 group

 Lower part of the Edwards

 Trinity—Trinity aquifer

 flux into Edwards 

 aquifer

Trinity aquifer flux into

Edwards aquifer

Specific storage

Specific storage

Horizontal hydraulic 

conductivity

Vertical hydraulic 

conductivity

Horizontal hydraulic 

conductivity

Effective porosity

Effective porosity

Initial dissolved-solids

concentration

Initial dissolved-solids

concentration

Initial dissolved-solids

concentration
Updated Conceptualization and Numerical Model of the Edwards Aquifer 41

Figure 22. Distributions of prior (before history matching) and posterior (after history matching) parameter uncertainty.—Continued

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing logarithm of 

parameter value 

Increasing density Increasing densityIncreasing density Increasing density

Increasing density

FLUX MULTIPLIER PARAMETERS

THRESHOLD PARAMETERS SCALING PARAMETERS

Posterior parameter uncertainty

 distribution for each parameter 

 within the group—The intensity

 of the shaded area reflects the

 overlapping of distributions

 Zone

 Recharge

 Brackish-water transition

 Saline

Prior parameter uncertainty 

 distribution—Uniform for all 

 parameters in this group

EXPLANATION

Dissolved-solids

concentration 

Dissolved-solids

concentration 

Annual recharge

Specific storage

Specific storage

Horizontal hydraulic 

conductivity

Horizontal hydraulic 

conductivity

Effective porosity

Effective porosity

Recharge (drier months) Recharge (wetter months)

Drier months (January–March, July, August, December)

Wetter months (April–June, September–November) 

High_thresh

Low_thresh
42 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

The initial (prior) distributions are the same for all individual 

parameters within a parameter type (for example, the distribution of hydraulic conductivity multiplier pilot point values was 

the same before history matching). The shaded distributions 

represent individual posterior (after history matching) distributions of all parameters within each parameter type (for example, the individual drier-month recharge multipliers for each 

surface-water basin). Therefore, the effects of history matching on individual parameters are depicted in figure 22 through 

comparison of the prior distribution and the posterior distributions for each parameter type. Where the prior and posterior 

variances of a given parameter are similar, it can be inferred 

that the observation dataset contributed little to inform that 

particular parameter. For example in figure 22, the effective 

porosity multiplier parameters for the upper and lower parts of 

the Edwards aquifer have posterior distributions similar to the 

prior, which indicates that the observation dataset contributed 

little to inform these parameters and did not decrease their 

uncertainties through history matching. The dissolved-solids 

concentration threshold parameters were well-informed by 

the observation dataset as their posterior distributions were 

much narrower (lower variance) than their prior distributions. 

Scaling parameters of hydraulic conductivity, effective porosity, and specific storage for the transition zone (fig. 22) were 

all informed by the observation dataset, as evidenced by the 

difference between the prior and posterior variances. Salinezone scaling parameters, on the other hand, were not informed 

by the observation dataset for effective porosity and specific 

storage. This is not surprising considering that few observation data have been collected that might inform the porosity 

or storage of the saline zone. However, uncertainty in the 

hydraulic conductivity scaling factor for the saline zone was 

reduced. Note that all of the parameter values shown in these 

distributions are log-transformed and that most parameters are 

multiplier values (or scaling parameters) and not necessarily 

the final values used as inputs to the model. For example, data 

values for hydraulic conductivity of the distribution for the 

upper part of the Edwards aquifer are multiplier pilot point 

values before scaling and kriging to the grid. The high_thresh 

and low_thresh concentration scaling threshold parameters 

are not multiplier values but are the actual history-matched 

threshold values (log-transformed) (in pounds per cubic feet 

instead of mg/L). 

Resulting posterior drier-month, wetter–month, and 

annual recharge multiplier parameter variances are important 

to understanding how reliably recharge is estimated and implemented in the model. The shifts of the posterior distributions 

to the left and right of the prior distribution indicate that some 

zones in the model needed less or more simulated recharge as 

a model input in order to achieve an acceptable history-matching result (fig. 22). Many of the best-fit recharge parameters 

(fig. 22) are not statistically different from the initial values 

specified in the history-matching effort; the widths of the distributions of the different recharge parameters did not decrease 

substantially as a result of history matching. Recharge from 

rainfall is the driving force behind groundwater flow and 

heads in the aquifer; therefore, an increased understanding of 

the complex Edwards aquifer recharge process would benefit 

model development by potentially decreasing the uncertainty 

of this parameter.

The history-matching effort was helpful in informing the 

parameters that control spring discharge, namely, the DRN 

altitude and conductance parameters for each spring. These 

distributions were not shown in figure 22 because the distributions were narrow, indicating these parameters were strongly 

informed by the observation dataset. 

Uncertainty Assessment of 1950–56 

Drought Conditions on Brackish-Water 

Movement within the Edwards Aquifer 

The best-fit parameters inferred through the historymatching effort are not likely to be “correct” parameters, but 

instead represent minimum error variance parameters. Therefore, in the predictive phase, results of hypothetical scenarios 

from the predictive model based on the best-fit parameters 

should be evaluated with explicit consideration of parameter uncertainty. Traditional groundwater modeling analyses 

would use the best-fit parameters to yield one estimate of the 

model predictions through a single forward run of the predictive model. However, any remaining posterior uncertainty in 

prediction-sensitive parameters means that many parameter 

sets can reproduce the observation dataset as well as the 

best-fit parameter set, but each of these parameter sets may 

yield different values of the model predictions. Therefore, it 

is imperative that posterior parameter uncertainty be propagated to the model predictions to evaluate uncertainty in the 

model predictions because rigorous assessment of prediction 

uncertainty provides important insights into the reliability of 

model predictions, compared to using only the single, best-fit 

parameters.

Uncertainty analysis was completed using equation 4 

combined with predictive sensitivity vectors (vectors that 

contain the sensitivity of the prediction with respect to all 

adjustable parameters) to quantify predictive uncertainty. To 

propagate parameter uncertainty to predictions of interest, the 

following relations (White and others, 2015; Doherty and others, 2010b) were used:

σs = y T∑ θ y ; σs = y T∑ θ y 2 2 (5)

where

y is the sensitivity vector of prediction s with 

respect to the adjustable parameters,

 and σs

2 and σs

2 are the prior and posterior variances of 

prediction s, respectively and all other 

terms are previously defined in equation 4.
Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement within the Edwards Aquifer 43

Although regional changes in the position of the transition zone during the predictive period are of ultimate interest, 

the three key quantitative targets for this predictive model 

and uncertainty analysis were (1) changes in dissolved-solids 

concentrations at 25 production wells near the transition-zone 

interface, (2) discharge at Comal and San Marcos Springs, and 

(3) the head at well J-17. These three quantitative targets are 

helpful for evaluating the effect of sustained drought on the 

San Antonio segment of the Edwards aquifer. 

Predictive Phase 

In the predictive phase, the model extent, boundary locations and types, aquifer geometry, grid layering, and discretization of the predictive model are the same as in the historymatched model (1999–2009). The only two aquifer stresses 

that were changed for the predictive model were recharge and 

groundwater withdrawals at production wells. The Trinity 

aquifer flux through the northern boundary of model layer 8 

(fig. 5) was the same as that of the history-matched model. 

Drains used to simulate major springs with the DRN were 

used with altitude and conductance values from the historymatched model as well. Spatial property distributions of 

hydraulic conductivity, effective porosity, and specific storage 

were also from the history-matched model. Monthly stress 

periods were used for the 7 years of drought (1950–56) and 

initial conditions for heads and dissolved-solids concentrations 

were obtained from the initial heads and dissolved-solids concentrations of the history-matched model. These initial conditions were used because this predictive model is not a hindcast 

of 1950s groundwater conditions in the Edwards aquifer, but 

is a hypothetical forecast prediction of effects of the combination of modern withdrawal rates and drought-of-record rainfall 

conditions on the three predictive targets.

Drought Recharge Conditions 1950–56

The history-matched model period of 1999–2009 does 

not best represent drought-of-record conditions in southcentral Texas. The current (2015) driest year on record was 

2011 (Winters, 2013). The drought of the 1950s has been used 

as an important simulation event or target with other models of the Edwards aquifer and other groundwater models in 

Texas (Lindgren and others, 2004; Winters, 2013). However, 

this period is especially important for the Edwards aquifer 

because the only period of zero flow from Comal Springs was 

recorded from June 13, 1956, to November 4, 1956 (Maclay, 

1995). Monthly recharge estimates of the 1950s have previously been used for the Edwards aquifer models (Lindgren and 

others 2004; Lindgren 2006; Slattery, 2004). These recharge 

estimates were used in the RCH of this predictive model 

using the same zonation scheme as the 1999–2009 historymatched model. The drier-month, wetter-month, and annual 

recharge multiplier parameters were also implemented in this 

predictive model to determine effects of uncertainty in these 

multipliers from the history-matched model on uncertainty in 

the predictions.

Withdrawals by Wells

Groundwater withdrawals for the predictive model were 

implemented using the WEL using the same well locations as 

those implemented in the history-matched model. However, 

the individual monthly withdrawal rate for each well was 

specified to be the 85th percentile of all withdrawal rates for 

a given month, which produced 12 withdrawal rates for each 

well. This same one-year long, monthly withdrawal time series 

for each well was used for each year of the 1950–56 period. 

Therefore, the predictive model represents a hypothetical scenario of the drought of the 1950s combined with hypothetical 

increased groundwater withdrawal rates. 

The focus of the predictive model was to evaluate 

changes in hydrologic conditions in a regional perspective 

across the San Antonio segment of the Edwards aquifer, but 

production wells near the transition-zone interface (fig. 10) 

were specifically targeted to estimate the change in dissolvedsolids concentration in water withdrawn from the wells from 

the beginning to the end of the model simulation. Concentration changes at these wells were used as model predictions 

because regional maps could not show the movement of the 

transition zone at a meaningful scale. These 25 wells were 

selected because they are near the transition-zone interface 

(fig. 10) and nonpotable water could be withdrawn during 

drought conditions.

Groundwater fluxes and cumulative volumes simulated 

by the predictive model (fig. 23), based on best-fit historymatched parameter distributions, indicate that drought 

conditions (1950–56) eventually lead to a depletion of spring 

discharge as groundwater heads decrease until the withdrawals 

by wells are eventually balanced by depletion in storage (fig. 

23). Compared to the history-matched water budget, predictive 

cumulative volumes of simulated recharge and spring discharge are much smaller, depletion of storage is much larger, 

and withdrawals by wells are similar. 

Prediction Uncertainty—Schur’s Complement 

for Linear-Based Conditional Uncertainty 

Propagation

Linear uncertainty analysis was completed in this assessment and is discussed in great detail in Doherty and others 

(2010b, appendix 4). Linear uncertainty analysis is based on 

the assumption that model-simulated equivalents (model outputs) are linearly related to parameters (model inputs) through 

the single Jacobian sensitivity matrix, which contains the 

sensitivities of each model-simulated equivalent to each model 

parameter. A second, implied assumption is that the parameter 

variability and interaction is characterized by a multivariate 

Gaussian distribution (or the log of the parameter variability 
44 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

Recharge

Storage

Well 

 withdrawals

Spring 

 discharge

EXPLANATION

History-matched groundwater flow (flux) History-matched cumulative volume

Predictive groundwater flow (flux) Predictive cumulative volume

Year

1950 1951 1952 1953 1954 1955 1956 1957

Year

1950 1951 1952 1953 1954 1955 1956 1957

1999

2000

2001

2002

2003

2004

2005

2006

2007

2008

2009

2010

1999

2000

2001

2002

2003

2004

2005

2006

2007

2008

2009

2010

–10,000

–5,000

5,000

0

10,000

Volume, in acre-feet times 1,000

–20,000

–10,000

10,000

0

20,000

–15,000

–5,000

15,000

5,000

Groundwater flow, in acre-feet per day

–10,000

–5,000

5,000

0

10,000

Volume, in acre-feet times 1,000

–20,000

–10,000

10,000

0

20,000

–15,000

–5,000

15,000

5,000

Groundwater flow, in acre-feet per day

Figure 23. Groundwater flows (fluxes) and cumulative volumes for the history-matched model (1999–2009) and the predictive model 

(1950–56), San Antonio region, Texas. Positive values represent groundwater gains, and negative values represent groundwater losses.
Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement within the Edwards Aquifer 45

is Gaussian). The goal of predictive linear uncertainty analysis is to quantify predictive uncertainty of the model based 

on parameter uncertainty using information gained from the 

observation dataset through the history-matching effort. This 

uncertainty analysis is done under the hypothesis that the 

history-matching effort reduces the uncertainty of the model 

predictions from the prior uncertainty by transferring information from the observation dataset to prediction-sensitive 

parameters.

The prior parameter uncertainty for the predictive 

analysis was inherited from the parameter uncertainty analysis. These prior uncertainty estimates were based on expert 

knowledge of a range of the potential values for the parameters before the history-matching effort (refer to the “Parameterization” section). Using equation 5, the prior uncertainty 

estimates were propagated to several key model predictions, 

including dissolved-solids concentration changes in water at 

wells near the transition-zone interface, spring discharges at 

Comal and San Marcos Springs, and head at well J-17. 

To investigate the effects of history-matching information 

on the prediction, equation 4 was combined with equation 5 to 

yield posterior prediction uncertainty estimates. The resulting posterior predictive uncertainty was used to estimate the 

upper and lower bound of the 95-percent credible interval for 

each model prediction. For the production wells of predictive 

interest (fig. 10), the upper bound represents a conservative 

estimate of the expected concentration change under the hypothetical recurrence of drought conditions. 

Predictive Uncertainty Results

Predictive uncertainty analysis results at the 25 production wells near the transition-zone interface indicate that the 

prior uncertainty of model input parameters (based on expert 

knowledge) yielded an upper bound of the 95-percent credible interval of dissolved-solids concentrations that exceeds 

the secondary drinking water standards of 1,000 mg/L (Texas 

Commission on Environmental Quality, 2013) for many of 

the wells (fig. 24). The initial concentration of dissolved 

solids at these wells was 325 mg/L; therefore, an increase 

of greater than 675 mg/L (the y-axis of fig. 24) represents a 

final concentration greater than 1,000 mg/L. However, the 

history-matching effort provided key information to inform 

prediction-sensitive model parameters, and, therefore, contributed to a substantial decrease of the upper bound of the 

posterior 95-percent credible interval (when compared to the 

prior) to below the secondary drinking water dissolved-solids 

concentration standards (fig. 24). Reductions in dissolvedsolids concentration change were on the order of 400 mg/L (at 

wells 5, 12, 16, and 25) to 1,300 mg/L (at wells 2, 8, 9, 11, and 

13). The reduction in uncertainty in regards to this prediction 

Figure 24. The upper bound of the 95-percent credible interval for the prior and posterior predictive uncertainty of dissolved-solids 

concentration change during the 1950–56 predictive model period at production wells of interest.

Production well of predictive interest* (see fig. 10 for location)

* Prediction of dissolved-solids concentration change during simulated 

 7-year (1950–56) drought of record conditions made at these production wells

Upper bound of 95-percent credible interval— Change during

 the 1950–56 predictive model period

 Prior (before history-matching) uncertainty of prediction

 Posterior (after history-matching) uncertainty of prediction

EXPLANATION

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

Change in dissolved-solids concentration,

in milligrams per liter

0

200

400

600

800

1,000

1,200

1,400

1,600

1,800
46 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

implies that this prediction can be made with some certainty 

and that those parameters that control this prediction are 

informed by the observation dataset. Even though predictive 

uncertainty was reduced for this prediction, dissolved-solids 

concentration changes were still greater than zero, indicating 

a minimal increase in concentration is likely at these 

25 production wells during the 7-year simulation period. 

This minimal concentration increase indicates a small potential for movement of the transition zone near these wells during the hypothetical recurrence of 1950–56 drought-of-record 

rainfall conditions with higher-than-average groundwater 

withdrawals by wells. 

Predictive results of total spring discharge during the 

seven–year period, as well as the prediction of head at well 

J-17, were much less reliable than the prediction of changes in 

dissolved-solids concentrations at the production wells. The 

predicted total spring discharges for Comal and San Marcos 

Springs were about 401,000 and 165,930 acre-feet (acre-ft), 

respectively. The predicted head at the end of the predictive model period for well J-17 was 460 ft above NGVD 29. 

More important than the predicted values are the widths of the 

95-percent credible intervals of these predictions. For example, the upper bound of the 95-percent credible interval of total 

discharge at Comal and San Marcos Springs at the end of the 

model simulation are about 6,080,600 and 1,244,750 acreft, respectively. The upper bound of the 95-percent credible 

interval of head at well J-17 at the end of the model simulation is 3,280 ft. These upper bounds are an order of magnitude 

larger than the actual predictions, which implies that (1) these 

predictions of spring discharge at Comal and San Marcos 

Springs and head at well J-17 simulated with this model are 

not reliable, and (2) parameters that control these predictions are not informed well by the observation dataset during 

history matching, even though the history-matching effort 

yielded parameters to reproduce spring discharge and heads at 

these locations during the history-matching period. This large 

posterior uncertainty is an example of how history matching is 

successful in yielding best-fit parameters, which allows 

for sufficient simulation of heads and spring discharge 

(see fig. 21), but is unsuccessful at reducing the uncertainty of 

parameters that control predictions made at these same locations. In other words, history matching is sufficient to yield 

best-fit parameter estimates but not necessarily sufficient to 

produce reliable predictions. Furthermore, spring discharges 

at these two springs and head at well J-17 represent a cumulative effect of upstream conditions over a larger extent within 

the Edwards aquifer (and a longer time). Therefore, many 

more parameters (with their own uncertainties) are potentially 

controlling the spring discharge and head predictions than the 

prediction of dissolved-solids concentration change at the prediction wells, and, therefore contributing to a large posterior 

uncertainty. Spring discharges and heads in the Edwards aquifer respond much more quickly and fluctuate with a greater 

magnitude than the dissolved-solids concentrations. The 

apparent lack of transition zone movement is the result of the 

stable density configuration of the Edwards aquifer transition 

zone and saline zone. That is, the brackish-water transition and 

saline zones inherently resist moving, with most groundwater 

flowing parallel to the transition-zone interface rather than 

normal to the interface. It should be noted that the predictions 

of spring discharge and head were absolute predictions, which 

are more difficult to accurately make than the difference prediction of dissolved-solids concentration change (difference 

from the beginning of the model simulation to the end of the 

model simulation) at the production wells. 

Model Limitations 

Models are inherently flawed and therefore, have limitations as they are simplified simulators of a much more complex natural system (Box and Draper, 1987). Model limitations 

need to be discussed and model results need to be understood 

in the context within which the model was designed and developed. Furthermore, any interpretations of model results should 

be in accordance with the assumptions made in developing 

that model, which must also coincide with the prediction of 

interest. 

Model error is a term that is used to describe the discrepancy between the model and the more complex, natural 

system it simulates (White and others, 2014). Sources of 

model error might include horizontal and vertical discretization error, parameter lumping and vertical averaging of 

properties, and the assumption of an equivalent porous media. 

Errors in simulated heads were expected in the recharge 

zone of the model because a confined approximation was 

used to simulate all layers of the model. The recharge zone is 

unconfined; using a confined approximation reduced numerical convergence issues in the model, decreased model run 

times, and most importantly, the focus of the modeling effort 

was on transition-zone movement within the confined part of 

the aquifer, not the recharge zone. Recharge from rainfall is 

the driving force behind groundwater flow and heads in the 

aquifer; therefore, an increased understanding of this process 

would benefit model development by potentially decreasing 

the uncertainty of this parameter. Decreasing the uncertainty 

of recharge parameters might be accomplished through a more 

rigorous simulation of the rainfall-to-recharge landscape partitioning process with a hydrologic model. Some information 

on groundwater withdrawals was only available in aggregated 

spatial and temporal scales, which may cause localized spatial 

and temporal error in simulated heads near these wells. Where 

well-construction information was not available, groundwater 

withdrawals were simulated from the entire thickness of the 

Edwards aquifer. In other cases, the screened interval was 

used to determine the layers in which groundwater withdrawal 

was simulated. Only major springs were simulated in the 

model, but many smaller springs are known to exist where 

Edwards aquifer water is potentially discharging. Localscale dissolution features and faulting at springs may not be 

fully represented because of the scale of the model and grid 
Summary 47

discretization, and these features may be better simulated 

in a local-scale, finely discretized model. Given the highly 

karstified character of the Edwards aquifer, the assumption of 

an equivalent porous media is likely to be violated at smaller 

spatial scales. Therefore, using the model to understand localscale groundwater flow and transport processes may not be 

appropriate. 

Few data from the transition zone and saline zone are 

available to characterize geologic structure, faults, aquifer 

properties, and distributions of heads and dissolved-solids 

concentrations. Some of these information deficits could be 

reduced through more stratigraphic determinations of Edwards 

aquifer tops and bottoms from existing geophysical logs and 

conversions of geophysical logs to dissolved-solids estimates. The assumption of the southern boundary as a no-flow, 

constant-concentration boundary condition is sufficient for the 

predictions made herein with this model, but small amounts of 

water might flow across this boundary changing the concentration of dissolved solids at this boundary over time.

This modeling analysis assumes simulated dissolvedsolids concentrations are only related to transition-zone 

and saline-zone processes. Furthermore, more complicated 

geochemical processes and reactions were not simulated 

within this model. In essence, water-density is dependent upon 

dissolved-solids concentration, which is simulated as a conservative constituent in the model.

Water density is also dependent upon temperature and 

regions of elevated temperature are known to occur within the 

subsurface at depths in the Edwards, resulting from geothermal heat flux (Groschen and Buszka, 1997). Furthermore, 

viscosity of groundwater can also be affected by temperature 

changes and dissolved-solids concentration changes (Langevin 

and others, 2008), but heat-transport and variable-viscosity 

effects were not simulated in this analysis. The effects of temperature increases on density and the effects of temperature 

increases and dissolved-solids concentration on viscosity were 

assumed to be negligible in the transition-zone water.

Perhaps the largest limitation of the current model is the 

relative short period simulated for the history-matched model. 

Parameter knowledge gained from the observation dataset 

(as visualized in changes in the posterior parameter uncertainty when compared to the prior) only represents conditions 

from 1999 through 2009 because this is the period for which 

higher-resolution model input datasets were available to support a complex, transient, variable-density modeling analysis. 

Construction and history matching of a full period-of-record 

model would encompass more historical observations (longer 

record for heads, major spring discharges, and dissolvedsolids concentrations) of the Edwards aquifer, but would 

require additional assumptions about the spatial and temporal 

distribution of groundwater withdrawal rates. However, the 

increased size of the observation dataset for a full period-ofrecord model has the potential to further reduce uncertainty of 

parameters—most importantly—potentially those parameters 

that control the predictions of spring discharge and head at 

well J-17. 

Uncertainty analysis presented herein, while comprehensive in the inclusion of uncertain model parameters and 

processes, is predicated on a linearity and normality assumption, which carries with it its own assumptions of multivariate 

Gaussian distributions of parameters and predictions as well as 

an assumed linear relation between observations and parameters which may not be completely accurate (Doherty and 

others, 2010b). Whereas it has been demonstrated that linear 

analysis in the context of variable-density modeling is capable 

of providing robust uncertainty estimates (Dausman and others, 2010b), more advanced uncertainty quantification techniques, such as null-space Monte Carlo (Doherty and others, 

2010b), may provide additional information and understanding 

of parameter and predictive uncertainty.

The predictive model in this assessment was based on 

a hypothetical recurrence of drought-of-record (1950–56) 

rainfall and higher-than-average groundwater withdrawals by 

wells. The predictions of interest are specific to these rainfall 

and withdrawal conditions at the 25 production wells, Comal 

and San Marcos Springs, and well J-17 where predictions 

were made. Results of the predictions of interest should not be 

extrapolated to other locations within the Edwards aquifer or 

to other hydrologic conditions. 

Summary

In 2010, the U.S. Geological Survey, in cooperation 

with the San Antonio Water System, began a study to assess 

the brackish-water movement within the Edwards aquifer 

(more specifically the potential for brackish-water encroachment into wells near the interface between the freshwater and 

brackish-water transition zones, referred to in this report as 

the transition-zone interface) and effects on spring discharge 

at Comal and San Marcos Springs under drought conditions 

using a numerical model. The quantitative targets of this study 

are to predict the effects of higher-than-average groundwater 

withdrawals and drought-of-record (1950–56) rainfall conditions on (1) dissolved-solids concentration changes at production wells near the transition-zone interface, (2) total spring 

discharge at Comal and San Marcos Springs and (3) the head 

at Bexar County index well J-17 (hereinafter referred to as 

“well J-17”). The predictions of interest, and the parameters 

implemented into the model, were evaluated to quantify their 

uncertainty so results of predictions could be presented in 

terms of a 95-percent credible interval.

The model area covers both the San Antonio and Barton 

Springs segments of the Edwards aquifer although the historymatching effort was focused on the San Antonio segment. The 

diffuse-flow model, which forms the basis for the model in 

this assessment, is primarily based on a conceptualization in 

which flow in the aquifer is predominately through a network 

of numerous small fractures and openings. In this assessment, 

the major changes from the diffuse-flow model that were 
48 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

incorporated into the history-matched model (1999–2009) are 

summarized: 

1. More active model cells were used to simulate the 

brackish-water transition zone and saline zone downdip 

of the freshwater zone.

2. The one-layer model was converted to an 8-layer model 

to better represent vertical aquifer heterogeneity and 

anisotropy as well as simulate variable-density flow and 

transport. Discussion of the model is in terms of the 

“upper” (layers 1–3), “middle” (layer 4), and “lower” 

(layers 5–8) parts of the Edwards aquifer based on the 

hydrostratigraphy.

3. SEAWAT version 4 was used to explicitly simulate 

variable-density flow and transport.

4. The history-matching period was updated to 1999–2009.

5. Parameter estimation, through history matching, was 

executed using highly parameterized inversion with 

PEST to calibrate to heads, discharges at major springs, 

and dissolved-solids concentrations at wells.

The parameterization of the numerical model used a 

combination of pilot points and piece-wise zones of uniform 

aquifer property values. The final parameterization was based 

on many iterations of parameterization and history matching, 

which required use of dissolved-solids concentration threshold 

and scaling parameters to best represent the different hydraulic properties associated with different water-quality zones of 

the Edwards aquifer. To represent the relations between water 

chemistry and karst flow-system properties in the historymatching effort, two threshold parameters were estimated, 

with each threshold parameter having three associated scaling 

parameters, one for hydraulic conductivity, one for specific 

storage, and one for effective porosity. The cumulative effect 

of this parameterization process is a flexible and automated 

hydraulic property adjustment procedure that approximates 

karstic alteration that has occurred within the Edwards aquifer 

along the brackish-water transition zone between freshwater 

and saline water. The best-fit estimates of these threshold 

and scaling parameters are consistent with the hydrogeologic 

and geochemical understanding of the existence of highly 

transmissive preferential-flow areas within the vicinity of the 

“bad-water line” of 1,000 mg/L dissolved-solids concentration (which falls within the range of high_thresh to low_thresh 

estimates of 329 to 2,074 mg/L).

In general, the best-fit parameters yield an acceptable 

fit with the observation dataset. Traditional measures of 

goodness-of-fit, such as the line-of-equality plots and correlation coefficients, indicate good agreement between observations and model-simulated equivalents and little or no bias 

over the range of observed values. Additionally, comparison 

of observed and model-simulated time series for several key 

observation series indicate the model is generally capable of 

reproducing important model outputs over a range of hydrologic conditions. 

In addition to best-fit parameter values derived from 

history matching, the uncertainty of model parameters was 

also estimated using linear uncertainty analysis. Comparison 

of prior (before history matching) and posterior (after history 

matching) variances of parameters indicates that the information within the observation dataset used for history matching 

informs many parameters. This is seen in two ways: 

(1) through a reduction in posterior variances (reduction of 

width and increase in height of distribution) or (2) shifts of the 

posterior distribution along the x-axis when compared to the 

prior distribution. The effective porosity multiplier parameters 

have posterior distributions similar to the prior, which indicates that the observation dataset did not inform this parameter well, and therefore did not decrease its uncertainty. The 

concentration threshold parameters were well-informed by the 

observation dataset as their posterior distributions were much 

narrower than their prior distributions. The transition-zone 

scaling parameters of hydraulic conductivity, effective porosity, and specific storage were all informed by the observation 

dataset, as evidenced by the difference between the prior and 

posterior variances. Saline-zone scaling parameters, on the 

other hand, were not informed by the observation dataset for 

effective porosity and specific storage. This is not surprising considering the lack of data within the saline zone. The 

history-matching effort was most helpful in informing the 

parameters that control spring discharge, namely, the altitude 

and conductance parameters for each spring. These posterior 

distributions were not shown because they were so narrow. 

The reduction in width implies that these parameters were 

strongly informed by the observation dataset. 

Resulting posterior drier-month, wetter-month, and 

annual recharge multiplier parameter variances are important 

to understanding how well recharge is estimated and implemented within the model. The shifts of the posterior distributions left and right indicate that there were zones where less or 

more water was needed in the model. The widths of the distributions were not decreased substantially, indicating that many 

of the best-fit recharge parameters are not statistically different 

from the initial values specified in the history-matching effort. 

Recharge from rainfall is the driving force behind groundwater 

flow and heads in the aquifer; therefore, an increase in understanding of this process would benefit model development by 

potentially decreasing the uncertainty of this parameter.

The predictive model in this assessment simulated a 

hypothetical recurrence of 1950–56 drought conditions and 

higher-than-average groundwater withdrawals by wells to 

provide insight into the potential effects of these conditions 

on dissolved-solids concentration changes at production wells 

near the transition-zone interface, discharges at Comal and 

San Marcos Springs, and groundwater head at well J-17. The 

only two aquifer stresses changed for the predictive model 

were that of aquifer recharge and groundwater withdrawals by 

wells. The 85th percentile of all withdrawal rates for a given 

month for each well was calculated, and a monthly time series 

was created and implemented into the model well package. 

This same one-year monthly withdrawal time series for each 
References Cited 49

well was used for each year of the 1950–56 predictive model 

period.

The uncertainty of the predictions made with the predictive model was quantified using linear uncertainty analysis. 

The goal of predictive linear uncertainty analysis was to 

quantify predictive uncertainty of the model based on parameter uncertainty using information gained from the observation 

dataset through the history-matching effort. The expectation 

is that the history-matching effort reduces the uncertainty of 

the model predictions from the prior uncertainty by transferring information from the observation dataset to predictionsensitive parameters. The resulting posterior predictive 

uncertainty was used to estimate the upper and lower bound of 

the 95-percent credible interval for each model prediction. For 

the production wells of predictive interest, the upper bound 

represents a conservative estimate of the expected concentration change under the prediction conditions. 

The uncertainty assessment of the predictive model 

yielded 95-percent credible intervals of dissolved-solids concentration changes at 25 production wells near the transitionzone interface, discharges at Comal and San Marcos Springs, 

and head at well J-17. Results at the 25 production wells near 

the transition-zone interface indicate that the uncertainty of 

model input parameters based on expert knowledge (prior 

uncertainty) yielded an upper bound of the 95-percent credible 

interval of dissolved-solids concentrations that exceeds the 

secondary drinking water standards of 1,000 mg/L for many of 

the wells. However, the history-matching effort provided key 

information to inform prediction-sensitive model parameters, 

and, therefore, contributed to a substantial decrease of the 

upper bound of the 95-percent credible interval to below the 

secondary drinking water standards. Reductions in dissolvedsolids concentration change were on the order of 400 mg/L 

to 1,300 mg/L. The reduction in uncertainty in regards to 

this prediction implies that this prediction can be made with 

some certainty using this model and that those parameters 

which control this prediction are informed by the observation dataset. Even though predictive uncertainty was reduced 

for this prediction, dissolved-solids concentration changes 

were still greater than zero, indicating a minimal increase in 

concentration is likely at these 25 production wells during the 

7-year simulation period. This minimal concentration increase 

indicates a small potential for movement of the brackish-water 

transition zone near these wells during the 7-year simulation 

period simulating a hypothetical recurrence of drought-ofrecord (1950–56) rainfall conditions with higher-than-average 

groundwater withdrawals by wells.

Predictive uncertainty results for total spring discharge 

at Comal and San Marcos Springs during the 7-year period, 

as well as the results for head predictions at well J-17, were 

substantially different than the dissolved-solids concentration change uncertainty results at the production wells. The 

upper bounds of the 95-percent credible intervals for spring 

discharge at Comal and San Marcos Springs and heads at 

well J-17 were an order of magnitude larger than the actual 

predictions (based on best-fit parameters) which implies that 

(1) these predictions made with this model are not reliable and 

(2) parameters that control these predictions are not informed 

well by the observation dataset during history matching, 

even though the history-matching effort yielded parameters 

to reproduce spring discharges and heads at these locations 

during the history-matching period. Furthermore, because 

discharges at Comal and San Marcos Springs and heads at 

well J-17 represent more of a cumulative effect of upstream 

conditions over a larger distance (and longer time), many more 

parameters (with their own uncertainties) are potentially controlling these predictions than the prediction of dissolved-solids concentration change at the prediction wells, and therefore 

contributing to a large posterior uncertainty. 

References Cited

Abbott, P.L., and Woodruff, C.M., Jr., 1986, eds., The Balcones escarpment—Geology, hydrology, ecology and 

social development in central Texas: Geological Society of 

America. 

Anderson, M.P., and Woessner, W.W., 1992, Applied groundwater modeling—Simulation of flow and advective transport: San Diego, California, Academic Press, Inc., 381 p.

Barker, R.A., and Ardis, A.F., 1996, Hydrogeologic framework 

of the Edwards-Trinity aquifer system, west-central Texas: 

U.S. Geological Survey Professional Paper 1421–B, 61 p.

Blondes, M.S., Gans, K.D., Thordsen, J.J., Reidy, M.E., 

Thomas, B., Engle, M.A., Kharaka, Y.K., and Rowan, E.L., 

2013, U.S. Geological Survey National Produced Waters 

Geochemical Database v. 2.1 (provisional): accessed on 

April 9, 2013, at http://energy.usgs.gov/EnvironmentalAspects/EnvironmentalAspectsofEnergyProductionandUse/

ProducedWaters.aspx#3822349-data.

Box, G.E.P., and Draper, N.R., 1987, Empirical modelbuilding and response surfaces: Oxford, England, John 

Wiley & Sons, p. 424.

Collins, E.W., and Hovorka, S.D., 1997, Structure map of the 

San Antonio segment of the Edwards aquifer and Balcones 

fault zone, south-central Texas—Structural framework of a 

major limestone aquifer—Kinney, Uvalde, Medina, Bexar, 

Comal and Hays Counties: Austin, University of Texas, 

Bureau of Economic Geology Miscellaneous Map 38.

Crain, E.R., 2013, Crain’s petrophysical handbook: accessed 

on May 23, 2013, at http://www.spec2000.net/01-index.htm.

Dausman, A.M., Doherty, J.E., Langevin, C.D., and Dixon, 

Joann, 2010a, Hypothesis testing of buoyant plume migration using a highly parameterized variable-density groundwater model at a site in Florida, USA: Hydrogeology 

Journal, v. 18, no. 1, p. 147–160.
50 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

Dausman, A.M., Doherty, J.E., Langevin, C.D., and Sukop, 

M.C., 2010b, Quantifying data worth toward reducing predictive uncertainty: Groundwater, v. 48, no. 5, p. 729–740.

Deike, R.G., 1990, Dolomite dissolution rates and possible 

Holocene dedolomitization of water-bearing units in the 

Edwards aquifer, south-central Texas: Journal of Hydrology, 

v. 112, p. 335–373.

Doherty, J.E., 2003, Groundwater model calibration using 

pilot points and regularization: Groundwater, v. 41, no. 2, 

p. 170–177, accessed January 13, 2015, at http://info.ngwa.

org/gwol/pdf/031075995.pdf?origin=publication_detail.

Doherty, J.E., 2005, PEST—Model-independent parameter 

estimation user manual (5th ed.): Brisbane, Australia, 

Watermark Numerical Computing, 336 p.

Doherty, J.E., 2011, PEST groundwater data utilities, Part A—

Overview: Watermark Numerical Computing, 68 p.

Doherty, J.E., 2015, Calibration and uncertainty analysis for 

complex environmental models—PEST: Complete theory 

and what it means for modelling the real world: Brisbane, 

Australia, Watermark Numerical Computing, 236 p.

Doherty, J.E., and Hunt, R.J., 2010, Approaches to highly 

parameterized inversion—A guide to using PEST for 

groundwater-model calibration: U.S. Geological Survey 

Scientific Investigations Report 2010–5169, 59 p.

Doherty, J.E., and Welter, D.E., 2010, A short exploration of 

structural noise: Water Resources Research, v. 46, no. 5, 

14 p. 

Doherty, J.E., Fienen, M.N., and Hunt, R.J., 2010a, 

Approaches to highly parameterized inversion—Pilot-point 

theory, guidelines, and research directions: U.S. Geological 

Survey Scientific Investigations Report 2010–5168, 36 p.

Doherty, J.E., Hunt, R.J., and Tonkin, M.J., 2010b, 

Approaches to highly parameterized inversion—A guide to 

using PEST for model-parameter and predictive-uncertainty 

analysis: U.S. Geological Survey Scientific Investigations 

Report 2010–5211, 71 p. 

Edwards Aquifer Authority, 2015, Bexar index well (J–17) 

historic data: accessed February 3, 2015, at http://www.

edwardsaquifer.org/aquifer-data-and-maps/historical-data/

historic-data-downloads. 

Esri, 2015a, ArcGIS Desktop—Release 10: Redlands, Calif. 

[Also available at http://www.esri.com.]

Esri, 2015b, Natural neighbor (spatial analyst): accessed 

February 9, 2015, at http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z0000006p000000.htm.

Evenick, J.C., 2008, Introduction to well logs and subsurface 

maps: Tulsa, Oklahoma, PennWell Corporation, 236 p.

Fetter, C.W., 1999, Contaminant hydrogeology (2d ed.): Long 

Grove, Illinois, Waveland Press, Inc., 500 p.

Fienen, M.N., D’Oria, Marco, Doherty, J.E., and Hunt, R.J., 

2013, Approaches in highly parameterized inversion—

bgaPEST, a Bayesian geostatistical approach implementation with PEST—Documentation and instructions: U.S. 

Geological Survey Techniques and Methods, book 7, 

chapter C9, 86 p.

Fienen, M.N., Doherty, J.E., Hunt, R.J., and Reeves, H.W., 

2010, Using prediction uncertainty analysis to design hydrologic monitoring networks—Example applications from the 

Great Lakes water availability pilot project: U.S. Geological 

Survey Scientific Investigations Report 2010–5159, 44 p.

Franke, O.L., Reilly, T.E., and Bennett, G.D., 1987, Definition of boundary and initial conditions in the analysis of 

saturated ground-water flow systems—An introduction: 

Techniques of Water-Resources Investigations of the United 

States Geological Survey, book 3, chapter B5, 15 p.

Freeze, R.A., and Cherry, J.A., 1979, Groundwater: Englewood Cliffs, New Jersey, Prentice-Hall, 604 p.

Garza, Sergio, 1962, Chemical analyses of water from observation wells in the Edwards and associated limestones, San 

Antonio area, Texas: San Antonio, Tex., Edwards Underground Water District Bulletin, 56 p.

Garza, Sergio, 1966, Ground-water resources of the San Antonio area, Texas—A progress report on studies, 1960-64: 

Texas Water Development Board Report 34, 36 p.

Gelman, Andrew, Carlin, J.B., Stern, H.S., Dunson, D.B., 

Vehtari, Aki, and Rubin, D.B., 2014, Bayesian data analysis 

(3d ed.): London, CRC Press, Taylor and Francis Group, 

675 p. 

Groschen, G.E., 1994, Analysis of data from test-well sites 

along the downdip limit of freshwater in the Edwards aquifer, San Antonio, Texas, 1985–87: U.S. Geological Survey 

Water-Resources Investigations Report 93–4100, 92 p.

Groschen, G.E., and Buszka, P.M., 1997, Hydrogeologic 

framework and geochemistry of the Edwards aquifer salinewater zone, south-central Texas: U.S. Geological Survey 

Water-Resources Investigations Report 97–4133, 47 p.

Guo, Weixing, and Langevin, C.D., 2002, User’s guide to 

SEAWAT—A computer program for simulation of threedimensional variable-density ground-water flow: Techniques of Water-Resources Investigations of the United 

States Geological Survey, book 6, chap. A7, 77 p. [Supersedes U.S. Geological Survey Open-File Report 01–434.]

Harbaugh, A.W., 2005, MODFLOW-2005, the U.S. Geological Survey modular ground-water model—The groundwater flow process: U.S. Geological Survey Techniques and 

Methods, book 6, chap. A16, [variously paged].
References Cited 51

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, 

M.G., 2000, MODFLOW-2000, the U.S. Geological Survey 

modular ground-water model—User guide to modularization concepts and the ground-water flow process: U.S. 

Geological Survey Open-File Report 00–92, 121 p.

Hill, M.C., and Tiedeman, C.R., 2007, Effective groundwater 

model calibration—with analysis of data, sensitivities, predictions, and uncertainty: Hoboken, New Jersey, John Wiley 

& Sons, Inc., 455 p.

Hovorka, S.D., Dutton, A.R., Ruppel, S.C., and Yeh, J.S., 

1996, Edwards aquifer ground-water resources—Geologic 

controls on porosity development in platform carbonates, 

South Texas: Austin, The University of Texas, Bureaus of 

Economic Geology Report of Investigations 238, 75 p.

Hovorka, S.D., Mace, R.E., and Collins, E.W., 1998, Permeability structure of the Edwards aquifer, south Texas—Implications for aquifer management: Austin, The University of 

Texas, Bureau of Economic Geology Report of Investigations 250, 55 p.

Hovorka, S.D., Ruppel, S.C., Dutton, A.R., and Yeh, J.S., 

1993, Edwards aquifer storage assessment, Kinney County 

to Hays County, Texas: Austin, Texas, The University of 

Texas, Bureau of Economic Geology, 115 p.

Isaaks, E.H., and Srivastava, R.M., 1989, An introduction to 

applied geostatistics: New York, Oxford University Press, 

561 p.

Klemt, W.B., Knowles, T.R., Edler, G.R., and Sieh, T.W., 

1979, Ground-water resources and model applications for 

the Edwards (Balcones fault zone) aquifer in the San Antonio region: Texas Water Development Board Report 239, 

88 p.

Konikow, L.F., and Bredehoeft, J.D., 1992, Ground-water 

models cannot be validated: Advances in Water Resources, 

v. 15, p. 75–83.

Kuniansky, E.L., Fahlquist, Lynne, and Ardis, A.F., 2001, 

Travel times along selected flow paths of the Edwards 

Aquifer, Central Texas, in Kuniansky, E.L., U.S. Geological 

Survey Karst Interest Group—Proceedings, St. Petersburg, 

Florida, February 13–16, 2001: U.S. Geological Survey Water-Resources Investigations Report 2001–4011, 

p. 69–77.

Lambert, R.B., Hunt, A.G., Stanton, G.P., and Nyman, M.B., 

2009, Water-level, borehole geophysical log, and waterquality data from wells transecting the freshwater/salinewater interface of the San Antonio segment of the Edwards 

aquifer, south-central Texas, 1999–2007: U.S. Geological 

Survey Data Series 403, [variously paged]. [Also available 

at http://pubs.usgs.gov/ds/403/.]

Lambert, R.B., Hunt, A.G., Stanton, G.P., and Nyman, M.B., 

2010, Lithologic and physicochemical properties and 

hydraulics of flow in and near the freshwater/saline-water 

transition zone, San Antonio segment of the Edwards aquifer, south-central Texas, based on water-level and borehole 

geophysical log data, 1999–2007: U.S. Geological Survey 

Scientific Investigations Report 2010–5122, 69 p.

Langevin, C.D., 2000, Simulation of ground-water discharge 

to Biscayne Bay, southeastern Florida: U.S. Geological 

Survey Water-Resources Investigations Report 2000–4251, 

127 p.

Langevin, C.D., Shoemaker, W.B., and Guo, Weixing, 2003, 

MODFLOW-2000, the U.S. Geological Survey Modular Ground-Water Model—Documentation of the SEAWAT-2000 version with the Variable-Density Flow Process 

(VDF) and the Integrated MT3DMS Transport Process 

(IMT): U.S. Geological Survey Open-File Report 03–426, 

43 p.

Langevin, C.D., Thorne, D.T., Jr., Dausman, A.M., Sukop, 

M.C., and Guo, Weixing, 2008, SEAWAT Version 4—A 

computer program for simulation of multi-species solute 

and heat transport: U.S. Geological Survey Techniques and 

Methods, book 6, chap. A22, 39 p.

Larkin, T.J., and Bomar, G.W., 1983, Climatic atlas of Texas: 

Austin, Texas Department of Water Resources, LP–192, 

p. 157, accessed January 12, 2015, at http://www.twdb.state.

tx.us/publications/reports/limited_printing/doc/LP192.pdf.

LBG-Guyton Associates, 1995, Edwards aquifer ground-water 

divides assessment, San Antonio region, Texas: San Antonio, Edwards Underground Water District Report 95–01, 

35 p.

Lindgren, R.J., 2006, Diffuse-flow conceptualization and 

simulation of the Edwards aquifer, San Antonio region, 

Texas: U.S. Geological Survey Scientific Investigations 

Report 2006–5319, 48 p.

Lindgren, R.J., Dutton, A.R., Hovorka, S.D., Worthington, 

S.R.H., and Painter, Scott, 2004, Conceptualization and 

simulation of the Edwards aquifer, San Antonio region, 

Texas: U.S. Geological Survey Scientific Investigations 

Report 2004–5277, 143 p.

Lindgren, R.J., Houston, N.A., Musgrove, M., Fahlquist, L.S., 

and Kauffman, L.J., 2011a, Simulations of groundwater 

flow and particle-tracking analysis in the zone of contribution to a public-supply well in San Antonio, Texas: U.S. 

Geological Survey Scientific Investigations Report 2011–

5149, 93 p.
52 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

Lindgren, R.J., Houston, N.A., Musgrove, M., Fahlquist, L.S., 

and Kauffman, L.J., 2011b, Hydrogeologic setting and 

groundwater-flow simulations of the South-Central Texas 

regional study area, Texas, section 3 of Eberts, Sandra M., 

ed., Hydrogeologic settings and groundwater-flow simulations for regional investigations of the transport of anthropogenic and natural contaminants to public-supply wells—

Investigations begun in 2004: U.S. Geological Survey 

Professional Paper 1737-B, p. 3-1–3-51.

Lindgren, R.J., Taylor, C.J., and Houston, N.A., 2009, 

Description and evaluation of numerical groundwater flow 

models for the Edwards aquifer, south-central Texas: U.S. 

Geological Survey Scientific Investigations Report 2009–

5183, 25 p.

Lozo, F.E., Jr., and Smith, C.I., 1964, Revision of Comanche 

Cretaceous stratigraphic nomenclature, southern Edwards 

Plateau, southwest Texas: Gulf Coast Association of Geological Societies Transactions, v. 14, p. 285–306.

Maclay, R.W., 1995, Geology and hydrology of the Edwards 

aquifer in the San Antonio area, Texas: U.S. Geological Survey Water-Resources Investigation Report 95–4186, 64 p.

Maclay, R.W., and Land, L.F., 1988, Simulation of flow in the 

Edwards aquifer, San Antonio region, Texas and refinements of storage and flow concepts: U.S. Geological Survey 

Water-Supply Paper 2336–A, 48 p.

Maclay, R.W., and Small, T.A., 1984, Carbonate geology and 

hydrology of the Edwards aquifer in the San Antonio area, 

Texas: U.S. Geological Survey Open-File Report 83–537, 

72 p.

Meyer, C.D., 2000, Matrix analysis and applied linear algebra: Philadelphia, Pennsylvania, Society for Industrial and 

Applied Mathematics, 718 p.

Miller, L.D., and Long, A.J., 2006, Statistical analyses of 

hydrologic system components and simulation of Edwards 

aquifer water-level response to rainfall using transfer-function models, San Antonio region, Texas: U.S. Geological 

Survey Scientific Investigations Report 2006–5131, 20 p. 

[1 appendix, online only].

Moore, Catherine, and Doherty, J.E., 2005, Role of the calibration process in reducing model predictive error: Water 

Resources Research, v. 41, no. 5, 14 p.

Musgrove, MaryLynn, Fahlquist, Lynne, Houston, N.A., Lindgren, R.J., and Ging, P.B., 2010, Geochemical evolution 

processes and water-quality observations based on results 

of the National Water-Quality Assessment Program in the 

San Antonio segment of the Edwards aquifer, 1996–2006: 

U.S. Geological Survey Scientific Investigations Report 

2010–5129, 93 p.

Oetting, G.C., Banner, J.L., and Sharp, J.M., Jr., 1996, 

Regional controls on the geochemical evolution of saline 

groundwaters in the Edwards aquifer, central Texas: Journal 

of Hydrology, v. 181, p. 251–283.

Oliver, D.S., Reynolds, A.C., and Liu, Ning, 2008, Inverse 

theory for petroleum reservoir characterization and history 

matching: New York, Cambridge University Press, 380 p.

Painter, Scott, Jiang, Yefang, and Woodbury, Allan, 2002, 

Edwards aquifer parameter estimation project final report: 

Southwest Research Institute [variously paged].

Pavlicek, Diane, Small, T.A., and Rettman, P.L., 1987, Hydrogeologic data from a study of the freshwater zone/salinewater zone interface in the Edwards aquifer, San Antonio 

region, Texas: U.S. Geological Survey Open-File Report 

87–389, 108 p.

Perez, Roberto, 1986, Potential for updip movement of saline 

water in the Edwards aquifer, San Antonio, Texas: U.S. 

Geological Survey Water-Resources Investigations Report 

86–4032, 21 p.

Puente, Celso, 1978, Method of estimating natural recharge 

to the Edwards aquifer in the San Antonio area, Texas: U.S. 

Geological Survey Water-Resources Investigations Report 

78–10, 34 p.

Railroad Commission of Texas, 2012, Online research queries—Imaged records: accessed August 22, 2012, at http://

www.rrc.state.tx.us/about-us/resource-center/research/

online-research-queries/about-imaged-records/.

Refsgaard, J.C., van der Sluijs, J.P., Højberg, A.L., and Vanrolleghem, P.A., 2007, Uncertainty in the environmental modelling process—A framework and guidance: Environmental 

Modelling & Software, v. 22, p. 1543–1556.

Rose, P.R., 1972, Edwards Group, surface and subsurface, 

central Texas: Austin, The University of Texas, Bureau of 

Economic Geology Report of Investigations 74, 198 p.

San Antonio Water System, 2012, Water statistics, year ending 

December 31, 2012: Infrastructure Planning Department, 

San Antonio Water System, accessed December 22, 2014, at 

http://www.saws.org/your_water/aquifer/docs/2012_StatsBook.pdf.

Schlumberger, 2015, Oilfield glossary: accessed February 17, 

2015, at http://glossary.oilfield.slb.com/.

Schultz, A.L., 1992, Using geophysical logs in the Edwards 

aquifer to estimate water quality along the freshwater/

saline-water interface (Uvalde to San Antonio, Texas): 

Edwards Underground Water District Report 92–03, 47 p.
References Cited 53

Schultz, A.L., 1993, Defining the Edwards aquifer freshwater/

saline-water interface with geophysical logs and measured 

data (San Antonio to Kyle, Texas): Edwards Underground 

Water District Report 93–06, 81 p.

Schultz, A.L., 1994, Review and update of the position of the 

Edwards aquifer freshwater/saline-water interface from 

Uvalde to Kyle, Texas: Edwards Underground Water District Report 94–05, 31 p.

Slattery, R.N., 2004, Recharge to the Edwards aquifer in the 

San Antonio area, Texas, 2003: accessed July 12, 2004, at 

http://tx.usgs.gov/reports/dist/dist-2004-01.

Texas Commission on Environmental Quality, 2013, Public 

drinking water, subchapter F—Drinking water standards 

governing drinking water quality and reporting requirements for public water systems, summary of secondary 

standards, chap. 290, section 105, accessed January 14, 

2015, at http://www.tceq.state.tx.us/assets/public/legal/

rules/rules/pdflib/290f.pdf.

Texas Water Development Board, 2004, Historical water-use 

information: accessed May 29, 2015, at https://www.twdb.

texas.gov/waterplanning/waterusesurvey/estimates/.

Texas Water Development Board, 2013, Groundwater database: accessed August 14, 2013 at http://www.twdb.texas.

gov/groundwater/data/Database%20in%20ASCII/All%20

Counties/gwdb.zip.

Texas Water Development Board, 2014, Water data for 

Texas—Drought in Texas: accessed October 29, 2014, at 

http://waterdatafortexas.org/drought/. 

Thomas, J.V., Stanton, G.P., and Lambert, R.B., 2012, Borehole geophysical, fluid, and hydraulic properties within and 

surrounding the freshwater/saline-water transition zone, 

San Antonio segment of the Edwards aquifer, south-central 

Texas, 2010–11: U.S. Geological Survey Scientific Investigations Report 2012–5285, 65 p., 3 apps.

Thorkildsen, D.F., and McElhaney, P.D., 1992, Model refinement and applications for the Edwards (Balcones fault 

zone) aquifer in the San Antonio region, Texas: Texas Water 

Development Board Report 340, 33 p.

U.S. Census Bureau, 2014, Annual estimates of the resident 

population for incorporated places of 50,000 or more, 

ranked by July 1, 2013 population, April 1, 2010 to July 1, 

2013: U.S. Census Bureau, Population Division, accessed 

February 18, 2015, at http://factfinder.census.gov/faces/

tableservices/jsf/pages/productview.xhtml?pid=PEP_2013_

PEPANNRSIP.US12A&prodType=table.

U.S. Environmental Protection Agency, 1975, Edwards 

underground reservoir—Notice of determination: Federal 

Register, v. 40, no. 242, p. 58344–58345.

U.S. Environmental Protection Agency, 1988, A portion of the 

Austin-area Edwards aquifer in parts of Hays and Travis 

Counties, Texas—Sole source aquifer, final determination: 

Federal Register, v. 53, no. 109, p. 20897–20899.

U.S. Environmental Protection Agency, 2015, Sole source 

aquifers: accessed February 18, 2015, at http://www.epa.

gov/region6/water/swp/ssa/index.htm. 

U.S. Geological Survey, 2013, National Water Information 

System—Web interface (NWISWeb) data: accessed August 

2013 at http://waterdata.usgs.gov/tx/nwis/nwis.

U.S. Geological Survey, 2014, SEAWAT—A computer program for simulation of three-dimensional variable-density 

groundwater flow and transport: accessed August 6, 2014, at 

http://water.usgs.gov/ogw/seawat/. 

Voss, C.I., 2011a, Editor’s message, Groundwater modeling 

fantasies—Part 1, adrift in the details: Hydrogeology Journal, v. 19, no. 7, p. 1281–1284. 

Voss, C.I., 2011b, Editor’s message, Groundwater modeling 

fantasies—Part 2, down to earth: Hydrogeology Journal, 

v. 19, no. 8, p. 1455–1458. 

Waite, L.E., 2008, Edwards (Stuart City) shelf margin of 

south Texas—New data, new concepts: AAPGDatapages, 

Search and Discovery Article no. 10177 (2009), accessed on 

December 9, 2014, at http://www.searchanddiscovery.com/

documents/2009/10177waite/. 

Welter, D.E., Doherty, J.E., Hunt, R.J., Muffels, C.T., Tonkin, 

M.J., and Schreüder, W.A., 2012, Approaches in highly 

parameterized inversion—PEST++, a Parameter ESTimation code optimized for large environmental models: U.S. 

Geological Survey Techniques and Methods, book 7, 

sec. C5, 47 p. 

White, J.T., Doherty, J.E., and Hughes, J.D., 2014, Quantifying the predictive consequences of model error with linear 

subspace analysis: Water Resources Research, v. 50, 

p. 1152–1173.

White, J.T., Welter, D.E., Hunt, R.J., and Doherty, J.E., 2015, 

Integrating linear-based uncertainty analyses into PEST++, 

in Maxwell, Reed, Hill, Mary, Zheng, Chunmiao, and 

Tonkin, Matt, MODFLOW and more 2015: Modeling a 

complex world—Proceedings, Golden, Colorado, May 31–

June 3, 2015, 559 p. 

Winston, R.B., 2015, Online guide to MODFLOW-2005: 

accessed February 17, 2015, at http://water.usgs.gov/ogw/

modflow/MODFLOW-2005-Guide/.

Winters, K.E., 2013, A historical perspective on precipitation, 

drought severity, and streamflow in Texas during 1951–56 

and 2011: U.S. Geological Survey Scientific Investigations 

Report 2013–5113, 24 p. 
54 Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement

Wong, C.I., Kroman, J.S., Hunt, B.B., Smith, B.A., and Banner, J.L., 2013, Investigating groundwater flow between 

Edwards and Trinity aquifers in Central Texas: Groundwater, v. 52, no. 4, p. 624–639.

Zheng, Chunmiao, and Bennett, G.D., 2002, Applied contaminant transport modeling: New York, Wiley-Interscience, 

621 p.

Zheng, Chunmiao, and Wang, P.P., 1999, MT3DMS, A modular three-dimensional multi-species transport model 

for simulation of advection, dispersion and chemical reactions of contaminants in groundwater system—Documentation and user’s guide: Vicksburg, Mississippi, U.S. Army 

Engineer Research and Development Center Contract 

Report SERDP-99-1, 202 p.

Publishing support provided by

Lafayette and Rolla Publishing Service Centers

Information regarding water resources in Texas is available at

http://tx.usgs.gov/

Brakefied and others—Updated Numerical Model with Uncertainty Assessment of 1950–56 Drought Conditions on Brackish-Water Movement—Scientific Investigations Report 2015–5081

ISSN 2328-0328 (online)

http://dx.doi.org/10.3133/sir20155081