Estimate of Groundwater Flow and Salinity Contribution to Great Salt Lake Using Groundwater Levels and Spatial Analysis

Groundwater discharge to Great Salt Lake (GSL) is difficult to quantify but represents a potentially significant source of water and salinity to the lake’s overall water budget and chemistry, respectively. Understanding groundwater and its role in the overall health of GSL is critical due to the current historically low lake levels. We compiled existing groundwater level data in wells in the basin-fill aquifer around GSL and used spatial analysis methods to 1) create potentiometric-surface maps in the areas adjoining GSL, 2) calculate groundwater contributions to GSL, and 3) estimate salinity inputs from groundwater to GSL. We observed groundwater-level declines in most of the basin-fill wells from the 1980s to 2010s. These declines are consistent with historical groundwater-level trends in the Salt Lake, Tooele, Curlew, and Weber Valleys and are a consequence of aquifer overdraft associated with less than average precipitation in the basin and increased groundwater withdrawals in the GSL watershed. Using the Darcy flux equation, we calculated a groundwater flux to GSL of 313,500 acre-feet per year, substantially greater than previous estimates derived from water balance studies but consistent with estimates derived from geochemical modeling of GSL water chemistry. We calculated a salt contribution from groundwater to GSL of 1.18 million metric tons per year, which represents about 10% of the solutes derived from surface flows to GSL in 2013.


Background and Objectives
Groundwater discharge is an essential component in limnological systems' hydrologic and chemical balances (Healy and others, 2007;Rosenberry and Winter, 2009).Despite its importance in developing accurate hydrochemical balances in lakes, groundwater contribution is often neglected or underestimated because it is difficult to quantify (Rosenberry and others, 2015).Understanding the groundwater component and its associated solute input is imperative for managing the environmental and economic resources of lakes affected by extensive anthropogenic water use and drought, such as Great Salt Lake (GSL) in northern Utah.Potentiometric gradient and aquifer hydraulic conductivity are key inputs to calculate groundwater flow.This report focuses on the potentiometric data input for estimating groundwater flow.Future work will better constrain hydraulic conductivity and geochemistry at the lake interface to improve groundwater flow estimates.
GSL is a hypersaline terminal lake and a sink for surface and groundwater across a large part of the eastern Great Basin (Spencer and others, 1985;Duffy and Al-Hassan, 1988;Arnow and Stephens, 1990).
Salinity inputs and evaporation impact GSL's ecosystems and mineral resources (Carling and others, 2013;Jagniecki and others, 2021).Surface-water flows to GSL and associated salt loading are well constrained (Shope and Angeroth, 2015).However, the quantity of groundwater discharge remains relatively unknown, and groundwater is a potentially significant source to GSL's overall salt load (Kirby and others, 2019;Bunce, 2022).The importance of understanding groundwater dynamics and its role in the overall health of GSL becomes prominant by the current, historically low lake levels.Groundwater inflows to GSL will become critically important as surfacewater discharges decrease due to increasing water demands (Null and Wurtsbaugh, 2020), rising air temperatures, and changing snow cover conditions in the basin (Hall and others, 2021).
Given the significance and uncertainty of groundwater to GSL's system, the objectives of the present study are to 1) compile historical groundwater levels and use them to create generalized potentiometricsurface maps, 2) roughly estimate groundwater flow to GSL using a combination of spatial analysis techniques and Darcy's Law, and 3) combine the results from the previous objective with existing groundwater chemistry data to estimate salinity inputs to GSL derived from groundwater.Water and salt dynamics play a fundamental role in shaping not only GSL's unique ecological, recreational, and mineral resources, but also the future development in the sprawling urban centers on the east shore of the lake.We provide the first systematic basin-wide assessment of groundwater levels in areas adjoining GSL to quantify groundwater contributions and their salt loading to the lake system.This information is required to constrain water and salt budgets needed by managing agencies to make informed decisions regarding the future health and productivity of ecosystems and industries in the lake.Our findings serve as a basis for future work to better define the role of groundwater in GSL's overall water volume and solutes budgets.

