Digital Surfaces, Thicknesses, and Extents of Geologic Model Units


Model output includes model-generated surficial geology over 53,030 mi² of the Columbia Plateau and the extent of the CRBG within this area (fig. 8), elevation of top of CRBG units and the Older Bedrock unit (figs. 9A–D), and thickness of CRBG units and the Overburden unit (figs. 10A–D). Histograms (figs. 11A–D) and an example of the spatial distribution (fig. 12) of mismatch between the surfaces and data points (residuals) provide a measure of error between the model-generated surfaces and the supporting data, which can be used to estimate the probability of encountering each surface within a prescribed interval in a well. Cross sections through the model domain (fig. 13) were generated. The Mabton and Vantage Interbed unit data were analyzed, but known discontinuity and uncertainty made mapping of these units explicitly infeasible. Instead, general thickness relations were developed for these units.


All model-generated surfaces and thicknesses are stored in GIS grid format, with each square grid cell being 500 ft on a side. The model-generated surficial geology shows the intersection of each of the model units with the land surface. The intersection of the extents of all CRBG units forms the model-generated extent of the CRBG (fig. 8), which shows good agreement with the mapped surficial geology (fig. 1A). The computed area of the model-generated CRBG polygon is 42,064 mi², which agrees well with the estimate of 44,000 mi² for the CPRAS by Kahle and others (2009).


Computed areas and extents for the model units generally compare favorably to the estimates of Kahle and others (2009). The Grande Ronde Basalt unit with a computed area of 41,866 mi² compares favorably to the estimated area of 42,000 mi². It underlies almost all of the 42,064 mi² computed CRBG extent. The only areas not underlain by the Grande Ronde Basalt unit are at the model periphery, where intracanyon flows of younger CRBG rock overlie the Older Bedrock unit directly. The computed Wanapum area of 24,379 mi² also is close to the estimated area of 25,000 mi². The Saddle Mountain unit has a computed area of 11,668 mi², contrasting with the estimate of 8,000 mi². Contrary to the estimated extent of Kahle and others (2009), the geologic model predicts that many of the topographic highs on the Palouse Slope may be Saddle Mountains Basalt, which is supported to a limited extent by identification of this unit in drillers’ logs and on geologic maps. Most of this area is geologically mapped as being mantled with sediment, so either interpretation is considered reasonable. Whether the rock is Wanapum or Saddle Mountains Basalt, a regional groundwater flow model likely will be insensitive to this distinction because the unit will be simulated as basalt either way. For the Overburden unit, the shapes of thick sedimentary sequences that are model-generated and estimated by Kahle and others (2009) are sufficiently different that the change in shape likely is more important for controlling groundwater flow and exchange with surface water than for estimating total area. For each unit, the model volumes are Grande Ronde (31,273 mi³), Wanapum (2,156 mi³), Saddle Mountains (663 mi³), and Overburden (244 mi³).


Model-Generated Elevations of Tops of Units


Each model-generated top of CRBG and Older Bedrock (figs. 9A–D) is exposed at land surface over part of the model area and buried over other parts. When the unit is buried, the model-generated top elevation is smooth, reflecting uncertainty in the model-generated trend surface. Contours show estimated elevations where buried. Where the unit is exposed at land surface, the elevation of the raster is the land-surface elevation resampled to the 500-foot grid. The Overburden unit is always at land surface (fig. 8), so no separate map of the top of this unit is provided.


Model-Generated Thicknesses of Units


Each model-generated thickness (figs. 10A–D) was computed as the difference between the unit top and the top of the uppermost underlying units. Whenever the top is at land surface, the top inherits the properties of the resampled DEM, resulting in thickness maps with significant rugosity. For this reason, only shaded maps (no contours) are provided. The top of the Older Bedrock unit is the lower bound of the model, so no thickness map is available for this unit.


Evaluation of Fit


The well and geologic data were compared to the final surfaces by computing the difference between the unit top elevation from the data and the elevation predicted by the model. This difference, called a residual, is positive if the data are higher than the model surface and negative if the data are lower than the model surface. The model surface is a good estimator of the data if the residuals are random and there is no spatial bias. These properties imply that the median, mean, and mode of the residuals are all approximately zero, the residuals are symmetrically distributed around the mean, and there are no persistent spatial trends in residuals. The histograms for the Columbia River Basalt units (figs. 11A–C) indicate that the trend model is a good estimator for these units, but there is some bias for the Older Bedrock trend surface (fig. 11D). Examination of the spatial distribution for all CRBG units indicates no strong spatial trends.


For example, the distribution for the top of Grand Ronde residuals (fig. 12) shows the model matches the data spatially. The top of Grande Ronde model elevation ranges from approximately -2,200 to 7,800 ft, with a total difference of about 10,000 ft. Most of the residuals are less than 100 ft (less than 1 percent of the trend surface being represented) and have no strong spatial trends. The sparse data in the center of the map are supplemented with guide points to aid in constructing reasonable geologic geometries (see section, “Using Thickness Maps for Quality Assurance and Model Revision”), but the guide points are not used when computing residuals to assess model fit.


