Inverse modeling using PS-InSAR for improved calibration of hydraulic parameters and prediction of future subsidence for Las Vegas Valley , USA

Las Vegas Valley has had a long history of surface deformation due to groundwater pumping that began in the early 20th century. After nearly 80 years of pumping, PS-InSAR interferograms have revealed detailed and complex spatial patterns of subsidence in the Las Vegas Valley area that do not coincide with major pumping regions. High spatial and temporal resolution subsidence observations from InSAR and hydraulic head data were used to inversely calibrate transmissivities (T ), elastic and inelastic skeletal storage coefficients (Ske and Skv) of the developed-zone aquifer and conductance (CR) of the basin-fill faults for the entire Las Vegas basin. The results indicate that the subsidence observations from PS-InSAR are extremely beneficial for accurately quantifying hydraulic parameters, and the model calibration results are far more accurate than when using only water-levels as observations, and just a few random subsidence observations. Future predictions of land subsidence to year 2030 were made on the basis of existing pumping patterns and rates. Simulation results suggests that subsidence will continue in northwest subsidence bowl area, which is expected to undergo an additional 11.3 cm of subsidence. Even mitigation measures that include artificial recharge and reduced pumping do not significantly reduce the compaction in the northwest subsidence bowl. This is due to the slow draining of thick confining units in the region. However, a small amount of uplift of 0.4 cm is expected in the North and Central bowl areas over the next 20 years.


Introduction
Las Vegas Valley (Fig. 1), encompassing an area of 4150 km 2 in southern Nevada, has experienced a long history of rapid population growth and consequent growing rates of groundwater pumping that have led to large water-level declines of as much as 90 m (Burbey, 1995) and subsequent extensive and complex subsidence patterns with four regional subsidence bowls, with the Northwest bowl experiencing more than 1.5 m of land subsidence (Bell et al., 2002) in this structural basin filled with more than a thousand meters of unconsolidated to semi-consolidated heterogeneous sediments (Fig. 1b) (Amelung et al., 1999).In recent decades, mitigation strategies including aquifer storage and recovery (ASR) (Brothers and Katzer, 1990) practices have led to widespread recovery of water levels and substantial but uneven rebound of the land surface in some localities, while other localities continue to experience subsidence (Hoffmann et al., 2001).PS-InSAR reveals that from 2002 to 2010 much of the valley has experienced uplift or at least stabilization of the land surface due to water-level recoveries, but the Northwest and Southern bowls continue to subside (Fig. 2) even though the major pumping centers do not coincide with these subsiding areas.Thus the complex relation between pumping locations and rates, water-level changes, and surface deformation patterns are a clear reflection of the highly heterogeneous hydrogeologic units underlying the valley, which include basin fill faults and their poromechanical properties.
Past numerical models of Las Vegas Valley (Jeng, 1998;Morgan and Dettinger, 1996) have not been able to adequately capture these complex subsidence patterns and therefore their effectiveness as groundwater management tools Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.In this investigation we use newly acquired basin-wide temporal PS-InSAR and water-level data from 2002 to 2010 as observations to inversely estimate aquifer parameters from a more finely discretized groundwater flow model, which extends the previous modeling efforts of Yan (2008) and Yan and Burbey (2008).The large number of model cells and the computational requirements of the coupled forward and inverse models necessitated the use of the high-speed com- puting facilities at Virginia Tech.The goal is to create a robust groundwater management model for the entire Las Vegas basin that accurately reflects the subsidence and recovery patterns observed during the past 100 years of aquifer development.

Numerical model for Las Vegas Valley
For this study we modify and extend the existing Las Vegas model developed by Jeng (Jeng, 1998) and later by Yan (Yan, 2008;Yan and Burbey, 2008) to account for data through 2010.The model simulates groundwater flow using MOD-FLOW_2005 (Harbaugh, 2005) and land subsidence using the SUB package (Hoffmann et al., 2003) for MODFLOW and considers the influence of basin-fill faults that are simulated using the HFB package (Hsieh and Freckleton, 1993) for the entire Las Vegas Valley using a 600 m × 600 m grid cell size.A four-layer model is used to represent the hydrogeologic conditions of the valley where the uppermost two layers represent the near-surface aquifer, layer three is the developed (principal) zone aquifer that extends 200-300 m below land surface and which is the primary focus of the calibration effort in this investigation because this is the layer from which all groundwater is pumped and virtually all subsidence occurs.Layer four is the deep-zone aquifer and extends to the basement rocks beneath the basin-fill sediments.The lateral boundaries are all prescribed as no-flow conditions because the basin fill is truncated by the surrounding carbonate and volcanic mountains that are considered to be of much lower permeability than the basin-fill sediments.A total of 87 stress periods is used to simulate the hydrogeologic conditions in the basin that includes nearly 100 years of aquifer development, from March 1912 to October 2010.