Study Area
The study area covers about 11,000 square miles in northern Utah and southern Idaho (figure 1).It is bounded on the east by the north-south-trending Wasatch Range and on the west by the Great Salt Lake Desert.Several mountain ranges (Hansel, Promontory, Oquirrh, Stansbury, Hogup, etc.) and associated valleys (Curlew, Hansel, Malald-Lower Bear River, Weber, Salt Lake, Tooele, Skull, etc.) bound the north and south flanks of the study area and drain into GSL (figure 1).GSL is the largest salt-water lake by area in the Western Hemisphere and the eighth largest in the world (Hammer, 1986).It is 75 miles long by 28 miles wide and covers approximately 1700 square miles with a maximum depth of 33 feet at the average water-surface elevation of 4200 feet above sea level (FASL).The Jordan, Bear, and Weber Rivers deliver on average 2.9 million acre-feet of water to GSL, approximately 95% of the total stream inflow (Stephens and Gardner, 2007;Mohammed and Tarboton, 2012).Previous studies (Arnow and Stephens, 1975;Waddell and Fields, 1976;Loving and others, 2000;Bunce, 2022) have estimated groundwater discharge to GSL ranges between 3% and 10% of the total inflow.The Jordan, Bear, and Weber Rivers also delivered an estimated 14.3 million metric tons of total dissolved solids (TDS) in 2013 (Shope and Angeroth, 2015).Groundwater potentially contributes a significant input to GSL's overall salt load (Hahl, 1968;Spencer and others, 1985;Loving and others, 2000), but the relationship between groundwater and GSL's salinity has not been well defined.
Salinity of the lake ranges from 5% to 29% and creates diverse opportunities for ecological, recreational, and mineral uses.GSL is part of the Pacific Flyway and provides important nesting and foraging habitat for over 250 species of birds as they travel between North and South America (The Nature Conservancy, 2022).Between 1.6 and 2.5 million metric tons of salt are commercially removed from the lake every year (Stephens and Gardner, 2007).Mineral extraction, brine shrimp cyst production, and recreation in GSL can generate an estimated economic value of $1.32 billion per year (Bioeconomics, 2012).
As a terminal lake, GSL loses water primarily through evaporation.Therefore, changes in streamflow conditions severely impact lake levels and salinity (Mohammed and Tarboton, 2012).Stream diversions for agricultural and municipal uses reduce the amount of water flowing into GSL by 39% (Null and Wurtsbaugh, 2020).These diversions and a warming climate led to the present-day lake-level decline (Wurtsbaugh and others, 2017;Wang and others, 2018).As of November 2022, the area covered by GSL was reduced to about 900 square miles at a historical low water-surface elevation of 4188.5 FASL.Consequently, salinity and the surface area covered by dry lakebed increased.Increased salinity levels stress microbialite, brine fly, and brine shrimp populations, jeopardizing the entire ecological community that depends on them.Dry lake beds are a major source of dust pollution and have the potential to accelerate snowmelt when dust is blown onto the snow (Reynolds and others, 2014;Skiles and others, 2018).