The top of Grande Ronde Basalt unit model fit was provided as the example because it shows how all of the different types of data were used together during the modeling process. The fits for the Saddle Mountains and Wanapum Basalt units are very similar to the fit for the Grande Ronde Basalt unit, but the Older Bedrock unit top has larger error and spatial bias (fig. 11D). The worse fit of the Older Bedrock unit is the result of increased uncertainty caused by little data being available over most of the model area, resulting in heavy reliance on inferred data points. Reidel and others (2002) provide a generalized map of total CRBG thickness that was used to infer the top of Older Bedrock using the model-generated CRBG unit tops. Use of this data increased the spread of the residuals, because well data did not always agree with the generalized trend.


The heavy tail of positive residuals on the histogram (fig. 11D) for the Older Bedrock unit is the result of most of the well picks for that unit (about 250 wells) being located near the northern and northeastern extent of the Columbia Plateau where deposition was dominated by intracanyon lava flows rather than by the sheet flows present over much of the study area. The Older Bedrock unit in this area forms topographic highs around which the CRBG lavas flowed and, even where it is buried, is expected to have highly variable elevation. Attempts were made to get a better fit to data in this area, but the result was that the CRBG units began pinching out anomalously. Data density is insufficient to define the true shape of the Older Bedrock top in this area.


Rather than forcing a better fit and causing geologically unreasonable pinchouts, the bias toward positive residuals was permitted in this area for two reasons. First, positive residuals imply that the modeled surface is lower than the data, but because wells commonly are drilled no deeper than necessary, more frequent encounters with the bedrock highs cause a bias of the data toward higher elevations. Sample bias toward the higher values indicates that the trend is lower than the data imply. Second, allowing the bias toward positive residuals facilitates the use of modeled geology as the foundation for a groundwater flow model. Allowing the biased residuals at the model periphery result prevents anomalous pinchouts of the overlying model unit, resulting in continuity of the model unit, allowing the possibility that this unit may transmit water during groundwater-flow simulation modeling. If sufficient groundwater-level data are available during the groundwater-flow simulation process, the degree of conductivity in this area will resolve itself. 


The total elevation range of about 25,000 ft for the top of the Older Bedrock unit implies that the larger 500 ft of error corresponds to about 2 percent, which still provides a good general representation of this unit. However, uncertainty in this trend surface elevation is more substantial than for the tops of the CRBG units because of the low density of data (less than 300 wells, with most representing only a small part of the study area).


In the absence of strong spatial trends in the residuals, the implicit assumptions underlying the kriging paradigm are satisfied. If local estimates of uncertainty are required, then kriging may be used as the exact interpolator for the residuals, so that the kriging variance would be an estimator of local uncertainty, allowing computation of error bars on interpreted surfaces. Furthermore, Sequential Gaussian Simulation may be used to construct stochastic realizations for use in other uncertainty assessment exercises as desired (Deutsch, 2002). More simply, because the mean value of the residuals is approximately zero, and assuming that the residuals are stationary (supported by absence of spatial trends), the trend surface and the distribution of residuals may be used to infer the likelihood that a contact is within any given range of elevations as supported by available local data. The estimate of contact elevation is more reliable if there are nearby wells where the unit has been identified. In locations where little data are available, the trend surface is the best estimator in an unbiased sense, but uncertainty increases with distance from data.


Geologic Cross Sections


The ten representative geologic cross sections in figure 13 show general trends and variability of geologic units across the CPRAS. The modeling assumption that all faults are vertical is best illustrated by cross sections passing through the Yakima Fold Belt structural region (fig. 8 and cross sections A-A’ and B-B’ in fig. 13). 


Mabton Interbed Unit


Generally, the Mabton interbed is everywhere that the overlying Saddle Mountains Basalt occurs, although locally, it may not be present because of local variability in the paleodepositional environment. In the western fault-bounded basins, thickness patterns are visually apparent and are correlated to thicker deposits of Saddle Mountains Basalt. Variability is relatively high, with anomalously thick deposits being reported in some wells with no corroborating evidence in nearby wells.


On a regional scale, it is assumed that interbeds transmit water much less efficiently than the CRBG aquifers and commonly are classified as “confining units.” Because interbeds are sedimentary units, however, interbeds may store an appreciable amount of water per unit volume that may be released during transient conditions. Therefore, if the volume of the unit is significant at a given location, it may supply an appreciable amount of water to the aquifer system. The volume of the Mabton interbed is the property needed for groundwater flow simulation, and the discontinuous nature of the deposits is less important.


To better understand the typical volume of the interbed, the thickness of the interbed must be estimated. It is assumed that thicker interbeds correspond to thicker deposits of the overlying basalt. A best linear fit between all Mabton interbed thickness data and modeled Saddle Mountains Basalt thickness (constrained to pass through the origin) was constructed (fig. 14) to estimate the relationship. Next, the 10 percent of the data with the worst fit was removed to minimize the outlier effect, followed by fitting an unconstrained line to the best 90 percent of the data. The equation of the best line is


