Analysis of geothermal energy resources in the Appalachian Basin of the eastern United States is of interest, given the region’s population- and climate-driven demand for thermal energy. This study provides a fuller picture of geothermal resources across New York and Pennsylvania than previous studies by providing a rigorous statistical analysis of temperature-depth data using records from nearly 8000 locations. The compilation of thousands of temperature-depth data enables a significant increase in the spatial resolution of geothermal resource assessment maps for this region. In addition, this project has contributed to the compilation of geothermal data at a national level through the National Geothermal Data System. These temperature-depth measurements are byproducts of historical and recent drilling for petroleum and natural gas in the sedimentary basin. Bottom hole temperatures (BHTs) were recorded before the wells reached thermal equilibrium and at a wide range of depths. To extract a comprehensive description of the thermal state of the Appalachian Basin strata required application of both a BHT correction scheme and a simple thermal model. The model results for individual wells were combined with geostatistical interpolation employing kriging to produce maps that reveal significant variations in subsurface thermal gradient and surface heat flow with markedly improved spatial resolution. An area in south-central New York State displays favorable geothermal resource potential, with heat flow estimates of 50–60 mW/m2. There are 2 elongate, 200–300 km long, northeast-trending bands of favorable geothermal resource potential in central and western Pennsylvania, with heat flow of 55–90 mW/m2.


A demonstration that current fossil fuel–based energy sources are potentially replaceable with relatively low grade geothermal resources in the eastern United States could have a large impact on the future energy development in the region. The most attractive application area for these low-grade resources would be to provide thermal energy for district heating and hot water supply to residential and commercial buildings, as well as for other low-temperature industrial processes. Models that combine resource approximations with infrastructure and utilization economics suggest that a combination of thermal direct use and combined heat and power applications may generate economically viable opportunities for the use of geothermal energy near some population centers in the northeastern United States (Fox et al., 2011; Reber et al., 2014). Whereas electricity generation in current competitive markets typically requires geothermal wellhead temperatures in excess of 150 °C, direct use of geothermal heat for space heating, water heating, and industrial and agricultural processes can be achieved economically with wellhead temperatures as low as 80 °C. In general, targeting low-temperature thermal resources reduces required drill depths, and subsequently risk and cost. The objective of this paper is to identify areas in the Appalachian Basin of Pennsylvania and New York where geothermal gradient and heat flow are most favorable for further exploration.

One strategy for geothermal exploration in the eastern United States is to focus on locations where there is low-thermal-conductivity rock in the shallow subsurface, thereby causing a local increase in the geothermal gradient. Sedimentary basins such as the Appalachian Basin have interbedded mudstone and therefore meet this condition. Furthermore, some parts of the Appalachian Basin have large population centers, which provide a potential market for direct-use geothermal heat (Fig. 1). If sufficiently high temperatures are found at relatively shallow depths in a sedimentary basin, circulation of fluids needed to extract the geothermal heat could be achieved either through natural higher permeability (e.g., hot sedimentary aquifer resources) or through enhanced geothermal systems (EGS) technologies. The plausible existence of adequate heat delivery technologies and the numerous potential users motivate this reevaluation of the geothermal energy resources within the Appalachian Basin. This paper builds on studies presented in Stutz (2012) and Stutz et al. (2012).

This study updates the 2004 Geothermal Map of North America (Blackwell and Richards, 2004) for Pennsylvania and New York; the 2004 map incorporated 21 temperature-depth measurements in Pennsylvania and 29 measurements in New York, a very sparse data set for a land area that totals nearly 240,000 km2.

The northern Appalachian Basin study area is composed of Paleozoic marine and terrestrial strata overlying Grenville basement (Colton, 1970; McLelland et al., 2010). The sedimentary thickness ranges from 0 to 10 km, increasing steadily to the southeast and reaching its maximum thicknesses along the western edge of the Appalachian Mountains (Milici, 1996; Ryder et al., 2012). Regionally, a basal interval of upper Cambrian sandstones is overlain by Cambrian and Ordovician shallow-marine dolostone and limestone. The onset of the Taconic orogeny is related to the initiation of a thick stack of marine siliciclastic deposits. The marine siliciclastic deposits began with the moderately deep water lowermost Upper Ordovician Utica Shale and continued coarsening and shallowing upward until nonmarine conditions were achieved in the latest Ordovician Queenston Shale. Silurian units reflect a variety of shallow-marine and nonmarine environments. The most significant lithological departure from the early Paleozoic siliciclastic and carbonate units is found in the halite and calcium sulfates of the upper Silurian Salina Group. The widespread Devonian and lower Mississippian strata found in Pennsylvania reflect subsidence and fill of a new phase of foreland basin tectonic activity, known as the Acadian orogeny. While the abundant Devonian siliciclastic strata and sparse limestones reflect the northwestward progradation of a delta system across the basin, there were repeated transgressions of organic-rich mudstone facies eastward and southeastward; this likely caused the interbedding of low-thermal-conductivity rocks in the basin fill. The final phase of Appalachian Mountains tectonic shortening, the Alleghany orogeny, caused a new phase of foreland basin subsidence and fill during the Pennsylvanian; the major products were nonmarine siliciclastic units (Sweezey, 2002).