Data Compilation
We compiled historical groundwater level data from various datasets including the U.S. Geological Survey (USGS) National Water Information System (NWIS, https://waterdata.usgs.gov/,U.S. Geological Survey, 2022a), the Utah Division of Water Rights (DWRi) well drilling records, the Utah Geological Survey (UGS) Geologic Hazards Program Subsurface Geotechnical Database, and the UGS Wetlands section field data.Data from the DWRi was derived from a combination of well information tables and the Water Rights points of diversion (WRPOD) feature class provided on the DWRi website (https:// waterrights.utah.gov/).The WRPOD feature class only includes wells with a Well Identification Number (WIN) and excludes wells without an assigned WIN.The UGS Geologic Hazards Program Subsurface Geotechnical Database consists of 5141 boreholes in watersheds that contribute to the lake.The UGS Wetlands section field data consists of 362 shallow boreholes in the wetlands proximal to GSL.In addition to groundwater level data, well properties (latitude, longitude, surface elevation, screen depth, and well depth) were also gathered where available.Data compilation is limited to wells in the basin-fill aquifer and boreholes within the study area (figure 2).The compiled groundwater level data show spatial and temporal complexities within the study area.Groundwater level measurements began in the 1900s but were limited to only a few sites (33).The number of data collection sites increased over time to reach a maximum of 3136 during the 1990s but remained variable (figure 3).Furthermore, while many sites were visited regularly, others were visited less often, or only once in many cases.We reprojected, combined, and aggregated the data using Python 3 and created two different datasets to deal with these spatial and temporal intricacies.The two datasets (tables S1 and  S2) are included as supplementary information in geodatabase format (contact authors for database).
The first dataset (table S1) contains data available from all sources (NWIS, DWRi, and UGS).Groundwater level data for each site were grouped by decade and a mean water level elevation value was calculated if more than a single measurement was recorded during the time frame.The second dataset (table S2) contains only USGS NWIS data for all wells available within the study area.Here, the groundwater elevation value is the mean of all available data for each site.This dataset also contains estimated aquifer properties including saturated thickness, hydraulic conductivity, and transmissivity for individual sites.Saturated thickness values were estimated by subtracting mean depth to water values from total well depth values.Hydraulic conductivity values were extracted from layer 2 of the USGS's groundwater model of the Great Basin carbonate and alluvial aquifer (Brooks, 2017).Finally, transmissivity values were calculated by multiplying saturated thickness times hydraulic conductivity.
We downloaded three high-resolution aerial photographs from the European Space Agency's Sentinel -2 satellites for August 25th, August 28th, and September 9th, 2022, using the USGS Global Visualization Viewer (GloVis, https://glovis.usgs.gov/).Sentinel-2 satellites carry an optical payload with visible, near infrared and shortwave infrared sensors encompassing 13 spectral bands: 4 bands at 10-meter, 6 bands at 20-meter, and 3 bands at 60-meter spatial resolution, with a swath width of 290 kilometers (European Space Agency, 2022).We created a single mosaicked image in ArcGIS Pro® and used it to trace the extent of the lake.We assigned a surface elevation of 4190 feet to this extent based on USGS lake stage observations at the Saltair boat harbor (USGS Site ID 10010000, https://webapps.usgs.gov/gsl/) on August 1st, 2022 (refer to figure 1 for location).
Lidar datasets are available for northern Utah but do not cover the entire study area.For consistency, we downloaded six digital elevation model (DEM) tiles from the Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Version 3 using the USGS Earth Explorer (https:// earthexplorer.usgs.gov/).These DEM tiles have a one -third arc-second resolution (~10 meters) and are referenced to the North American Vertical Datum of 1988 (NAVD 88).We created a single mosaicked image in ArcGIS Pro® and used it to extract surface elevations in well sites where this information was missing.We also used the ASTER DEM as an explanatory variable input in the data interpolation process (discussed below).

Potentiometric Surface Interpolation
We used the Empirical Bayesian Kriging Regression Prediction (EBK-RP) tool in ArcGIS Pro® to interpolate the water level elevation values in table S1.EBK-RP is a geostatistical interpolation method that uses EBK with explanatory variable rasters known to affect the value of the data being interpolated.The tool combines kriging with regression analysis to make predictions that are more accurate than either kriging or regression can achieve on their own (ESRI, 2022a).As the potentiometric surface generally follows topography, and because it is a component in the calculation of groundwater level elevation, we included the ASTER DEM as an explanatory variable.We created two generalized potentiometric-surface maps: one for the 1980s (figure 4) and one for the 2010s (figure 5).These two decades were chosen because 1) not enough information is available to create a potentiometric surface for the 2020s, thus the 2010s data represent the most recent groundwater conditions for the study area, 2) average decadal GSL surface levels were close to the historical average (4200 FASL) in the 1980s, and 3) both decades contain about the same number of sites available for interpolation (figure 3).We also used table S1 to create a water   level difference map by subtracting water level elevation values (figure 6).This calculation was possible only for sites where data was available for both decades.
We used the Kriging tool in ArcGIS Pro® to interpolate the water elevation, saturated thickness, and transmissivity values in Table S2 (Figure 7).The three output rasters are floating point, have the same areal extent, and have a cell size of 100.Units are consistent for time (day) and length (feet) for all data.We used these rasters as input for the Darcy Flow tool in ArcGIS Pro® to estimate groundwater seepage velocity values around GSL (discussed below).

Groundwater Seepage Velocity and Darcy Flux
We applied two different approaches to estimate groundwater flow to Great Salt Lake.Both approaches are basic Darcy Flow estimates and should be considered rough approximations.The first approach used tools available in GIS software and the second approach used basic linear discretization of the valleys around the lake.

