Aquifer Storage Calculations Using GIS and Modflow, Los Angeles County, California

Theodore A. Johnson and Wanjiru M. Njuguna

GIS and Modflow were used to determine the additional storage capacity of the Central and West Coast basin aquifers in coastal Los Angeles County, California, as part of conjunctive use planning. The Water Replenishment District of Southern California wants to store excess surface water in the underlying aquifers for future extraction during drought. Modflow layers containing aquifer elevations, storage coefficients, water levels, and a land surface Digital Elevation Model (DEM) were used in conjunction with ArcView Spatial Analyst to create model cells and perform grid and volume calculations. Results show that more than 600,000 acre-feet of additional groundwater can be stored in the aquifers.


1.0 INTRODUCTION

Groundwater management agencies are responsible for ensuring that a sufficient supply of high quality groundwater exits in the underlying aquifers to meet the current and future needs of the overlying users. With increasing populations, decreasing water availability and the susceptibility of surface water systems such as aqueducts, reservoirs, and pipelines to natural disasters or terrorist acts, the protection and optimization of a healthy groundwater basin is crucial. Old and new conjunctive use technologies are being employed to ensure the effective management of surface water/groundwater systems, including aquifer storage and recover (ASR), computer modeling linked to optimization algorithms, artificial recharge in spreading grounds, injection wells, or in-lieu pumping, recycled water use, and water banking/storage programs.

The Water Replenishment District of Southern California (WRD) is the groundwater management agency for the Central Groundwater Basin (CB) and West Coast Groundwater Basin (WCB) of coastal Los Angeles County, California (Figure 1). Of the 730,000 acre-feet per year (afy) of total water used by the overlying 4 million people in 43 cities, about one-third is obtained from over 400 groundwater pumping wells that are located throughout the area. The remaining two-thirds of the water supply is imported from Northern California and the Colorado River via aqueducts and pipelines since the natural and augmented groundwater supply can not meet the full annual demands.

Project Location Map

Because the cost of the imported water continues to rise, the availability of Colorado River water is reduced as Nevada and Arizona take more of their entitled share, and to ensure adequate supplies in times of drought, WRD must find ways to optimize use of the groundwater in the CB and WCB. One of the ways to optimize groundwater in the basins is to fill the aquifers with less expensive, excess surface water when it is available so that it can pumped out months or years later when it is needed. This conjunctive use strategy of temporarily storing groundwater is known as water banking.

To determine the amount of groundwater that can be stored in the basin aquifers requires a detailed understanding of the three-dimensional hydrogeologic regime of the basins. WRD has determined the amount of additional groundwater that can be stored in the CB and WCB by studying the basins’ hydrogeology, building and calibrating a Modflow computer model to simulate groundwater movement, and use of ArcView Spatial Analyst to perform volume calculations using the model files. The remainder of this paper describes the methodology used by WRD to determine the available storage space in the CB and WCB.

2.0 GROUNDWATER STORAGE CONCEPTS

The concepts of groundwater storage used in the analysis are important to understand. Although the final determination is basically a volume calculation of an aquifer’s available storage space (length x width x height), hydrogeologic aspects must be brought into the equations for a thorough analysis. These include the top and bottom elevations of the aquifers or model layers, the land surface elevation (LSE), the groundwater levels in each of the aquifers, whether the aquifer is confined or unconfined, and values for the specific yield and confined storage coefficients.

Figure 2 displays the concepts used in the analysis. As it would be undesirable to store excessive groundwater in the basins so that water levels rose to the land surface causing potential flooding, WRD maximized the rise of groundwater levels to a depth no higher than 50 feet from the land surface (LSE-50). The amount of available storage, therefore, was calculated from the current groundwater levels up to LSE-50.

Aquifer Types