,     (1)


where units are in feet. The R-squared value for this equation is 0.34. Constraining the line to pass through the origin yields a best-fit slope of 10 percent (meaning that the thickness of the Mabton is expected to be 10 percent of the thickness of the overlying Saddle Mountain Basalt) (fig. 14), with an R-squared value of 0.25. Although the correlation between thick Mabton interbeds and thick Saddle Mountains Basalt is evident when considering the spatial distribution of Mabton interbed deposits across the study area, the poor R-squared values result from the high variability of the interbed thicknesses at shorter spatial scales than is well defined by the trend models of the buried units. This variability in sedimentary interbed thickness likely was controlled by local hills and valleys in the paleotopography.


Because all data used for this analysis have reported Mabton interbed thicknesses, wells where the Mabton interbed does not occur are not represented, resulting in bias towards overestimating the average Mabton interbed thickness. This bias implies that the relations above are upper bounds for the average behavior of the Mabton interbed thickness. The assumption that the Mabton interbed is 10 percent (larger of the two computed line slopes in fig. 14) of the total distance between the top of the Wanapum Basalt and top of Saddle Mountains Basalt (fig. 10B) yields a maximum thickness of less than 300 ft, which is still lower than the largest reported values.


Vantage Interbed Unit


Similar to the Mabton interbed, the Vantage interbed generally is everywhere that the overlying Wanapum Basalt occurs, although locally it may not be present because of local variability in the paleodepositional environment. Again, the typical thickness of the interbed is the property needed for the hydrogeologic framework of the groundwater flow model, so a similar approach to that used for the Mabton Interbed unit was followed. In contrast to the Mabton interbed, the Vantage interbed shows no strong patterns in map-view or strong correlations to the thickness of Wanapum Basalt unit (fig. 15). Variability is high, with anomalously thick deposits being reported in wells with little or no corroborating evidence in nearby wells. The same procedure used for the Mabton Interbed unit to remove data outliers was used for the Vantage Interbed unit, with the equation of the best line being

,     (2)

where units are in feet. The R-squared value for equation (2) is 0.004, supporting the observation that there is little or no correlation between the thickness of the Vantage interbed and the overlying Wanapum Basalt thickness. Because there are no apparent regional spatial patterns of Vantage interbed thickness, the best estimators of thickness are assumed to be the mean value of 19.5 ft and the median value of 15 ft, which are well supported by the near-constant value of the best-fit line. Again, zero-thickness values were not used in the analysis, indicating that this estimate may be an overprediction on average. Again, the poor R-squared value is an indication of the inability of the geologic model to represent paleotopographic features that controlled the deposition of thick and thin sedimentary deposits.


Applicability and Limitations


The geologic model integrates a large number of different types of data to create a fully consistent three-dimensional representation of the large-scale geometry of the CPRAS. As such, it is the synthesis of the best current understanding gained from previous studies at various scales. Its error is quantifiable (fig. 11A-D), with the error being a function of the amount of data. The following list of limitations should be considered when using the results of this study:


• The data used to construct the surfaces contain errors. This fact is known because some studies did not agree on lithologic picks. These errors may cause bias in the interpreted surfaces. Many of the lithologic picks made for particular units are inferential on the basis of previous work, which may be the source of some errors.

• In most cases, data are not closely enough spaced to resolve features at a smaller scale than the trend model.

• Generally, error increases with depth because of limited amounts of data in many areas. Increased error with depth is evidenced in the residual histograms (fig. 11A–D). The top of the Older Bedrock unit has some bias and may be underestimated in the north and northwest where CRBG depositional style has changed from flood basalt to intracanyon flows.

• Error increases with distance from data. The error implied by the residual histograms may not be representative at points far from supporting data because the implicit assumption is that the trend surface is correct on average.

• Error in heavily folded or faulted areas is greater than for smoothly varying areas with moderate data support. The dataset was insufficient to resolve these features.

• The surfaces and thickness are generated using smooth trend surfaces representing unit tops and consequently result in smooth tops where buried. The model-generated values represent expected values at any subsurface location, but the true value likely will be somewhat different. Confidence that a unit will be located within a given range of elevations may be computed from the residuals.


The resulting geologic model provides important insights into the following:


• The model preserves important groundwater modeling characteristics.—If the trend model residuals are random and stationary, then, on average, the geometry of the system is correct. Good estimates of volume allow an improved understanding of aquifer storage parameters. Preservation of connectivity of large-scale features allows an improved understanding of potential groundwater flow paths within the system and connection with surface water.

• The error and uncertainty of the model-generated units can be estimated.—The model was created by creating trend surfaces. If the residuals are random and stationary, then the error of the model-generated surfaces is assumed to have the same statistical properties as the residuals, which forms the basis for quantifying the error. 


First posted February 25, 2011