GIS-Based Calculations of Darcy Flow
For the first approach, we used the Darcy Flow tool in ArcGIS Pro® to estimate groundwater seepage velocity around the lake.This method uses Darcy's Law to model two-dimensional, vertically mixed, horizontal, and steady state flow, where groundwater head is independent of depth.Darcy's Law states that Darcy velocity in porous material is calculated from the hydraulic conductivity and the hydraulic gradient as: (1) where: q = Darcy velocity, Darcy flux, or specific discharge (V/T/A or L/T) K = hydraulic conductivity (L/T) ∇i = hydraulic gradient (dimensionless) ∇H = change in hydraulic head over length (L) ∇L = change in length (L) Hydraulic conductivity (K) may be calculated from the transmissivity and thickness as: (2) where: T = transmissivity (L 2 /T) b = aquifer thickness (L) The specific discharge (q) is defined as the volume of water flow per unit time through a cross-sectional ar-ea normal to the direction of flow (Bear, 1979).Specific discharge is directly proportional to hydraulic conductivity.The aquifer flux is defined as: (3) where: U = aquifer flux (V/T/L) T = transmissivity (L2/T) ∇i = hydraulic gradient (dimensionless) The aquifer flux (U) represents the discharge per unit width of the aquifer.The average fluid velocity within the pores, or seepage velocity, is the Darcy velocity (q) divided by the effective porosity of the medium: (4) where: V = groundwater seepage velocity (L/T) q = Darcy velocity, Darcy flux, or specific discharge (V/T/A or L/T) n = effective porosity (%) The groundwater seepage velocity (V) is calculated on a cell-by-cell basis in the Darcy Flow procedure (ESRI, 2022b).For cell {i, j}, the aquifer flux (U) is calculated through each of the four cell walls, using the difference in heads between the two adjacent cells and the harmonic average of the transmissivities (Konikow and Bredehoeft, 1978), which are assumed to be isotropic (ESRI, 2022b).
The Darcy Flow tool requires four raster datasets as input: groundwater head elevation (FASL), saturated thickness (feet), formation transmissivity (square feet per day), and effective formation porosity.We created the first three datasets by interpolating the aquifer property data in table S2 (figure 7), and we assigned the effective formation porosity a value of 0.35 for the basin-fill aquifer.Two raster datasets result from this calculation: an output magnitude raster and an output direction raster.In the output magnitude raster, each cell value represents the magnitude of the seepage velocity vector (average linear velocity) at the center of the cell and is calculated as the average value of the seepage velocity through the four faces of the cell (ESRI, 2022b).In the output flow direction raster, each cell value represents the direction of the seepage velocity vector (ESRI, 2022b).
We extracted the mean groundwater seepage velocity around GSL from the output magnitude raster using the Zonal Statistics as Table tool in ArcGIS Pro® with the GSL perimeter shapefile (both 4190 and 4200 FASL, table 1) as the Feature Zone Data (dataset that defines zone of interest).We conducted zonal statistics using both polylines and polygons of GSL at both elevations to compare how the cells' statistics were aggregated.We also divided the 4190 FASL GSL perimeter shapefile into three different   S2.
sections: west/south, north, and east (figure 8).We extracted the groundwater seepage velocities from these different sections (using the Zonal Statistics Tool as above) to estimate the potential groundwater contribution by area (table 2).Finally, we multiplied the groundwater seepage velocity by the mean aquifer thickness (232 feet, calculated from Table S2) and the perimeter length of GSL (obtained from the GSL perimeter shapefiles) to obtain an estimated groundwater contribution (Q) in acre-feet per year (tables 1 and 2).

Straight Line Calculations of Darcy Flow
For the second approach, we estimated the groundwater discharge from each adjoining basin-fill valley using Darcy's Law: (5) For this approximation, we used a mean K value of 12.43 feet per day (from table S2).We created several polygons to roughly constrain the areal extent of the basin fill in the valleys (figure 8).We calculated the cross-section area (A) for the flux calculations by multiplying the width of the polygon (perpendicular to the flow direction) by the mean aquifer thickness of 232 feet (from table S2).We calculated the hydraulic gradient (∇i) by extracting H1 and H2 values from the groundwater elevation raster, near the upper and lower ends of the polygons (figure 8) and dividing the difference between the two piezometric heads by the length of the polygon (parallel to groundwater flow).The geometrical properties and estimated groundwater discharge for each basin-fill area, (Q) in acre-feet per year, are presented in table 3.

Quantifying Error
To better understand the amount of variability that can be introduced by hydraulic conductivity, we iteratively calculated the groundwater flow for each area using the straight-line approach.We created a lognormal distribution of hydraulic conductivity, us-ing 1.094 (log 10 of 12.43 ft/day) as the mean and 1.3 as the standard deviation, which is the standard deviation of the hydraulic conductivity of the upper basinfill aquifer from the USGS groundwater model.We randomly sampled hydraulic conductivity 100,000 times from this distribution and used the resulting summary statistics to constrain variability of groundwater flow estimates from the calculations.

