A Geographic Information Systems ( GIS )-Based Approach to Derivative Map Production and Visualizing Bedrock Topography within the Town of Rutland , Vermont , USA

Many state and national geological surveys produce map products from surficial and bedrock geologic maps as a value-added deliverable for a variety of stakeholders. Improvements in powerful geostatistical exploratory tools and robust three-dimensional capabilities within geographic information systems (GIS) can facilitate the production of derivative products. In addition to providing access to geostatistical functions, many software packages are also capable of rendering three-dimensional visualizations using spatially distributed point data. A GIS-based approach using ESRI’s ® Geostatistical Analyst ® was used to create derivative maps depicting surficial overburden, bedrock topography, and potentiometric surface using well data and bedrock exposures. This methodology describes the importance and relevance of creating three-dimensional visualizations in tandem with traditional two-dimensional map products. These 3D products are especially useful for town managers and planners—often unfamiliar with interpreting two-dimensional geologic map products—so they can better visualize and understand the relationships between surficial overburden and potential groundwater resources.


Introduction
Many state geological surveys produce derivative maps from surficial and bedrock geologic map sheets for value-added deliverables because they request matching funds from the State Geologic Survey Mapping (STATEMAP) Component.This program encourages the production of derivative products to ensure mapping activities support issues of social relevance [1,2].There is usually a gap between the time necessary for fieldwork and the subsequent creation of a geologic map as well as the availability of resources required to create various derivative maps.Derivative maps may include the spatial distribution of radioactive elements, slope stability, aggregate resources, surficial sediment characteristics, aquifer potential, or potentiometric surface.These derivative maps were historically created using hand-drawn contours and spatially distributed data collected at points; as outcrops, soil pits, or well locations.The creation of a continuous raster data layer covering the entire mapping area using these point data can be expedited and automated using geographic information system (GIS) based interpolation and/or extrapolation techniques.
The use of interpolation and extrapolation functions with point data offers a powerful opportunity for visualizing subsurface geologic relationships.These techniques have been used to characterize subsurface features for many decades; however, improvements in GIS software offer access to powerful geostatistical exploratory tools and robust 3D capabilities that can facilitate the production of three-dimensional visualizations using spatially distributed data [3].The GIS-based approach described in this project is particularly useful for creating derivative maps that depict surficial overburden, bedrock topography, and potentiometric surfaces [4][5][6][7][8][9][10][11][12][13][14][15][16][17].

Study Area
The town of Rutland is located in south-central Vermont within a north-south trending valley primarily underlain by Paleozoic rocks of the Shelburne, Danby, Winooski, Monkton, Dunham and Dalton formations.There are rare outcroppings of Ordovician Hortonville and Bascom formations, and Proterozoic gneiss of the Mt Holly Complex [18,19].The surficial geology is primarily dominated by glacial till; however, there are extensive deposits of lacustrine sediment, kame terraces, isolated exposures of kame deposits, alluvium and colluvium [20].
Elevations within the town range between approximately 476-1,430 feet (145-435 m) above mean sea level with the greatest relief along the northwest edge of the town limits (Figure 1).However, the town is primarily characterized by lower elevations and rolling hills with subdued topography.East Creek and its tributaries drain southward to Otter Creek and then northward to Lake Champlain.Otter Creek traverses approximately 10 miles (15 km) through town, flowing from the southeast to the northwest with a relief of approximately 50 feet (15 m).
Residential development within the town is tied to an expansion of bedroom communities, vacation homes, and commercial lodging that supports tourism associated with local ski resorts and derivative maps derived geologic maps will likely play an important role in helping manage town-related water resources.