Initial conditions and observations
The key parameters that affect both water-level changes and land subsidence distributions throughout the valley and which are calibrated using inverse modeling are the transmissivities (T ), elastic (S ke ) and inelastic (S kv ) skeletal storage coefficients and fault conductances of the developed-zone (or principal) aquifer (layer 3), which are believed to be largely responsible for the observed deformations and water level distributions.A total of 72 T parameters (zones) and six S ke and S kv parameters (zones) from the Yan model (Yan, 2008) were used as initial parameter estimates for this investigation.The zones and parameter values for all the other layers were kept the same as the Yan model.Fault conductance values were initially set to be equivalent to the adjacent aquifer conductance (or transmissivity) values.The large discrepancy between the number of transmissivity and storage zones was based on the limited clay thickness and compressibility data necessary for estimating storage values.
Water levels from both domestic and municipal wells were provided by the Las Vegas Valley Water District and the U.S. Geological Survey covering more than 50 years are used as observations for both the forward and inverse calibration procedures.New land surface deformations acquired from In-SAR and PS-InSAR data processing were obtained from 58 ENVISAT satellite acquisitions (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) in descending track mode and were processed using ROI_PAC software (Rosen et al., 2004) to supplement the existing InSAR data from previous investigators (Bell et al., 2002).These new data (Fig. 2) were processed to obtain deformation phase data for each permanent scatterer (PS) to detect the average velocity fields from the time-series data.The PS-InSAR data were then filtered to remove topographic errors, atmospheric errors, the phase noise introduced by the filtering operation, to correct for elevation-dependent atmospheric effects and the orbit tilts.The final deformations were land-truthed with GPS data over a seven-year span.A 40 m resolution is used in both the azimuth and range directions to describe the subsidence distribution for Las Vegas Valley.These 40 m pixels were then averaged over each model grid cell to provide a widespread set of temporal and spatial subsidence (or uplift) values used as observations for the inverse model.Detailed monthly InSAR data are used to investigate how the seasonal variations in water levels are reflected in subsidence and rebound patterns.These elastic stress-strain signals are used to estimate the S ke , which depends more heavily on the seasonal response of the developed zone aquifer that includes clay interbeds.proc-iahs.net/372/411/2015/Proc.IAHS, 372, 411-416, 2015 Iterative parameter zone reconstruction, which consists of either coarsening or refining the existing parameter zone number and distributions, was accomplished during the calibration process using an adaptive multi-scale algorithm (Ameur et al., 2002).Parameter inversion was accomplished using UCODE_2005 (Poeter et al., 2005).Prior information was used for parameters T and S se and S sv at some localities obtained from Morgan and Dettinger (1996).All simulations were conducted at Virginia Tech's ARC High Performance Computing Center due to both the size and scale of the iterative problem being solved.Transmissivity zone distributions from Yan's model were refined through 15 iterations of successive approximations during the inverse procedure to calibrate the final T values for layer 3. Computed CSS (composite scaled sensitivity) values are all larger than 0.021 times the largest CSS, indicating that the model and the observed data provide sufficient information to estimate the transmissivities.Similarly, the S ke and S kv of layer 3 were calibrated at the basin scale using the multi-scale algorithm, which was accomplished through five iterations of successively refining and/or coarsening the zone domains.CSS values were all found to be larger than 0.11 times the largest CSS with S kv and suggests that the inelastic skeletal storage coefficient is highly sensitive to the observations and perhaps the most important parameter for achieving accurate calibration of the model.
Fault conductance values near the Eglington Fault in NW Las Vegas Valley (Fig. 1) were calibrated using eight zones and assuming a 100 m wide fault zone.The calibrated transmissivity values ranged over two orders of magnitude (from 1 to 100 m 2 day −1 ) in this region, indicating that the fault zones can possess different hydraulic characteristics than the adjacent aquifers.