Salt Loading
We calculated mean TDS values around GSL using data available from previous studies in the area (Kirby and others, 2019).From Kirby and others' (2019) data, we obtained three different TDS values using 1) their calculated TDS values to compute a mean TDS for the whole study area, 2) their TDS raster and the Zonal Statistics as Table tool in ArcGIS Pro® to estimate the mean TDS value around GSL at 4190 FASL (polyline), and 3) their TDS raster and the Zonal Statistics as Table tool in ArcGIS Pro® to estimate the mean TDS value around GSL at 4200 FASL (polyline, table 4).We used the sectioned 4190 FASL GSL shapefile (west/south, north, and east) in figure 8 to extract mean TDS values by segment from the TDS raster using the Zonal Statistics as Table tool in ArcGIS Pro® (table 5).This approach helped to estimate the salt loading contribution by section.We combined the mean TDS values, in milligrams per liter (mg/L), with the groundwater discharge (Q, acrefeet) values in tables 1 and 2 to estimate a salt loading to the lake in metric tons per year (tables 4 and 5).

Generalized Potentiometric Surface
The generalized potentiometric surface maps for the 1980s and 2010s show that, at the scale of the study area, groundwater flow patterns are relatively constant over time (figures 4 and 5).Groundwater flows from the high-elevation mountains surrounding the study area towards the adjacent valleys and into GSL.For both decades, the steepest groundwater gradients occur in the vicinity of the Raft River Mountains, across the west flank of the Oquirrh Mountains (southeastern Tooele Valley), across the east flank of the Oquirrh Mountains (western Salt Lake Valley), and along the Wasatch Range (figures 4 and 5).There are noticeable changes in the potentiometric surface between the two decades.For example, in Salt Lake Valley, the 4300-foot contour has moved farther south (upstream), and we observed cones of depression north of the Oquirrh Mountains and in the Taylorsville and Bountiful areas (figures 9 and 10).We also observed some recovery in the Weber River delta area where groundwater recharge projects have been operating since the 2000s (Hurlow and others, 2011;figures 11 and 12). Figure 6 also shows groundwaterlevel recovery in some wells.However, groundwaterlevel recovery is an exception rather than the rule because 129 out of 147 wells having data available in both decades show a groundwater-level decline (figure 6).Groundwater-level declines are particularly high in localized areas of the north Malad River Valley in Idaho (-86 feet), Curlew Valley (-35 to -49 feet), near the Bear River (-51 feet), and Salt Lake Valley (-44 feet).The data in figure 6 show there has been a mean decrease in water levels of 10.14 feet from the 1980s to 2010s.

Seepage Velocity and Darcy Flux
Mean groundwater seepage velocities into GSL are 0.15 (polyline) and 0.12 (polygon) feet per day using the 4190 FASL shapefiles (table 1).Mean groundwater seepage velocities along the GSL perimeter are 0.09 (polyline) and 0.11 (polygon) feet per day using the 4200 FASL shapefiles (table 1).These mean seepage velocities result in groundwater fluxes ranging from 310,000 to 411,000 acre-feet per year (table 1).The average of these estimates is 356,000 acre-feet per year.The mean groundwater seepage velocities by section (4190 FASL) are 0.05 (west/ south), 0.30 (north), and 0.10 (east) feet per day (table 2).These mean seepage velocities result in groundwater fluxes of 50,291 (west/south), 274,660 (north), and 85,829 (east) acre-feet per year for a total of 411,000 acre-feet per year (table 2).
Using the Darcy flux equation on the linear traverses in adjoining valleys (figure 8), groundwater contributions to GSL range from 3500 (Skull Valley) to 63,900 (Tooele Valley) acre-feet per year for a total of 313,500 acre-feet per year (table 3).The groundwater discharge from the Wasatch Range (sum of Salt Lake, Weber, and Malad Valley and Brigham City) is 146,400 acre-feet per year, 47% of the total groundwater contribution to the lake.

Salt Loading
The mean TDS values along GSL are 2116, 2538, and 3594 mg/L with an average value of 2749 mg/L using three different methods as explained in the "Methods" section (table 4).These TDS values and the calculated groundwater contributions in tables 1 and 3 result in yearly dissolved solid (salt) fluxes ranging from 810,000 to 1,820,000 metric tons, and an average of 1,060,000 metric tons, into GSL (table 4).Table 5 reports the mean TDS values by perimeter section as 4627 (west/south), 3055 (north), and 3023 (east) mg/L.These TDS values and the calculated groundwater contribution by section in table 3 result in yearly salt fluxes of 290,000 (west/south), 1,040,000 (north), and 320,000 (east) metric tons (table 5), resulting in 1,650,000 metric tons per year.