In an unconfined aquifer such as portions of the CB and WCB, large amounts of groundwater can be stored in the open pore spaces between the clay, silt, sand and gravel grains. The amount that can be stored and later extracted is based on the effective porosity or specific yield (SY) values of the aquifers, which vary depending on grain size, sorting, roundness of the grains, and interconnectedness of the pore throats. Typical values of specific yield in CB and WCB type sediments of clays, silts, sands, and gravel, range from 0.03 to 0.38 (Todd, 1980).

In a confined aquifer, such as the majority of the CB and WCB, much less water can be stored because the aquifer is completely full. The only availability to force more water into the aquifer is dependent on the storage coefficient (SC) of the confined aquifer, which describes the compressibility of the mineral skeleton of the aquifer matrix and the expansion of the water. Values of SC are much less than SY, typically ranging from about 0.005 to 0.00005 for CB and WCB type aquifers (Todd, 1980).

As shown on Figure 2, to determine the amount of additional groundwater that can be stored in the CB and WCB aquifers, the basin had to be evaluated for areas that were fully unconfined (Area A), mixed confined and unconfined (Area B), fully confined (Area C), and where the basin was already too full (current groundwater elevations were higher than LSE-50, Area D). Information on the aquifer top and bottom elevations, the model-calculated groundwater elevations, and the SY and SC values that were used for the available storage calculations were obtained from the groundwater flow model described in the next section.

3.0 GROUNDWATER FLOW MODEL

The WRD, in cooperation with the United States Geological Survey (USGS), has developed a groundwater flow model of the coastal plain of Los Angeles County, which includes the CB, WCB, Santa Monica Basin, and the Hollywood Basin (Reichard, et al, 2002, in press). The model was a six year effort of extensive data collection and analysis managed in an ArcInfo GIS incorporated into the USGS Modflow code, which is a commonly used groundwater modeling program in the groundwater industry. The model was constructed to better understand the three-dimensional hydrogeologic system of the Los Angeles coastal plain, to simulate current and future groundwater occurrence and movement, to determine the groundwater balance, and to provide WRD with a tool to evaluate numerous groundwater management scenarios related to artificial recharge, seawater intrusion control optimization, and conjunctive use strategies.

The model domain is shown on Figure 3. For the model construction, the domain was discretized into ½ mile by ½ mile grids in 67 rows, 70 columns, and 4 layers to produce 18,760 individual cells. As shown on Figure 4, the vertical layers were constructed to simulate the four main aquifer systems, including the Recent Holocene (Gaspur aquifer), Upper Pleistocene Lakewood Formation (including the major Gage and Gardena aquifers), the Lower Pleistocene Upper San Pedro Formation (including the major Lynwood and Silverado aquifers), and the Lower Pleistocene Lower San Pedro Formation (including the Sunnyside aquifer). Because the aquitards were not implicitly modeled between the layers, the model is considered a quasi-3D model. Groundwater movement between layers is controlled by the vertical conductance term in Modflow.

Model Grid

Basin Cross Section

Physical properties for each aquifer or layer in each grid cell must be entered into the model. The following data grids were required: Top of layer elevation, bottom of layer elevation, hydraulic conductivity, SY, and SC. Much of this information came from California Department of Water Resources (1961) with updates by Reichard, et. al (2002, in press). Boundary conditions were established as either no-flow, head-dependent, or general heads. For the Pacific Ocean, freshwater equivalent general heads based on the bathymetry were used. Structural barriers to groundwater flow such as faults (Newport-Inglewood, Overland, Charnock, and Palos Verdes) and folds were added using Modflow’s horizontal flow barrier package. The recharge components of spreading ground percolation, seawater barrier injection wells, aerial recharge, and mountain front recharge were added. Production wells were the only discharge component added. If a well was perforated across multiple aquifer layers, the percent of total water pumped from a particular layer was determined based on the transmissivity ratios of the perforated aquifers.