The Appalachian Basin strata are progressively more strongly folded and more influenced by thrust faults toward the southeast. The late Paleozoic Allegheny orogeny led to intense thrust faulting and folding in the Valley and Ridge Province, which deformed the eastern margin of the sedimentary basin (Shumaker, 1996). There was very little oil and gas exploration success in the Valley and Ridge Province, with the result that this study uses very few temperature-depth data within the highly deformed Valley and Ridge, and these few occur only near its western margin.

The highest quality geothermal gradient data are derived from borehole temperature readings collected in wells that are in thermal equilibrium (Roy et al., 1972). These equilibrium well temperature data are not available with enough data density coverage, so petroleum industry archival bottom hole temperature (BHT) data from well log reports were used to expand the temperature-depth data set for depths ranging between 600 and 4000 m. As of early 2011, temperature readings from ∼8000 wells deeper than 600 m had been extracted from largely predigital well log reports archived by New York and Pennsylvania state agencies. This 160-fold increase in observation locations represents a significant contribution to the analysis of geothermal resources in the area and should illuminate spatial trends and variability that were not resolved in the Geothermal Map of North America (Blackwell and Richards, 2004). In addition to the analysis presented here, the newly digitized data were contributed in 2012 to the National Geothermal Data System (NGDS; geothermaldata.org/), which serves as the central location for publicly available geothermal data.

The drilling technology for these deep wells sometimes involved a drilling fluid other than air to maintain hydrostatic pressure, lubricate the drill bit, and remove cuttings. When non-air drilling fluids are used, then near-field temperatures become significantly disturbed by circulation of these fluids. To record thermal equilibrium conditions, an amount of time significantly longer (e.g., several times multiplied) than the time span during which a non-air fluid circulated at the well bottom ought to have elapsed between end of circulation and the BHT measurement. However, only in a minority of the well records was it clarified which fluid was used, and for most of the Appalachian Basin wells there is incomplete information in the available records regarding the time elapsed between circulation and measurement. Therefore, the conservative assumption is that thermal equilibrium conditions were not routinely attained. Thus the BHT data for the study area are of uncertain, but likely variable, quality. An additional limitation of the reliance on well data from the petroleum industry is its omission of information about the basement rocks underlying the strata. These well data do not constrain the amount of radiogenic heat production in the underlying basement rocks that contributes to the spatial variability of heat flow.

The Google-sponsored national assessment of geothermal energy potential, completed by the Southern Methodist University Geothermal Laboratory in 2011, resulted in development of a set of methods that allowed regional assessment of heat flow and of depths to a temperature of interest based on consistent, rapid analysis of thousands of wells (Blackwell et al., 2007, 2011). The data set challenges are illustrated conceptually in Figure 2. Those methods relied on (1) simplified corrections to equilibrium conditions for measured well temperatures, (2) condensing lithologic variability through use of identified areas of uniform stratigraphy (Childs, 1985), (3) generalizing basement-heat input parameters, and (4) conversion of the depth-dependent temperature information to surface heat flow. Whereas this simplified approach facilitates a data-rich regional overview, the numerous details of geological heterogeneity that are overlooked when using these methods contribute to uncertainty of the results at all well locations.

Due to the nature of hydrocarbon resources, the BHT data points used in this study are not spatially homogeneous, but cluster in areas where oil and gas drilling were more prevalent, such as west-central Pennsylvania and western New York State, and at certain depths targeted for specific hydrocarbon plays. Therefore, the spatial distribution of the temperature-depth data is not ideal for sampling across the geologic features in New York and Pennsylvania. Nonetheless, variations across the two states were mapped using ordinary kriging (e.g., Gaussian process regression) (Appendix 1) interpolation methodologies with available BHT data. Efforts to use co-kriging employing information on rock type did not produce useful results.

Specific to the Appalachian Basin geology and data, this paper evaluates two classes of uncertainties. (1) An approximate correction of BHT to equilibrium thermal conditions that utilizes the Harrison et al. (1983) correction developed from another sedimentary basin is compared to data from available equilibrium profiles within the New York and Pennsylvania region (Frone and Blackwell, 2010). (2) The uncertainty related to interpolation between data points is investigated using appropriate statistical methods.


Estimation of the Thermal Field at Single Locations