Error
The uncertainty of the Darcy flux calculations on the linear traverses through adjoining valleys is high.The 5th percentile for flow values is 2350 acre-feet/yr and the 95th percentile is 44.5 million ac-ft/yr.See table 3 for a complete list of variations associated with potential variability in hydraulic conductivity.

Generalized Potentiometric Surface
Groundwater levels declined in most of the basinfill aquifer from the 1980s to the 2010s (figure 6).These declines are consistent with observed historical  (Burden and others, 2005).Groundwater levels declined by an average of 27 feet from 1953 to 1985 in the Ogden area (Hurlow and others, 2011).Our data shows that, in general terms, this trend continued into the 2010s.These declines are a consequence of aquifer overdraft associated with less than average precipitation in the basin and increased withdrawals for municipal, industrial, and agricultural use (Burden and others, 2005).Young and others (2021) estimated that the GSL basin lost 8.8 ± 2.3 million acre-feet of groundwater storage whereas GSL lost 4.5 ± 0.8 million acre-feet of surface water during the 2012-2016 drought.
Overdraft conditions in basin-fill aquifers can cause several problems.Groundwater removal through pumping for anthropogenic use (and associated groundwater level declines) may lead to reduction of water in streams and lakes, land subsidence or ground failure due to soil compaction, increased costs for users due to higher pumping lifts, and deterioration of water quality from saltwater intrusion (U.S. Geological Survey, 2022b).Consumptive water uses in the GSL watershed have already depleted surface inflows to the lake by ~39% (Null and Wurtsbaugh, 2020).These inflow reductions are in part responsible for the recent ~10-foot drop in GSL's surface water level (figure 1; Null and Wurtsbaugh, 2020).The downward trend in water-surface elevation is expected to continue as human population and water consumption increase under changing climatic conditions in the state.Land subsidence and earth fissures due to long-term groundwater pumping in excess of recharge have been reported in Cedar Valley in southern Utah, where average basin-wide subsidence is estimated to continue at a rate of 0.04 to 2.4 inches per year under the current rates of groundwater decline (3 feet per year; Lund and others, 2011).Several instances of land subsidence have been reported in Woods Cross City (see figure 10 for location), but current subsidence rates are unknown.Several longterm monitoring sites on GSL's east shore show a significant increase in TDS over time (Clark and others, 1990), but saltwater intrusion into the freshwater aquifer has not been explicitly documented.However, groundwater-level declines on the east shore may create the conditions required to induce saltwater intrusion.Such phenomena have been observed in Lake Urmia, a terminal saline lake in Iran, where groundwater-level declines of 13 feet in the freshwater aquifer induced saltwater intrusion (Ahmadi and others, 2022).The general trends in groundwater levels shown in figure 6 and the 100-foot contour intervals shown in figures 4 and 5 allow visualization of groundwater conditions in the study area.Unfortunately, 100-foot contours do not provide good spatial resolution on fine-scale responses to groundwater withdrawals in the basin-fill aquifers.Similarly, mean groundwater levels on a decadal timescale provide a glimpse of the hydrological conditions at specific sites, but do not offer a detailed temporal resolution of groundwaterlevel responses to withdrawal or recharge.Groundwater-level maps at finer spatial and temporal scales than presented here are needed for each individual valley adjoining GSL to track withdrawal responses on a yearly (or seasonal) basis.Thorough and more comprehensive groundwater-level maps are particularly needed along the east shore of GSL to monitor potential saltwater intrusion to the freshwater aquifer.