After construction of the model, it was run for a steady state calibration for water year (October 1 through September 30) 1970/71, which was a relatively stable year prior to the initiation of the Dominguez Gap barrier well system. Physical parameters of hydraulic conductivity and aerial recharge were adjusted within set bound limits until a suitable calibration was reached when compared to actual water levels for that year. Transient calibrations were then run with annual stress periods from 1970/71 through 2000/2001, adjusting SY and SC values primarily to obtain good matches between model generated water levels and hydrographs of 50 actual wells located throughout the model domain with sufficient historical record to compare against the model results.

The calibrated model has been used for numerous groundwater management studies over the past several years and has been instrumental in providing a better understanding of the water balance and three-dimensional groundwater flow systems. It has been used to determine the optimum groundwater pumping patterns near the coast to minimize the adverse effect of seawater intrusion. It has been used to study the impacts of increasing groundwater production in the CB by 20%, and has been used to determine the amount of underflow occurring across other groundwater basins. The next significant use of the model will be to evaluate the potential effects of storing additional groundwater in certain areas of the basin based on the results of the storage calculations described in the next section.

4.0 DETERMINATION OF AVAILABLE GROUNDWATER STORAGE

Grids created in Modflow (described above) were imported into ArcView 3.2 with the intention of calculating basin storage volumes for the WCB and CB. In ArcView, the extension Spatial Analyst version 1.1 was used to create additional grids and perform algebraic functions using these grids as variables. The desired result of these calculations is to be a visual description of the various layers in the basin (i.e the represented aquifer surfaces) and a numerical estimate for the amount of acre feet available in the basin in which to store water.

4.1 Grids

The specific grids generated in Modflow used for the available storage calculation are described below in the following paragraphs. The first grid is the groundwater elevation grid (GWE - Figure 5). This grid depicts the model calculated groundwater elevation in each of the 4 model layers in the LA basin. The cell size for this grid is ½ mile squared (804.672 meters squared) and each cell contains an elevation value for groundwater. In situations where the model cell was dry during the simulation, a value of -888.88 was assigned to the cell.

Groundwater Elevation Grid

The next grid generated in Modflow is the top of layer grid (TOL - Figure 6). This grid depicts the upper surface of the model layer using elevation data to generate the grid. TOL grids exist for layers 2, 3, and 4, as Layer 1 is fully unconfined and has no top elevation except for the groundwater surface. The cell size for this grid is also 804.672 meters squared and is of the same extent as the GWE grid. Related to this layer is the Bottom of layer grid (BOL - Figure 7). Similarly, this surface is a grid that depicts the base elevation of the model layer. This grid is composed of cells sized 804.672 meters squared.

Top of Layer Grid

Bottom of Layer Grid

The next two grid surfaces generated in Modflow for the storage volume calculations are the specific yield (SY) and storage coefficient (SC) grids (Figure 8). The specific yield surface contains cells that hold numerical values ranging from 0.075 to 0.25. These values represent a ratio that describes the volume of water capable of draining from a saturated area as a result of the force of gravity on that area. This ratio, multiplied by the area of unconfined aquifers results in the amount of acre feet of available storage space in that particular region of the basin. Similar to the specific yield grid is the grid for the confined storage coefficient (SC) values throughout the basin. The values in the grid range from 0.00015 to 0.002952.

Storage Coefficient Grids

In contrast to the grid layers describing subsurface features, the layer that was used to describe the land surface elevation is the 24k Digital Elevation Model obtained from the United States Geological Survey. For the purposes of this project, we generated a new elevation grid 50 feet below the ground surface known as Land Surface Elevation minus 50 (LSE-50). For this project, it was decided that LSE-50 would be the uppermost level at which groundwater could rise without having shallow water problems such as liquefaction potential, flowing artesian wells, or risk damage to substructures. So the space in the basins available for additional groundwater storage is the space between the groundwater elevation in each model layer and the LSE-50 surface.

4.2 Masks