Archived oil and gas well logs were assembled from collections held by the Geothermal Laboratory of Southern Methodist University, the Pennsylvania Geological Survey (Pennsylvania Internet Record Imaging System), the New York State Museum (Empire State Oil and Gas Information System), and the New York State Department of Environmental Conservation (NYSDEC, 2011). Initially, the BHT data were available only in scanned paper records, from which temperature data were read and recorded in spreadsheets. Only wells with BHT recorded at depths >600 m were compiled, to minimize the effects of groundwater movement and near-surface temperature variations on thermal gradient calculations (Frone and Blackwell, 2010). The resulting data set contained 8919 data points from 7969 wells. For the 745 wells for which temperatures were reported at multiple depths, only the deepest reported temperature was retained for analysis.

Estimation of the in situ temperature profiles using oil and gas BHT data is challenging (see Appendix 1). We use a depth-dependent correction developed in the Anadarko Basin (Harrison et al., 1983). The application of an Anadarko Basin empirical correction to another basin (e.g., Appalachian Basin) assumes that the BHT perturbations caused by fluid circulation were similar in the two basins at similar depths. A comparison of equilibrated temperature logs to BHT data for several wells in West Virginia (Blackwell et al., 2010) and 14 wells in New York and Pennsylvania (Fig. 3A) (Shope, 2012; Shope et al., 2012) demonstrates a reasonably good fit using the Anadarko Basin correction. However, the use of this empirical correction for our entire data set introduces a second set of assumptions, namely that most wells in the Appalachian Basin were drilled in roughly similar fashions and will have comparable magnitudes of deviation from thermal equilibrium. This second assumption is at odds with the fact that we do not know the drilling fluids or times since circulation for most wells in our data set.

The Harrison et al. (1983) correction (Equation 1) expresses ΔT (°C), the difference between equilibrium BHT and the observed BHT on a geophysical log header at depth zw (m), as a second-order polynomial: 

Although the deepest thermally equilibrated well in the northern Appalachian Basin, located in West Virginia, has a depth of ∼3600 m (Spicer, 1964; Blackwell et al., 2010), we extrapolated the correction to greater depths for the 13 wells in our data set deeper than 4000 m (maximum total depth of 6482 m) (Fig. 3B). Although it was noted in Shope et al. (2012) that the Harrison et al. (1983) correction may result in some underestimation of BHT for depths <1000 m (Fig.3A), the results presented here apply the correction uniformly across all BHT depths.

We calculate an average thermal gradient (dT/dz) at the location of each data point based on the approximation in Equation 2. The TBHT (°C) is the Harrison et al. (1983) corrected BHT, TS (°C) is the average annual surface temperature, and zw is the depth of the well: 
Details are given in Appendix 1.
A surface heat flow, Qs (mW m–2), at a given borehole can be calculated as the product of the thermal gradient and an average thermal conductivity value, forumla (W m–1 °C–1) (Equation 3): 
This necessitates a method (Appendix 1; Stutz et al., 2012) to estimate the average thermal conductivity, which depends on lithology, temperature, and pressure. The aggregate thermal resistance, ΣjRj (W–1 m2 °C) (see Appendix 1), for all the stacked lithologies in the borehole (or to an intermediate depth of interest) is used to calculate the average site thermal conductivity forumla from the ground surface to the BHT depth (Equation 4): 

Instead of interpreting the lithologies based on the available logs for each of the nearly 8000 wells that met the depth criteria, the appropriate lithologic succession was estimated from the nearest type column from the American Association of Petroleum Geologists Northern Appalachian Correlation of Stratigraphic Units of North America (COSUNA) section (see AAPG, 1985) (conceptually illustrated in Fig. 2). The appropriate COSUNA column was rescaled to fit the local depth-to-basement for each location, which was extracted from the map of the basement of North America (AAPG Basement Rock Project Committee, 1967). This approach incorporates regional-scale changes in relative thicknesses of the constituent formations but omits local-scale architectural changes within the formations or groups. The COSUNA lithology-thickness data above the depth of interest and lithology-specific thermal conductivities (Carter et al., 1998; Beardsmore and Cull, 2001) were key inputs to Equation 4.

Spatial Statistical Treatment

Prior to performing a spatial assessment of the geothermal resource potential within the region, local and global outlier analyses were conducted to detect potentially rogue observations in the data set (Appendix 1). As a result, 153 of the 7969 wells were removed as outliers. An additional 37 wells were removed from the data set because they were located in counties with fewer than 5 wells along the data-sparse borders of the selected study region. A total of 7779 wells remained in the data set for kriging interpolation.