Methodology
To better characterize the subsurface hydrologic characteristics within a GIS, 552 private, municipal and exploratory boring wells were spatially rectified and coordinates for 946 bedrock outcrops were collected using a Trimble ® JunoST running TerraSync™.The surficial geology of the town was mapped and 15 different surficial units were identified using traditional field and digital mapping techniques.Lastly, information about well depth, overburden, underlying bedrock geology and depth to water table was extracted from rectified well logs.Bedrock exposures were used with well log data to develop a representation of surficial overburden, bedrock topography, and local potentiometric surface.All raster analyses were performed using a 30-m digital elevation model (DEM) obtained from the Vermont Center for Geographic Information (VCGI).A raster layer representing the thickness of surficial sediments (also called an overburden or isopach map) was created using: (1) surficial overburden information from well logs; (2) bedrock exposures identified during field mapping; and (3) bedrock outcrops mapped by previous workers [18,21].To facilitate the process of estimating overburden and providing a raster layer that covers the entire map area and not just those areas with wells, the ESRI ® ArcGIS ® Geostatistical Analyst extension was used to determine which kriging function was best suited for these data.Ordinary kriging and universal kriging techniques produced similar error prediction values; however the final visual product varied considerably.For both techniques, the J-Bessel variogram model was used because it provided the best fit based on prediction error results.The ordinary kriging technique, while more commonly cited in the literature [5,10,17,22], produced more jagged contours regardless of the chosen variogram model or extent of smoothing.The universal kriging technique produced a smoother surface and resulted in less abrupt contours [23,24].Contour lines were manually smoothed using the smooth function, part of the ArcGIS ® Advanced Editing tools, to produce the final isopach map (Figure 2).Some authors argue that certain techniques may be more "visually faithful to reality" even though their statistical behavior is not always the best.Both models were statistically similar, so universal kriging was used because it produced the least visually jagged product [25].
Training and testing subsets for 11 variogram models were split 80/20 for cross validation and to evaluate validation prediction results.A variogram model that provides accurate predictions should report mean prediction error values close to zero if the predictions are unbiased, root-mean-square standardized prediction error values should be close to one if the standard errors are accurate, and root-mean-square prediction error values should be closer to one [3].If average standard error values are greater than root-mean square prediction error values, the model overestimates variance in the predicted values.If average standard error values are lower than root-mean square prediction error values, the model underestimates variance in the predicted values [3].Tables 1 and 2 summarize the results of cross-validation of the 11 variogram models used to decide which model to use in creating the isopach surface.A raster layer representing bedrock topography was created by subtracting the derived isopach layer from the DEM and 100-foot (30 m) contours were created from the resulting bedrock DEM using the contour function within the ArcGIS ® Spatial Analyst extension.
A potentiometric surface was created using static water levels recorded in well logs and a 30-meter DEM.Depth to the static water in each well was subtracted from DEM cell directly beneath the well location and joined to the attribute table to identify the elevation of water within each well.To facilitate the process of potentiometric surface production and provide a surface covering the entire map area and not just those areas with wells, this surface was created using inverse distance weighting (IDW) and the resulting raster layer was used to create 50-foot (15 m) contours following [4,9,13,16,22].Groundwater flow lines, which indicate potential flow along hydraulic gradient within an aquifer, were manually drawn perpendicular to potentiometric contours.

Isopach Map
The data is not normally distributed; it is strongly weighted towards thin overburden and bedrock exposures and there is also a trend indicating thicker overburden in the eastern part of town decreasing to the west (Figure 3).Prediction values exhibit most variance in the western half of the town because of limited well data and the fact that the region contains abundant bedrock outcrops and areas of thicker till.Bedrock topography generally mimics the landscape, especially in areas of thin overburden.However, it is often difficult to identify differences between the DEM representing surface topography and the DEM representing bedrock topography unless the data is visualized as a three-dimensional model (Figure 4).

Potentiometric Surface
The interpolated potentiometric surface is consistent with well data, surficial geology, and bedrock topography.Typically a potentiometric surface does not characterize the actual surface of the water table but is a proxy for the potential hydraulic head within an aquifer [10,26].The hydraulic gradient is gentle and the contours are widely spaced and therefore uncertainty exists in the inferred flow direction in some areas of the map.However, the general trend of flowing east to west and northeast to southwest is readily apparent (Figure 5).A potentiometric surface gradient depicted using filled 50-foot (15 m) contours and resulting flowlines, which were created using the ArcMap ® Geostatistical Analyst and manually drawn respectively.