Groundwater Flow and Salt Loading
The total groundwater flux to GSL calculated in this study (an average of 356,000 acre-feet per year from the seepage calculations, 313,500 acre-feet per year from the Darcy flux calculations) is substantially higher than previous estimates derived from water balance studies (75,000 acre-feet per year; Waddell and Fields, 1976;Loving and others, 2000).Arnow and Stephens (1990) estimated that between 6250 and 100,000 acre-feet of groundwater enter the lake per year.However, our results are consistent with recent estimates derived from geochemical modeling of GSL water chemistry (10% of the total inflow or ~300,000 acre-feet per year), assuming surface water contributions of 2.9 million acre-feet per year to the system (Stephens and Gardner, 2007;Bunce, 2022).Estimates of groundwater seepage using seepage meters averaged 0.77 cm/day from July 8-15, 2010, (Anderson, 2012) at locations of suspected groundwater seepage.Extrapolating that rate to the entire lake at its current coverage area would suggest that about 560,000 acre-feet of groundwater seeps into the lake per year.However, the seepage estimates could vary spatially, seasonally, and temporally, and were likely biased towards higher seepage rates due to how the measurement sites were selected.
Based on our preliminary estimates, the largest groundwater contribution originates in the north and east sections of GSL (table 5) where steep hydraulic gradients occur (figures 4 and 5).The groundwater flow derived in the north section of GSL was unexpectedly high (table 2).This high number is likely explained by the high transmissivity values calculated for the areas around the Park and Curlew Valleys (figure 7d).Groundwater flow conditions in Park Valley are poorly known.Extensive groundwater pumping in Curlew Valley has substantially reduced groundwater levels and discharge from the Locomotive Springs complex during the past 40 years (Hurlow and Burk, 2008), therefore the groundwater flux from this valley to GSL has likely declined significantly.
The average salt contribution from groundwater to GSL calculated in this study (1.18 million metric tons per year) represents about 10% of the solutes delivered by the Jordan, Bear, and Weber Rivers to GSL in 2013 (14.3 million metric tons;Shope and Angeroth, 2015).The highest TDS concentrations are found in the west/south sections of GSL (table 5) where hydraulic gradients are shallow, evaporation rates are high, and recharge likely occurs at a slow rate.
Our new estimates of groundwater discharge and its salinity contribution will likely require revision of GSL's water and salt budgets.However, two im-portant considerations limit how these should be evaluated and used.
1. Error on the estimates of groundwater flux is large.Based on sensitivity analyses from iterative calculations, the estimates are most sensitive to values of hydraulic conductivity.Hydraulic conductivity is lognormally distributed and can range by orders of magnitude over a study area as large as ours, resulting in estimates of flux that range over orders of magnitude.Of the iterative calculations conducted, 90% of resultant estimated flux fell between 2350 and 44,000,000 acre-ft/yr, indicating a need to better constrain aquifer properties.2. Because we only considered wells completed in the basin-fill aquifer, our calculations do not include flow paths that are entirely within bedrock (but do include groundwater that discharges from bedrock to basin fill in the subsurface).Significant discharge from bedrock springs occurs in the southeastern part of the Malad-Lower Bear River Valley, along the margins of the Promontory Mountains, and in the northwest part of Tooele Valley.These springs contribute groundwater flow to the GSL playa and, perhaps, different salt loading having different compositions and concentrations than groundwater in the basin-fill aquifer.Other aquifer properties, including porosity, saturated aquifer thickness, and cross-sectional area, also influence flux estimates.Further information is needed to constrain the differences in porosity values around GSL's shorelands and in the basin-fill aquifer.These data can potentially improve the seepage velocity estimates around GSL. Additionally, the aquifer thickness values we used (difference between depth to water and total well depth) have two sources of uncertainty.First, most of the wells in the basin fill only partially penetrate the saturated thickness of the aquifer.Thus, aquifer thickness values around GSL are likely larger than we estimated and could result in larger groundwater fluxes than presented here.Second, most of the wells used in this study were drilled to target the most productive aquifer depths.For some areas, this would result in larger groundwater and salinity contributions than expected due to bias toward higher aquifer-property values.Groundwater levels change over time due to natural and anthropogenic influences, resulting in variable hydraulic gradients and saturated thicknesses.Using previous work and ongoing groundwater-level observations, we can constrain these estimates fairly well.Cross-sectional area of groundwater flow paths is more complicated, especially if one assumes that the area matches that of the lake margin.The lake perimeter varies dramatically depending on lake level.Increases in lake level will increase the lake perimeter length and the resultant cross-sectional area of estimation.Hypothetically, an increase in lake level would result in a decreased groundwater gradient, but due to lack of lake-margin groundwater data, we are unsure of this relationship.
Regarding salinity inputs from groundwater, the spring systems around GSL need further consideration.Springs in the area have been measured to reach TDS concentrations of up to ~76,000 mg/L (Bunce, 2022) and are point-sources of solutes to GSL.It is also possible that there are density differences in the groundwater system around and below GSL.These density variations have the potential to create flow boundaries that we did not account for (Rosen, 1994;Sheibani and others, 2020).Further studies are needed to understand spring dynamics and density-driven flows in order to provide further insight on their overall role in the water and salinity budgets in GSL.