Kriging was used to interpolate between well locations the values of various geothermal variables of interest, and to quantify the standard errors, a measure of the precision of the predicted values. In this study, kriging assumes a constant but unknown mean for the geothermal variable in the extent of the estimation window (often called ordinary kriging; details are described in Appendix 1). Variations from the mean are modeled with the use of a semivariogram that describes the behavior of geothermal variables as a stationary stochastic field (Goovaerts, 1997) (Appendix 1). In this study a maximum of n = 25 points was used to estimate the value at each grid location. A minimum of 20 points within a distance of 53 km were required to estimate the value. With these criteria, kriging was appropriate for interpolation with our geothermal products, which are limited in spatial extent but highly clustered in some areas and dispersed in others (Aguirre, 2014). The Gaussian semivariogram parameters for the nugget, sill, range, and anisotropy (see Appendix 1) for each geothermal variable of interest are summarized in Table 1. One expects the difference between observations that are closer to one another to be smaller than the difference between observations that are more distant from each other. For a set of data that fulfills that expectation, the nugget, range, and sill are terms that refer to the parameters of the semivariance curve (semivariogram) with increasing distance for all point-pairs in the dataset. If the data were themselves without error, if the thermal field remained constant through time, and if the quantity under study (e.g., thermal gradient) was perfectly consistent with all of the heat-flow assumptions implied by our simplified treatment, then in theory at infinitesimal distances between two wells the values calculated from the observed thermal field would be identical. The nugget, which is not the theoretical zero (Fig. 4), expresses the observed difference between the calculated thermal values over an infinitesimal distance. The range measures the distance within which there is spatial correlation between data points and beyond which there is no modeled spatial correlation, and the sill is the modeled value of semivariance at distances greater than the range.