Results and predictions
The final calibrated transmissivity and storage zones and values are provided in Fig. 3. Final simulated water levels are compared with observed values across the entire Las Vegas Valley in Fig. 4 after the years 1982 and 2006.The Nash-Sutcliffe efficiency index (NS) (Nash and Sutcliffe, 1970), mass balance error (m) and the Schulz criterion (D) (Schulz et al., 1999) are used to evaluate model fit.These criteria suggest that the model fit is heads is excellent.The final simulated land subsidence distributions are shown in Fig. 5, which shows a similar pattern as developed by Bell (2002) based on observed values, except in the eastern part of the basin where larger simulated subsidence values occur.The NS and m criteria indicates that the fit between the simulated and observed land subsidence is very good.
Basin fill faults are shown to have an influence on the hydrogeologic setting of the basin.All four main subsidence bowls (Northwest, North Las Vegas, Central and Southern     Bell (2002).Bell (2002), who implied that fault gouge or secondary carbonate cementation of the fault zone may exist, or perhaps some other mineralization may be occurring along the fault to cause it to act as a subsidence barrier.However, no significant hydraulic head discontinuity is found to exist across the fault indicating that the influence of the fault as a hydraulic barrier may not be significant or some more complex process such as strain partitioning may be contributing to the large subsidence gradients adjacent to several subsidence bowls.
The North bowl located not far from the pumping/recharge center has recently undergone land uplift and hydraulic head recovery (Fig. 2).At the Central and Southern bowls, groundwater levels have generally recovered, so the occurrence of land subsidence in the southern bowl suggests the existence of residual compaction as a result of slowly draining clay interbeds.The lack of significant subsidence/uplift at or near the center of active pumping/recharge could be explained by the absence of thick clay interbeds.The interbed thickness is large in the eastern part of the valley.The lack of significant subsidence in this area can be explained by the general lack of appreciable water level decline and by the fact that it juxtaposes a bedrock boundary.
Future trends in both water level and land subsidence were investigated by extending the current pumping and artificial recharge rates and patterns from 2011 to 2030.Secondary recharge rates, which have grown substantially since 1972, are based on a single average projected growth rate for 2011-2030 by extrapolating the expected rate from the 1972-2010 rates.Figure 5a shows the map of the predicted groundwater levels changes for the period 2011-2030, which generally shows an overall recovery of water levels of as much as 50 m.Even the water levels in the Northwest subsidence bowl show an increase in water levels.For example, location h1 (Fig. 6a) undergoes a head recovery of 14 m during this 20-year period.Figure 6b shows the predicted subsidence changes for this same 20-year period, which suggests that residual compaction from slowly draining confining units and clay interbeds will continue into the future leading to continued predicted subsidence of 11.3 cm in the Northwest subsidence bowl in spite of the water-level recovery in this area.However, uplift of approximately 0.4 cm will occur in the North and Central bowls.At h1 (Fig. 6a) the subsidence is estimated to be 6.3 cm in spite of the large recovery.

Conclusions
An updated numerical groundwater and subsidence model was created by updating observations of water level and land subsidence through the processing of new InSAR scenes for 2002-2010 and refined discretization.Early InSAR data revealed that subsidence was widespread across the region into the mid-90's and two subsidence bowls were shown to be controlled by basin fill faults and perhaps differences in interbed thickness.Recovery of water levels due to decreased pumping and increased artificial recharge are responsible for the observed uplift in recent years.However, residual compaction is occurring in areas where thick clay units are still expelling their water, such as the in the Northwest subsidence bowl (Fig. 2).Calibrated aquifer parameters of T , S ke and S kv are similar to Yan's model except near basin-fill faults where skeletal storage coefficients are generally lower than those of Yan.This may be due to abrupt changes in clay thickness or transmissivity differences along the fault zone.
Overall, we find that high resolution temporal and spatial subsidence observations from PS-InSAR are extremely beneficial for quantifying the complex patterns and values proc-iahs.net/372/411/2015/Proc.IAHS, 372, 411-416, 2015 of aquifer properties in the Las Vegas basin, which provide considerably more information than water levels alone and therefore results in a more robustly calibrated model.Water levels are expected to recover into the future as are land surface elevations except in the Northwest subsidence bowl area where the draining of thick clay layers will continue to result in compaction in spite of notable water level recovery in the area.

Figure 1 .
Figure 1.Site map showing location of basin fill faults (top) and A-A cross section (bottom) through Las Vegas Valley.

Figure 2 .
Figure 2. PS-InSAR interferogram showing rate of surface uplift and subsidence for the period 2002-2010.
bowls) are bounded by basin-fill faults and are significantly offset from the center of artificial recharge and center of municipal pumping.For example, on the upthrown block of the Eglington fault, calibrated S ke andS kv values are higher than those fromYan (2008), so it is inferred that thicker interbeds (clay) may exist at this locality and may represent the primary factor controlling the formation of the Northwest subsidence bowl.Calibrated transmissivity values for the fault barriers indicate that a conductance barrier may exist along the upthrown block of the Eglington fault as suggested by

Figure 6 .
Figure 6.(a) is the predicted water level from 2011 to 2030 and (b) represents the predicted subsidence or uplift for the period 2011-2030.