The next step in this process was to import the grids into ArcView 3.2. This was accomplished using a Modflow pre and post processor called GMS. First, 2D grids were created in GMS using the 3D modflow data. These grids were exported as an ascii file (*.asc). In Spatial Analyst, importing the file as an ascii raster data source opened these grids. Masks were then created on the various grids to isolate the various scenarios that are present in the basins. These scenarios include areas where the basin is unconfined, confined and mixed (both confined and unconfined), areas where storage is available in the basin and areas where cells are dry versus wet (non-dry).

Prior to generating the masks, the cell size of the Modflow grids was converted from 804.672 meters squared to 100.584 meters squared. This was done to standardize the Modflow grids with the ArcView grids and to provide a sufficiently small cell size for accurate calculations. We found through trial and error that by using the coarser grid sizes, errors were obtained using map calculator when one grid (groundwater elevation for example) was a different size or origin than another grid. Therefore, we compromised with a smaller grid size with consistent origins for all grids without getting too fine with the grid spacing to avoid unnecessary file sizes upon grid calculations.

4.3 Analysis Process

In ArcView, the cell sizes for the surfaces generated in Modflow (GWE, TOL, BOL, SY and SC) were all converted to a cell size of 100.584 meters squared. This was done by setting the analysis properties to "100.584" and using Map Calculator to "evaluate" the particular surface. This process generates a new surface with the desired cell size. With this standardizing process accomplished we began the analysis for the unconfined, confined and mixed situations. A summary of the process is described below.

4.3.1 Unconfined Aquifer

The first aquifer scenario that we generated was the surface depicting the available storage for the unconfined portions of the basin (Area A on Figure 2) where the Modflow cells were wet (non-dry). To begin, we generated a mask that isolated the dry cells versus the non-dry cells. This was done by assigning a value of "1" to the cells which were non-dry, and a No Data value to the cells that were dry. Once we had isolated the unconfined, non-dry cells the next procedure was to locate the areas within these cells that could potentially hold more water. Calculating the areas in which the GWE was less than LSE-50 provided this result. The result of this calculation was a surface that showed unconfined, non-dry (wet) cells that had storage space available. The amount of acre feet of additional groundwater storage space that this represented was calculated using the equation below.

Unconfined non-dry area storage = (((LSE-50)-GWE)*SY)*108,897.352/43,560

The value 108,897.352 is a constant that converts the one dimensional cell value into a volume (cubic feet) by multiplying the value by the area of the cell (100.584 meters squared) and then converts from square meters to square feet. Finally, the result, which is in cubic feet, is converted to acre feet by dividing by 43,560.

The next task was to calculate possible storage for the dry cells (areas where Modflow determined that the GWE dropped below BOL, so the cell was dry). We did this by reclassifying Dry cells versus Non-dry cells in the GWE grid to assigning the dry cells a value of "1" and the non-dry cells a value of "No Data". This mask was then designated to be the Analysis mask. To calculate the areas of the dry cells where storage was available, we used Map Calculator to calculate the areas where the Bottom of Layer (BOL) is greater that the Land Surface Elevation minus 50 (LSE-50). Situations where this is true indicate that the BOL is too high for any possible storage. False values indicate areas where the BOL is lower than the LSE-50 grid and storage is available.

To calculate the actual amount of storage available, we reclassified the BOL>(LSE-50) and assigned the False values (where storage is available) with a value of "1" and then created a mask that isolated these areas. Using the mask of this area, the available storage was determined using map calculator and the following equation:

Unconfined dry area storage = (((LSE-50)-BOL)*SY)*108,897.352/43,560

4.3.2 Confined Aquifer

The second aquifer condition for which we calculated available storage is the confined areas (Area C in Figure 2). The initial step was to determine the areas where storage is available. These areas are the areas where the GWE lies below the (LSE-50) grid. Reclassifying the data and assigning these areas with a value of "1" and the areas with GWE above (LSE-50) with a value of No Data isolated the areas for which this was true. Once areas of available storage had been isolated, the totally confined portions of this area need to be identified and isolated. To find the totally confined areas we calculated areas where the GWE lies above the TOL. Once again we reclassified the values and isolated the areas that had groundwater levels above the TOL. Storage in acre feet for the totally confined areas was calculated using the following equation:L