Discussion
The interpolated isopach map is consistent with both well data and the spatial extent of known surficial units.The areas of thickest overburden in the map occur along the eastern edge of the field area and in the southwestern corner of town.These areas are covered by thick kame terrace deposits in the east and thick glacial till and glaciofluvial deposits in the southwest.Prediction values exhibit most variance in the western half of the town because of limited well data and the presence of both abundant bedrock outcrops and areas of thicker till.Well data is also limited in the central portion of the study area because that section utilizes municipal water and not privately owned wells.Areas of thin overburden in the central and northwestern corner of town are also validated through field evidence as a thin veneer of till riddled with abundant bedrock outcrops.There are a number of areas lacking well data where abundant bedrock outcrops substantially improved the final distribution of overburden and ultimately bedrock topography.Similar to many regions in Vermont and New England, bedrock topography within the town generally mimics surface topography with rare exceptions.Both universal and ordinary kriging methods quickly produce a raster representation of overburden and an acceptable contour map when compared with well data and field relationships.It is possible that contours derived using either method could be further smoothed to create an isopach map, however in this study universal kriging resulted in the most aesthetically pleasing contours.Similarly, IDW produces a raster surface that can be used to quickly characterize hydraulic potential throughout a given field area.The advantage of using these techniques is that they provide an expedited alternative to traditional hand-drawn techniques and can be automated using ESRI's ® Model Builder or through Python scripting [13].
This automated approach is not without assumptions and limitations; these include the resolution of the underlying DEM, distribution of well log data, and date of well installation.I used a 30-m DEM in this study, which was the highest possible resolution for this area of Vermont.However, 10-m DEMs are available in many other states and access to DEMs with higher resolution produced using Light Detection and Ranging (LiDAR) technology is increasing in urban areas throughout the industrialized world.Higher resolution DEMs will better characterize subtle changes in topography and water levels within closely-spaced well log data, which are often homogenized using 10 and 30-m resolution data.Geostatistical approaches based on well data rely on the distribution of private and municipal wells, which are typically concentrated near population centers and not well-distributed throughout a given area of interest.The age of well log data is also problematic because state and town databases often include wells covering time periods spanning 40-60 years.Using static water level heights from wells drilled decades apart and drilled during different seasons to characterize the hydraulic gradient in areas with high groundwater withdrawals rates could produce erroneous renderings.Therefore a systematic analysis of derivative maps produced using a GIS versus those hand-drawn by geologists should be undertaken to better characterize the usefulness of these products.This would entail finding appropriate paper maps and georeferencing them to quantify the spatial accuracy of interpolated products using an error matrix analysis.
However, the final products produced in this study are consistent with field evidence and detailed inspection of well log data.The production gap between time spent in the field and publication of derivative maps will likely drive more geological surveys towards automated production workflows.The task of explaining two-dimensional geologic map products is often complicated by variations in the spatial visualization skills given the diversity in regional planners and town managers [27][28][29].However, many authors believe it is possible to strengthen and develop facility with spatial visualization, rather than thinking it is a fixed cognitive skill, through increased use of three-dimensional representations [11,14,[29][30][31].Although the methodology in this paper uses ESRI's ArcScene ® software to render the three-dimensional perspective in Figure 4, numerous other software packages offer similar capability, including open-source options such as Quantum GIS+GRASS using the NVIZ module and gvSIG.Therefore, many of these maps are best viewed in a three-dimensional environment so that planners and town managers unfamiliar with interpreting two-dimensional geologic map products can better visualize the relationships between surficial overburden and subsurface hydrologic characteristics [32,33].

Conclusion
The availability of powerful geostatistical exploratory and three-dimensional tools within modern GIS software facilitates the production of derivative products.In addition to traditional derivative maps depicting overburden, bedrock topography, and potentiometric surfaces, the production of three-dimensional visualizations is especially useful to help non-experts better visualize and understand the relationships between surficial and subsurface features.Although the methodology described in this article specifically addresses groundwater resources it could be applied to visualizing subsurface faults, geochemical variations or ore-body orientation using public or private well data.Creating and distributing derivative map products that non-experts are able to utilize may encourage additional national or state level investment as the community of stakeholders grows.

Figure 1 .
Figure 1.A map depicting the study area in south-central Vermont, which is located in the northeastern United States.

Figure 2 .
Figure 2.An illustration of using the isopach layer created using universal kriging to produce contours for the northern half of Rutland Town (A) before smoothing and the results of (B) manually smoothing the data using the built in smooth function of ArcGIS ® .

Figure 3 .
Figure 3. (A)A histogram of overburden thickness from well data and bedrock outcrop locations in north Rutland created using the ArcMap ® Geostatistical Analyst.This illustrates a non-normal distribution influenced by thin till cover and over-sampling of bedrock outcrops to increase control on overburden.(B) A trend analysis performed using the Geostatistical Analyst Extension indicates a trend that overburden is thicker in the east and decreases to the west.

Figure 4 .
Figure 4.A three-dimensional representation of bedrock topography relative to surface topography created using ArcScene ® software.The green represents the overburden layer and brown represents bedrock topography.The area most strongly influenced by extensive bedrock outcrops is visible in the center of the image A and two areas of thickest overburden are visible in the south and east-central areas of town B and the northeast corner of town C. The red points represent the location of wells distributed throughout the town of Rutland used to estimate overburden.The vertical gap between the two layers represents the thickness in overburden.Note: vertical exaggeration is 1.25.

Figure 5 .
Figure 5.A potentiometric surface gradient depicted using filled 50-foot (15 m) contours and resulting flowlines, which were created using the ArcMap ® Geostatistical Analyst and manually drawn respectively.

Table 1 .
Summary of prediction error values reported using the cross-validation function of the ArcGIS ® Geostatistical Analyst for each variogram type using the training dataset.Note: RMS = root mean square error.

Table 2 .
Summary of prediction error values reported using the validation function of the ArcGIS ® Geostatistical Analyst for each variogram type using the testing dataset.Note: RMS = root mean square error.