CONCLUSIONS AND RECOMMENDATIONS
We provide the first systematic, basin-wide assessment of groundwater levels in areas adjoining GSL to quantify groundwater contributions and their salt loading to the lake system.We observed groundwater-level declines in most of the basin-fill aquifer from the 1980s to the 2010s (figure 6).These declines are consistent with historical groundwater level trends in the Salt Lake, Tooele, Curlew and Weber Valleys and result from aquifer overdraft associated with less than average precipitation in the basin and increased withdrawal for human consumption.We calculated a mean groundwater flux to GSL of 356,000 acre-feet per year using a seepage velocity method in ArcGIS and 313,500 acre-feet per year using the Darcy flux equation for linear traverses through adjoining valleys.Both estimates are substantially greater than previous estimates derived from water balance studies, but are consistent with estimates derived from previous in situ seepage measurements and geochemical modeling of GSL water chemistry.We calculated a salt contribution from groundwater to GSL of 1.06 million metric tons per year which represents about 10% of the solutes derived from surface flows to GSL in 2013.These estimates have very large uncertainty, and the input parameters need to be better understood and constrained.Groundwater monitoring wells and a formal groundwater model are recommended to constrain groundwater parameters.
The data presented here have the potential to improve current water and salt budgets for GSL's system.However, further work is needed to improve these estimates and better delineate surface/ groundwater dynamics in the area.In order to do so, we recommend the following: 1. Estimates of the hydraulic properties of the basin-fill aquifer should be refined by compiling results from high-quality well tests and aquifer tests and, perhaps, conducting new aquifer tests.2. A monitoring-well network should be established with local/state/federal participation within each valley adjoining GSL.This system of wells should be thorough and accessible to visit and measure water level fluctuations on a seasonal or, at least, yearly basis.This information could be used to create detailed, yearly potentiometric surface maps and track/compare changes in groundwater levels over the years.3. Nested piezometers and/or monitoring wells should be installed along different sections of GSL.These piezometers at different depths could be used to calculate hydraulic gradients and to monitor water level/salinity trends in areas of the aquifer susceptible to brine intrusion.Core or cuttings recovered during the installation of these piezometers/monitoring wells could be used to estimate porosity values in the subsurface near GSL. 4. Sample springs around GSL and measure their flow.Geochemical and isotopic data on springs can provide information regarding sources, flow paths, residence time of groundwater and help to better understand their role in the water/salt budget of GSL.

Figure 1 .
Figure 1.Physiographic overview of GSL study area.

Figure 2 .
Figure 2. Geographical location of historical groundwater-level sites.

Figure 3 .
Figure 3. Number of groundwater-level sites by decade.

Figure 4 .
Figure 4. Potentiometric surface elevation in study area for the 1980s.

Figure 5 .
Figure 5. Potentiometric surface elevation in study area for the 2010s.

Figure 6 .
Figure 6.Groundwater level difference between the 1980s and 2010s.

Figure 7 .
Figure 7. Raster datasets used as input to calculate groundwater seepage velocity around GSL in ArcGIS Pro®: a) USGS water level sites, b) potentiometric surface, c) saturated thickness, and d) transmissivity.All data needed to create these raster datasets are available in tableS2.

Figure 9 .
Figure 9. Potentiometric surface elevation in the Salt Lake Valley area for the 1980s.

Figure 10 .
Figure 10.Potentiometric surface elevation in the Salt Lake Valley area for the 2010s.

Figure 11 .
Figure 11.Potentiometric surface elevation in the Weber Delta area for the 1980s.

Figure 12 .
Figure 12.Potentiometric surface elevation in the Weber Valley area for the 2010s.

Table 1 .
Groundwater flux estimates for GSL in acre-feet per year.
Figure 8. Polygons used in Darcy's Law equation to estimate groundwater discharge from each basin-fill valley adjoining GSL.Colors around GSL show divisions used to estimate groundwater and salt contributions by section.

Table 2 .
Groundwater flux estimates for GSL by section at 4190 FASL in acre-feet per year.

Table 3 .
Groundwater contribution from adjoining areas in acre-foot per year H.A. Zamora and P. Inkenbrandt Estimate of Groundwater Flow and Salinity Contribution to Great Salt Lake

Table 4 .
Salt flux estimates for GSL in metric tons per year.Average values are in boldface.

Table 5 .
Salt flux estimates for GSL by section at 4190 FASL in acre-feet per year