Confined area storage = (((LSE-50)-GWE)*SC)*108,897.352/43,560

Because the SC grid has very small values, the resulting calculations for acre-feet of storage in each cell were also small, sometimes less than 0.5. We found that during the process of having to convert the floating point grids to integer grids before final processing, much of the cell data were lost (values less than 0.5 were set to 0 as an integer) and the results were under reported. To remedy this situation, we multiplied the resulting floating point grid by a value of 1,000 before converting to an integer grid, and then divided the results by 1,000 as a final step.

4.3.3 Mixed Aquifer Situation

The final scenario for which we calculated aquifer storage was the mixed confined and unconfined aquifer situations (Area B of Figure 2). Areas of unconfined aquifer storage were masked off, and volumes calculated as described in Section 4.3.1. Then, the areas of confined aquifer storage were identified and masked off, and volumes calculated as described in Section 4.3.2. The results of the two steps were combined to produce the total available storage in these mixed areas.

5.0 RESULTS

The procedures described in Section 4 were repeated for a total of four separate layers representing aquifer systems that lie within the Central and West Coast basins. A value of available storage was calculated in each grid cell for each model layer (Figure 9).

Basin Storage by Layer

It was then desirable to summarize the available storage by sub basin regions. For the amount of storage available for the sub basin regions, a mask was created of each of the sub basins and used in conjunction with the total storage grid to arrive at the acre feet value of storage for each sub basin. The Central Basin was found to have 464,840 acre feet of available storage, and the West Coast Basin was found to have 155,329 acre feet. Figure 10 and Table 1 present the results of this analysis.

Total Available Storage

Table 1 - Results

6.0 CONCLUSIONS

Effective groundwater management involves optimizing conjunctive use programs to store excess surface water when available for later use. Storage of surface water in the underground aquifers is a proven conjunctive use management operation. Before this can occur, the amount of storage space available for use in the aquifers must be determined.

This paper presented the methodology used to determine the amount of available storage space in the Central and West Coast groundwater basins of coastal Los Angeles County, California. One conclusion drawn from the analysis is that there are over 600,000 acre feet of potential storage availability in the basins.

However, it is recognized that this amount would be difficult, if not impossible, to fully implement due to groundwater levels not being able to rise equally everywhere to the maximum height of 50 feet from the ground surface. Instead, water levels rise sharpest at their recharge points (injection wells or spreading grounds), and it may take considerable time for water to rise substantially away from these points. Therefore, the 600,000 acre foot value is considered the maximum value. Additional work can now be done through computer modeling and other studies to refine and optimize specific conjunctive use projects to maximize groundwater storage in the basins.

7.0 ACKNOWLEDGEMENTS

Several people contributed significantly to the work done for this project and their efforts are sincerely appreciated, including Tracy Parr (formerly with WRD) and Melody Kneale (Komex H2O Science) for their GIS applications, Eric Reichard and Steve Crawford of the U.S.G.S. for their Modflow work, and Jason Weeks of WRD for his technical review of this paper.

8.0 REFERENCES

Todd, D. K., 1980, Groundwater Hydrology, second edition. John Wiley & Sons

California Department of Water Resources, Bulletin No. 104, 1961. Planned Utilization of the Ground Water Basins of the Coastal Plain of Los Angeles County, Appendix A - Ground Water Geology.

Reichard, E., Land, M., Crawford, S., Schipke-Paybins, K., Nishikawa, T., Everett, R., Johnson, T., 2002, Geohydrology, Geochemistry, and Ground-Water Simulation-optimization of the Central and West Coast Basins, Los Angeles County, California; U.S.Geological Suvey Water Resources Investigations Report 02-xxx (unpublished to date).


Theodore A. Johnson
Senior Hydrogeologist
Water Replenishment District of Southern California

Wanjiru M. Njuguna
Assistant GIS Specialist
Water Replenishment District of Southern California