Kriging was implemented in the R Foundation for Statistical Computing package gstat along with sp, spdep, lattice, and maptools (R Development Core Team, 2011). ArcMap 10.0 (ESRI Inc., http://www.esri.com/) was used to generate the maps of geothermal variables over the region.


Figure 4 shows the fit to the geothermal gradient data using a Gaussian semivariogram for binned point pairs (black circles) separated by increasing distances. The fitted Gaussian semivariogram indicates that all points within 10 km are modeled as having roughly the same semivariance, equal to the value of the inherent errors in the data set, such as measurement error and error in rock conductivity (e.g., Table 1, nugget; see Appendix 1). A degree of spatial dependence appropriate to the use of kriging as a spatial interpolator is illustrated by the good fit of the Gaussian semivariogram to the data to a distance of ∼55 km and by the increase in semivariance for point pairs separated by increasing distances toward a sill value. At large distances between point pairs the Gaussian fit is no longer adequate, as the point pairs reach a peak (the sill) and gradually decrease in value in a manner that does not match the Gaussian fit. This decreasing semivariogram behavior is indicative of no spatial dependence, or that other factors not considered in this analysis affect the structure of the semivariance (e.g., spatially varied geological features). Therefore, as there is no benefit to be gained by a broader estimation window, 53 km was adopted as the estimation window. The geothermal variables of interest in the Appalachian Basin display directionality (θ, ϕ). A stronger spatial dependence occurs in the southwest-northeast direction than in the southeast-northwest direction for the region as a whole (Aguirre, 2014).

Kriging the partial and final products of the single-well thermal model for nearly 8000 boreholes generated maps of the geothermal gradient (Fig. 5A) and heat flow (Fig. 6). Figure 5B illustrates the standard error of the predicted mean that resulted from the kriging interpolation for the geothermal gradient. For heat flow, the uncertainties associated with kriging interpolation (Fig. 6B) have a pattern like that in Figure 5B because the data density governs the precision of the interpolation (Aguirre, 2014).

The thermal gradient and heat flow maps produced for this study (Figs. 5 and 6) provide improved detail for evaluations of potential geothermal resources in the Appalachian Basin of New York and Pennsylvania. The average thermal gradient calculated from the data set is 23 °C/km with an average surface heat flow of 53 mW/m2. These values are similar to the Intergovernmental Panel on Climate Change (Fridleifsson et al., 2008) global average geothermal gradient between 25 °C/km to 30 °C/km and with the Pollack et al. (1993) global average heat flow for continental Paleozoic sedimentary rocks of 61 ± 30.2 mW/m2. On a regional basis the Appalachian Basin portion of Pennsylvania displays higher thermal gradient and heat flow than in New York State.

The average predicted geothermal gradient for the interpolated area of New York State is 23 °C/km, and the average predicted heat flow is 48 mW/m2. The precision in the estimation of any point within the interpolated area of New York (Fig. 5B) ranges from as low as ±0.5 °C/km in the areas of highly clustered data (as a result of n limited to 25) to as high as ±4 °C/km on the edges of the interpolation space where the data are sparse or nonexistent (n equaled only 20 sites within 53 km). The highest values of geothermal gradient and heat flow occur in several counties along the southern tier and central New York (Figs. 5 and 6); in these counties, the gradient ranges from 26 °C to 28 °C/km with an interpolation precision within ±1 °C/km (Fig. 5); the heat flow ranges from 50 to 60 mW/m2 with a precision within ±2 mW/m2 (Fig. 6).

The average predicted geothermal gradient for the interpolated region of Pennsylvania is 24 °C/km and the average heat flow is 60 mW/m2, slightly greater than in New York. Pennsylvania displays two northeast-trending zones of higher than average thermal values, one a broad swath through the central part of the state, and one along the western border of the state. In the western Pennsylvania zone, the geothermal gradient is 24–28 °C/km and heat flow is 55–75 mW/m2. Near the south end of the central Pennsylvania northeast-trending anomalous zone, the geothermal gradient is between 29 °C/km and 32 °C/km (Fig. 5A) and the heat flow is between 70 and 90 mW/m2. The interpolation precision within this area is generally ±1 °C/km for gradient (Fig. 5B) and ±2 mW/m2 for heat flow (Fig. 6B).


The northern Appalachian Basin in Pennsylvania and New York is filled with a thick series of relatively low conductivity Paleozoic strata above crystalline basement rocks. These low-conductivity sedimentary rocks lead to slightly increased temperature gradients in the upper crust. Results from a simplified analysis of the BHTs archived from oil and gas exploration in the basin allow a revised and more spatially resolved assessment of the geothermal resources in the basin.

This analysis focused on converting a widely available but low-quality data set, BHTs from thousands of wells, into maps of regional variability of the parameters of interest for a generalized geothermal resource assessment. The single-well thermal model estimates heat flow, but the accuracy of those predictions is not well constrained. Studies using similar model constraints (Frone et al., 2015; Erkan, 2015) determined the heat flow error to differ according to thermal conductivity values. The thermal conductivity error of one measured sample on a divided bar can be 5% (Frone et al., 2015) and, depending on if the sample was run wet or dry, as much as 25% (Erkan, 2015; Blackwell and Steele, 1989). For the Appalachian Basin of New York and Pennsylvania, a quality determination of the accuracy of the modeled temperatures at depth is not possible until more equilibrium temperature profiles at depths below 2000 m are acquired, to supplement the few equilibrium wells described in this paper (Fig. 3).

The North America geothermal map (Blackwell and Richards, 2004) is based on very few data for New York and Pennsylvania (a total of 50 points) and yet those data were known to be of high quality. That analysis estimated a maximum heat flow of ∼65–70 mW/m2 for New York and Pennsylvania, with those maxima located in south-central New York and western Pennsylvania. For the same regions as the maxima in Blackwell and Richards (2004), the new analysis estimates heat flow of 60–65 mW/m2 for south-central New York and 55–60 mW/m2 for western Pennsylvania. Most of the Appalachian Basin in New York and Pennsylvania was estimated by the data-sparse 2004 study to have heat flow values of 45–60 mW/m2; the new analysis predicts that values between 45 and 75 mW/m2 are widespread. The accuracy is not known for either estimate.

Even if the accuracy is unknown, the use of nearly 8000 BHT data and interpolation by localized kriging allowed for an improved representation of the spatial variability of geothermal resources as compared to previous studies in the region. For example, maps showing the standard error, a measure of interpolation precision, of the predicted geothermal resources communicate part of the uncertainty associated with assessing these resources. We can better distinguish between those regions with high temperature gradients that are statistically of high reliability due to numerous control points and those regions with equally high gradients but statistically low reliability.

Regions of highest heat flow are arrayed in an arc parallel to the Appalachian structural grain and occur within the proximal to medial sectors of the Appalachian Plateau in New York, and within the proximal and distal sectors of the Appalachian Plateau in Pennsylvania. The plateau (horizontally ruled area, Fig. 7) is a region of subdued thrust-related structural relief, detached regionally above the Silurian salt deposits (Jacoby and Dellwig, 1974). Although folds in localized parts of the proximal Appalachian Plateau have structural relief exceeding 200 m (cross-hachured areas, Fig. 7), over the great majority of the plateau area structural relief is much less than 200 m. In central Pennsylvania some of the most promising regions for further geothermal energy exploration are within the area of complex and high-relief folds. Because the fold-related geological complexity is not captured by the simplified one-dimensional thermal model, there is more uncertainty on the predicted temperature at each well in this strongly folded area than for sites farther west and north in the Appalachian Plateau. The thermally favorable areas of western Pennsylvania and south-central New York occur in the regions of subdued deformation, which may simplify thermal exploration.

In Shope (2012) it was noted that the southeastern flank of the northeast-trending Rome trough (Fig. 7) is a zone of relatively high geothermal gradient. This observation led to an examination of the possible impact of the Rome trough, a late Precambrian rift-like subbasin that underlies part of the plateau (Wilson, 2000; Ryder 2008), on the distribution of geothermal resources in the Appalachian Basin (Figs. 5 and 6). Neither the COSUNA stratigraphic data (AAPG, 1985) nor the AAPG Basement Rock Project Committee (1967) depth-to-basement map included the strata that compose the deep part of the Rome trough. It is possible that a variation in mantle or basement heat input beneath the Rome trough is responsible for an increase in the heat flow within the sedimentary basin strata (Shope, 2012). In Shope (2012) it was proposed that this hypothetical enhanced heat flow may generate hotter temperature in the deep horizons of basin fill within the trough than the values predicted in this paper (Fig. 7). However, testing of this hypothesis awaits additional well temperature data or more refined treatment of the crustal basement lithologies.

Even with the advances in knowledge of the thermal characteristics of New York and Pennsylvania presented here, there is room for improvement in the characterization of geothermal resources. Resource characterization can be improved through additional understanding of the caveats of the BHT data, and by utilizing borehole-specific conductivity and thickness data. In addition, inclusion of secondary geological information regarding faults and fracture networks may increase the precision of the estimates in kriging interpolations, especially in data-sparse locations.


The compilation and analysis of nearly 8000 temperature-depth data have led to a significant increase in the spatial resolution of geothermal resource assessment in New York and Pennsylvania. The analysis reveals that heat flow exceeds 50 mW/m2 in much of the Appalachian Basin of New York and Pennsylvania (Fig. 6). Areas of favorable geothermal gradient and surface heat flow represent areas with higher potential for heat extraction from either natural porosity and permeability or utilizing EGS. An area in south-central New York that has sufficient BHT data to permit construction of high-precision maps displays favorable geothermal resource potential, with heat flow estimates of 50–60 mW/m2. Two elongate, 200–300-km-long, northeast-trending bands of favorable geothermal resource potential occur in central and western Pennsylvania, with estimated heat flow of 55–90 mW/m2. For the western band the simplified thermal model based on one-dimensional heat flow is a good approximation of the subsurface geology, a condition that increases the likelihood that the thermal resource predictions of 55–75 mW/m2 are accurate. However, the highest heat flow coincides with a region of few wells and thus higher interpolation uncertainty (Fig. 6B). Additional caution is needed when considering the eastern band of apparent high heat flow because its large-amplitude folds increase the uncertainty on the results from the simplified thermal model and, in part of that region, the interpolations have low precision due to sparse BHT data (Fig. 6B).

A combination of the geothermal resource analysis presented here and a set of infrastructure and utilization economic models (Fox et al., 2011; Reber et al., 2014) suggests that heat from deep within the Appalachian Basin may be of sufficient quality (temperature) to provide direct thermal energy for district heating and hot water supply to residential and commercial buildings in New York and Pennsylvania. Where population centers have the potential for economically viable utilization of geothermal energy (Reber et al., 2014) and thermal gradient estimators are less precise, further exploration is warranted to increase the accuracy associated with the maps presented here of the geothermal resources. Two especially important foci for additional study include the utilization of the full well log data suites in areas denoted as favorable by this initial study and the acquisition of thermally equilibrated temperature logs in those areas.

This work was part of the 2010–2014 U.S. Department of Energy Project DE-EE0002852 (Geothermal Data Aggregation: Submission of Information into the National Geothermal Data System), led by Southern Methodist University (Dallas, Texas) with collaboration from Cornell University (Ithaca, New York), Siemens Corporation, the Bureau of Economic Geology at the University of Texas–Austin, the Geothermal Resource Council, MLKay Technologies, Texas Tech University (Lubbock), and the University of North Dakota (Grand Forks). Partial support was also provided by the National Science Foundation (NSF) Integrative Graduate Education and Research Training program (IGERT) through grant DGE-0966045 to Cornell University, NSF GK12 training grant DGE-1045513 (Grass roots: Advancing education in renewable energy and cleaner fuels through collaborative graduate fellow/teacher/grade-school student interactions), and the Cornell Atkinson Center for a Sustainable Future. We are grateful to the New York State Geological Survey and the Pennsylvania Topographic and Geologic Survey for access to well data. We thank two anonymous reviews and Associate Editor John Holbrook for their numerous constructive suggestions.


Thermal Model at Single Locations

There are three approaches to estimation of the equilibrium thermal field from bottom hole temperature (BHT) data. The first is based on modeling the thermal disturbance due to circulation of drilling fluids, but it is not viable for this study because in most cases the conditions leading up to and at the time of measurement are not documented. A second approach, which is physically robust, requires temperature readings for a given well at a succession of times following cessation of well drilling (Horner, 1951; Beardsmore and Cull, 2001). The latter approach cannot be applied to wells with a single BHT measurement (Deming, 1989), a common situation for the data in this study. A third empirical approach is widely used, followed in this study, is to estimate the nonequilibrium thermal deviation as a function of depth. If within a single field or basin, most wells were drilled in a similar fashion, through similar geological units and to similar depths, then if a correction to equilibrium conditions can be established for a small but representative set of wells, a corresponding correction can be applied to the entire data set.

An empirical BHT correction specific to the Anadarko Basin in Oklahoma was developed by Harrison et al. (1983), who compared borehole temperatures of wells considered to have reached thermal equilibrium to BHT values measured in recently mud-drilled wells. Harrison et al. (1983) interpreted that their data justified a correction term that increased with depth from 600 m to 3000 m, and that no useful relationship for greater depths was apparent. The simplified form of this correction (Equation A1) expresses ΔT (°C), the difference between equilibrium temperature and the observed temperature on a geophysical log at depth zw (m), as a second-order polynomial, 
(Blackwell and Richards, 2004). Gallardo and Blackwell (1999) analyzed the projection of the Harrison et al. (1983) correction to greater depth, and found a good fit if the correction below 3000 m depth is treated as a constant. The resultant ΔT value is a correction factor that can be added to the BHT from a geophysical log header to yield a corrected equilibrium BHT.

The degree to which an Anadarko Basin–based empirical correction is suited to the northern Appalachian Basin can be investigated by comparing a Harrison et al. (1983) corrected set of Appalachian Basin wells to temperature profiles for equilibrated wells located nearby (Spicer, 1964). In Frone and Blackwell (2010) it was shown that application of the Harrison et al. (1983) correction led to reasonably accurate results to a depth approaching 4000 m in the West Virginia sector of the Appalachian Basin. In Shope et al. (2012) corrected BHTs (Harrison et al., 1983) were compared to 14 equilibrium temperature wells in Pennsylvania and New York (Spicer, 1964); the maximum depth of the equilibrium temperature wells is <2400 m. In Shope et al. (2012) it was found that for the depth range of 1000–2400 m, the Harrison et al. (1983) correction (Equation A1) yielded a good fit to the equilibrium well profiles (Fig. 3). Given these published test results, our exclusion of wells shallower than 600 m, and a lack of constraints on accurate temperature measurements deeper than 4000 m anywhere in the northern Appalachian Basin, we applied the Harrison et al. (1983) correction to all of the BHT data.

The other steps in the thermal model were provided in Stutz et al. (2012). For calculation of an average thermal gradient (dT/dz) at the location of each data point (Equation 2), the vertical depth to the BHT measurement (zw) was assumed to be the lesser of either the logging depth, as measured by the well logger, or the true vertical depth, as reported by the driller. The value of TS (average annual surface temperature) was estimated to be 9 °C based on the Gass (1982) map of U.S. near surface (15–46 m) groundwater temperatures.

The average thermal conductivity of a column of rock depends on lithology, temperature, and pressure. Because temperature and pressure generally increase with depth, conductivity asymptotically approaches a constant value at sufficient depth. According to Blackwell et al. (2007), for sedimentary rocks this depth is at or near 4 km. Therefore, regardless of lithology, the rocks penetrated by any borehole at a depth >4 km within the sedimentary strata are assigned a constant conductivity, kds.

For shallower depths, conductivity values of various lithologies can be found elsewhere (e.g., Blackwell and Steele, 1989; Gallardo and Blackwell, 1999; Beach et al., 1987; Joyner, 1960). The thermal conductivity values adopted for the major lithologies in this study are 3.7 W m–1 °C–1 (sandstone), 2.8 W m–1 °C–1 (siltstone), 1.7 W m–1 °C–1 (shale), 2.7 W m–1 °C–1 (limestone), 4.3 W m–1 °C–1 (dolomite), and 5.6 W m–1 °C–1 (halite). Organic-rich shale is assigned a conductivity of 0.9 W m–1 °C–1 (Cercone et al., 1996). Because nearly all of the formations differentiated in the COSUNA (Correlation of Stratigraphic Units of North America; AAPG, 1985) columns are mixes of several lithologies, the assigned conductivity is an average weighted by the thickness of two principal lithologic constituents. For crystalline basement rocks, it may be appropriate to assume a single thermal conductivity for an entire borehole, forumla. However, sedimentary successions are typified by highly variable lithologies and justify a more detailed approach.

Utilizing the thickness of a lithologic unit (j) and its corresponding thermal conductivity, a thermal resistance, Rj (W–1 m2 °C) can be determined in which hj (m) is the lithologic unit thickness and kj (W m–1 °C–1) is the unit conductivity: 
The resistance for each lithology unit is then summed (ΣjRj) to calculate the total well resistance from the surface to the bottom of the well (or to an intermediate depth of interest).

Outlier Analysis

Prior to determining the thermal gradient field across the study region, local and global outlier analyses were conducted. The local analysis was performed using a 32 km by 32 km grid to define local regions. No outlier analysis was conducted for local grid boxes with fewer than 25 points. The global analysis was performed using the complete data set.

For a sample (i.e., local grid box data or global data set), the outlier identification criterion was defined as a measurement outside of the median ±3 times the upper quartile-to-median and lower quartile-to-median distances. Reported measurements outside of these distances in both the global and local analyses were removed from the data set for kriging interpolation. The data in local boxes with fewer than 25 points were also included in the kriging interpolation (for more information see Aguirre, 2014).

Kriging Interpolation

Kriging was used to compute a prediction and the standard error, a measure of precision of that prediction, for geothermal variables of interest. Kriging is a weighted interpolation using observed values and semivariogram functions describing the behavior of the random field (Matheron, 1963; Stein, 1999; Chilès and Delfiner, 1999, 2012). We assumed a constant but unknown mean for the geothermal variable in the extent of the estimation window, a condition often called ordinary kriging. The kriging model is presented in Equation A3, where G(x) is the unknown value of a geothermal variable of interest at location x (x1, x2), µ is the constant mean within the estimation window, ε(x) is the deviation of the geothermal variable’s value from the mean at location x, and δ(x) is the measurement or positioning error (e.g., depth in well or geographic position) associated with an observation at x. For this analysis a radial estimation window of 53 km was used to perform the kriging estimation. Because independent data on fluid circulation within the boreholes are lacking in many subregions of our study area, an additional assumption is needed. It was convenient to assume that the mean within the radial estimation window was constant because the data were gathered from the northern portion of the Appalachian Basin, a region over which extreme deviation from the mean was not expected (nor was a trend observed; Aguirre, 2014). Thus the model was 
Kriging predictions are linear functions of a data set of n values, as shown in Equation A4, where the kriging estimator of a geothermal variable Ĝ(x) at location x is a weighted average using weights wi on the observed values g(xi) at location xi: 
The selected set of points (x1,…xn) are the n points closest to x. The kriging weights wi sum to one, which is necessary for the estimator to be unbiased. In this study a maximum of n = 25 points was used to estimate the value Ĝ(x); more points could be used but additional precision was not needed in areas with reasonable densities of points. Furthermore, a minimum of at least 20 points within a distance of 53 km were required to estimate the value Ĝ(x). We used a maximum interpolation distance of 53 km to define when a sufficient local data base was available to estimate G at a point x (Aguirre, 2014).
Kriging algorithms can compute the variance of the expected mean, σ2EM, and the variance of prediction, σ2pred, for any point of interest. The two variances are defined in Equations A5 and A6. Kriging was performed under the assumption that the E{δ(x)} = 0 and Var{δ(x)} = σ2n, where σ2n is the nugget effect, which is attributed to BHT measurement error, error in the assumption of rock conductivities or lithologic thicknesses, positioning error, or other problems within the data set. Additional assumptions include that the variations in the true value of a geothermal variable were ε(x) and the variation was independent of local measurement errors, δ(x). 
The variance of the predicted mean describes how precise the estimate of the mean value of a geothermal variable is with the use of the available data. The variance of the prediction describes how accurately the value of a future measurement can be predicted. The variance of the predicted mean, σ2EM, is the error of interest because the mean value at an estimation location is the quantity of interest. For this assessment, the precision was limited by the use of a maximum of 25 points, which resulted in a lower bound on the precision of an estimation. For this study, the standard error (square root of σ2EM) was used as the measure of precision.

A semivariogram (Goovaerts, 1997) is used to model variations from the mean, and requires the adoption of an appropriate semivariogram function to describe the spatial relation of a geothermal variable of interest. In expectation, the difference between observations that are closer to one another is smaller than the difference between observations that are more distant from each other. The semivariance describes the behavior of those differences and thus increases with increasing distance between points. Often there is a distance beyond which no spatial dependence is observed; this distance is known as the range. At distances further than the range the semivariance remains constant at a value called the sill. For our BHT data set, the semivariogram does not go to zero at zero distance because of the nugget effect discussed here. The sill was defined as the nugget plus a partial sill, equal to the sill variance minus the nugget variance.

Directionality of spatial variability of a property can also be extracted from semivariogram analysis. The directionality is termed the anisotropy of the region and is defined by the angle, θ, of rotation from Cartesian axes, and the distance ratio, ϕ, of the direction of stronger dependence to the direction of lesser dependence (i.e., the ratio of the directional semivariogram ranges). Estimates of these parameters define the shape for the empirical semivariogram used in the kriging analysis.

The shape of the empirical semivariogram was determined by fitting a variety of classic and statistically permissible semivariogram functions to the data. The best fit was the Gaussian semivariogram model in Equation A8, where the semivariance assuming isotropic conditions at the Euclidian norm distance de = |d| = |xi – xj| between an observation of interest (i) and another observation at some distance away (j) is defined in Equation A7. 
In Equation A8 for the semivariance estimator, c is the partial sill, a is the range, and σ2n is the nugget. Figure 4 shows the fit to the geothermal gradient data provided by the Gaussian semivariogram. The fitted Gaussian semivariogram indicates that all points within 10 km have roughly the same semivariance, equal to the value of the nugget; the semivariance increases for point pairs separated by greater distances. The Gaussian fit matches the data throughout the kriging interpolation distance of 53 km.

Figures 5B and 6B illustrate the standard error of the geothermal gradient and heat flow, respectively. Areas with higher data density have lower standard error and higher precision.