Recent heat flow studies indicate that the Appalachian Basin in West Virginia may represent an important location for high heat flow and future geothermal energy development. Currently, however, only limited one-dimensional (1-D) heat flow studies exist in this region, making it difficult to assess the potential for geothermal development. Here, we develop the first high resolution 2-D basin model for a portion of West Virginia. The model uses 2-D finite difference heat conduction, basin cross sections, equilibrium temperature, and oil and gas bottom-hole temperature data to quantify heat flow at the surface and at the base of the sedimentary basin. The temperature data show elevated temperature gradients in the eastern portion of the basin. A 2-D advection-diffusion model, created using available lithologic and structural data, was designed to test whether variations in crustal properties, structure, erosion, or fluid advection can account for the observed temperatures in the basin. Thermal properties were populated using measured values as well as published averages. A linear heat flow vs. heat production relationship was used to determine heat flow at the base of the model. The model constrains the heat flow at the base of the sedimentary basin to 49–55 mW/m2. Analysis of modeling results suggests that heat flow at the base of the sedimentary basin is nearly uniform. Variations in basin temperatures are most likely due to variations in sediment thermal properties, complex structures, and/or localized fluid advection.
The thermal regime of the lithosphere is important for understanding formation and evolution of the Earth’s crust (Pollack, 1982; Nyblade and Pollack, 1993; Jaupart and Mareschal, 1999). Heat flow provides insight into the thermal regime of the lithosphere on a million year time scale. Stable continental igneous terrains have been shown to have an approximately linear relationship between measured surface heat flow and crustal heat production (Birch et al., 1968; Lachenbruch, 1968; Roy et al., 1968; Costain et al., 1986; Decker et al., 1988). This relationship implies that changes in surface heat flow within a stable crustal block are related primarily to lateral variation in average crustal heat production. There have been a limited number of heat flow studies in West Virginia (Joyner, 1960; Frone and Blackwell, 2010). This has resulted in a poor understanding of the heat flow in the basin and the nature of the thermal properties in the crust below it. Both of these properties have implication for oil and gas maturation and for geothermal energy utilization.
In this paper 6 equilibrium temperature logs and 148 corrected bottom-hole temperatures (BHT) are used to estimate basin temperatures and the surface heat flow along a transect of the West Virginia portion of the Appalachian Basin. Two temperature corrections are tested and compared against wells with equilibrated temperature logs to evaluate the validity of the corrections. Numerical modeling of a basin cross section is used to test the validity and potential causes for a previously observed heat flow anomaly, and to predict the equilibrium temperature distribution with the basin. The results provide insight into the nature of the crust below the basin, constraints on the surface heat flow, and a map of current basin temperatures.
The Appalachian Basin is a foreland basin located in the eastern United States (Fig. 1). The basin is ∼480,400 km2 in area and underlies portions of New York, Pennsylvania, eastern Ohio, western Maryland, West Virginia, eastern Kentucky, western Virginia, Tennessee, northwestern Georgia, and Alabama (Ryder, 2006). Sedimentary rocks contained in the basin range in age from early Cambrian (540 Ma) to early Permian (280 Ma) (Fig. 2). Three major orogenic events affected the sedimentation and structure in the basin during the Paleozoic Era. From oldest to youngest these are the Taconic orogeny (Ordovician), the Acadian orogeny (Devonian to Early Mississippian), and the Alleghanian orogeny (Mississippian to Permian) (Ryder, 2006). The Rome Trough forms the deepest part of the basin (Fig. 3). The trough is interpreted to be a major graben formed by passive margin extension during the early Cambrian along the east coast of the North American craton (Gao et al., 2000). It extends northeast from eastern Kentucky through West Virginia and western Pennsylvania. Grenville age (1.4 Ga) metasedimentary and igneous rocks underlie most of the basin (Shumaker and Wilson, 1996).
The West Virginia portion of the basin contains sedimentary rocks deposited from the Cambrian through the Pennsylvanian (Fig. 2). Cambrian and Ordovician lithologies are generally passive margin deposits consisting of marine sandstones and carbonates deposited in a shallow sea (Milici, 1996). Thicker sections of the Rose Run Sandstone and Black River–Trenton carbonates within the Rome Trough suggest that the trough was actively subsiding during the Cambrian and Ordovician (Patchen et al., 2006). The onset of Taconic orogeny sedimentation is marked by the Knox unconformity (Milici, 1996). The unconformity is the result of Middle Ordovician uplift and erosion of the Early Ordovician carbonate platform (Ryder, 2006). During the Taconic orogeny, the Blue Ridge crystalline basement rocks were thrust to the surface on the eastern margin of the Appalachian foreland. This resulted in erosion and widespread siliciclastic sediment deposition within the basin from the Late Ordovician to the early Silurian. Early Silurian sandstones and conglomerates were overlain by marine shales (Milici, 1996). Deposition of siliceous units related to the Taconic orogeny gave way to carbonate shelf and reef limestones and dolomites in the middle and late Silurian (Swezey, 2002). During the late Silurian the subsidence of the Salina salt basin resulted in deposition of as much as 400 m of evaporate minerals in the central portion of the basin (Milici, 1996). From the late Silurian into the Early Devonian, carbonate shelf deposits dominated deposition. The beginning of the Acadian orogeny in the Late Devonian resulted from oblique convergence of Laurentia with Avalonian continental fragments (Ettensohn, 1985). Devonian shales (Marcellus Shale) were deposited on the margins of the Catskill Delta (Ryder, 2006). Devonian black shales are overlain by Late Devonian–Middle Mississippian sandstones and interbedded sandstones and shales (Swezey, 2002). A period of carbonate platform deposition occurred in the Middle–Late Mississippian. In the Late Mississippian the Alleghanian orogeny resulted from the closing of the Iapetus Ocean basin and the assemblage of Pangea. This event is marked by a shift from marine limestones and shale sedimentation to coal-bearing siliciclastic units (Milici, 1996). Sandstone and interbedded sand and shales units eroded from the Alleghanian orogeny dominate the Pennsylvanian and Permian Periods (Swezey, 2002).
Basin Heat Flow
Heat flow is calculated from the product of the formation temperature gradient and thermal conductivity. Standard measurement requires an equilibrated borehole temperature log, with measurements typically taken every 1–10 m, and core or cuttings from the same well interval where the temperatures were measured (Blackwell and Steele, 1989; Blackwell and Spafford, 1987). Temperature gradients are calculated from the temperature log and the thermal conductivity values are measured from the well core or cuttings. Ideally, the wellbore is cased from surface to total depth and has been shut-in or not producing for an extended period of time to allow for near-field temperature disturbances cause by drilling and production to dissipate. This type of measurement can produce accurate results, but it can be difficult to obtain access to multiple wells within a given area. The difficulty in acquiring detailed equilibrium temperature logs has resulted in development of other methods of estimating heat flow. The same types of data are needed (temperature gradient and thermal conductivity), but typically they are of lower quality. While continuous or multiple temperature measurements in a single well are rare, a single temperature value in a well is very common. In petroleum wells a single maximum recorded temperature (or BHT) is typically measured as a part of the well logging process after drilling. This type of data is inherently noisy and results in large uncertainties in the temperature gradient (Blackwell and Steele, 1989). Despite the error associated with heat flow measurements derived from BHT, they are typically within 15% (Chapman et al., 1984). The large number of BHT data points that typically exist in sedimentary basins (often hundreds to thousands in oil and gas production regions) makes these data extremely valuable for broad scale [two-dimensional (2-D) and 3-D] heat flow studies in basins (Andrews-Speed et al., 1984; Chapman et al., 1984; Deming and Chapman, 1988).
The first heat flow study in West Virginia resulted in three heat flow values for the state (Joyner, 1960). These values range from 50 to 53 mW/m2 ± 25% for northern West Virginia. A more recent study of heat flow in West Virginia using BHT data shows a complex heat flow pattern across the state with an apparent high heat flow anomaly (∼75 mW/m2) in eastern West Virginia (Frone and Blackwell, 2010). This study used a generalized 1-D basin model for determining thermal conductivity values for individual wells. The results indicated heat flow values of 40–50 mW/m2 on the western edge of the state with values increasing eastward, with the highest values (70–75 mW/m2) generally east of the Rome Trough. These results are significantly different from previous heat flow values and maps (Joyner, 1960; Blackwell et al., 1990; Blackwell and Richards, 2004) that show a heat flow of 50–55 mW/m2 for practically the entire state. Continental heat flow globally is 65 ± 1.6 mW/m2 and 61 ± 1.2 mW/m2 for Paleozoic sedimentary and metamorphic rocks (Pollack et al., 1993).
Heat Flow–Heat Production Relationship
Two types of temperature data are used for this study: equilibrium temperatures and BHTs. Both data sets are from published and/or public data sources, therefore no new temperature measurements were completed for this study.
A subset of 15 temperature logs located in West Virginia from a national collection was used (Spicer, 1964). Most of the data date back to the 1920s and 1930s and were collected using maximum reading thermometers (Van Orstrand, 1935). One of the wells in this collection was drilled in 1891 and at the time was the second-deepest well in the world; Hallock (1897) provided some insight to the data collection and provided evidence that these data represent equilibrium temperatures. In addition, comparison of wells in the El Durado field in Kansas from Spicer (1964) and from McKenna and Blackwell (2001) show consistent results between the maximum reading thermometer data set and modern electric logging tools. The logs consist of temperature measurements collected at 250–500 m intervals over the depth of the wells. Well depths range from 1200 to 2200 m, and are located primarily in the northern and western portion of Kansas (Figs. 3 and 4).
Oil and gas BHT data were collected from various publications (American Association of Petroleum Geologists, 1994; Hendry et al., 1982) and from the West Virginia Economic and Geological Survey online database of well files. Data were limited to a minimum depth of 600 m for the study in order to minimize the effects of shallow groundwater flow and the effects of surface temperature changes on drilling mud temperatures and thermal gradient calculations (Blackwell et al., 2010). The final data set consisted of 1454 temperature-depth pairs, the locations of which are shown in Figure 3.
To investigate the potential causes for the observed heat flow and temperature variations in northern West Virginia, a 2-D steady-state temperature model was developed. The model was constructed using the detailed west-northwest stratigraphic cross section by Ryder et al. (2009) (A–A′, Fig. 3) and simulating heat flow advection-diffusion and the basement heat production. The goal of the model is to match equilibrium and BHT data along the cross section by varying the values for Qo, b, and A in the lithosphere. The final model is 228.5 km long and 9.2 km deep at 100 m resolution and includes basic facies information, well locations, faults, and other structures in the basin. The model was written in MATLAB (www.mathworks.com) using the finite-difference method.
The model heat production values incorporate the digital natural gamma logs used by Ryder et al. (2009) to construct the cross section. Formation top data were used to determine the rock type represented by the natural gamma log. The heat production results from Equation 4 were averaged for each rock type (Table 1).
Temperature values within a distance of 25 km from either side of the cross section were projected onto the 2-D cross section (red and blue dots in Fig. 3). This resulted in 143 BHT data points and 6 equilibrium temperature logs (Figs. 4 and 5).
For the model, the crust below the sedimentary cover had a constant thickness that was defined by b in Equation 1. This is similar to the way the lithosphere below foreland basins is modeled as an elastically deformed plate (Flemings and Jordan, 1989; Garcia-Castellanos et al., 1997). This model assumes constant heat flow and heat production from the basement rocks. An additional model was tested where the crust was thinned based on the thickness of overlying sedimentary formations. In this model, the heat flow entering the basin varied and results in lower heat flow and temperatures in the deepest part of the basin (Rome Trough). The data do not show a decrease in temperatures in the Rome Trough area, so this model of the crust was not examined further. The model used for the crust is a silicic upper crust and a mafic lower crust (Drury, 1989; Christensen and Mooney, 1995). The terms mafic and silicic in this case refer only to the relative amount of heat producing elements found in each layer, with the mafic layer having lower and the silicic having higher heat production. Each model iteration was run to steady state with varied input values for reduced heat flow (Qo), upper crust heat production (A), and initial upper crustal thickness (b). The heat flow at the sediment-basement contact was constant across the model due to the constant values used for Qo, A, and b.
Using constraints from Thakur (2011) for Qo (33 ± 1 mW/m2) and b (7.5 ± 0.2 km) a range of models was tested with different values for A. The basement A value was varied between 1.5 and 3.5 μW/m3, the range for granitoid rocks (Vilà et al., 2010).
The model results showed that for the constrained Qo and b values for the eastern United States from Thakur (2011) and basement heat production values between 2.0 and 3.0 µW/m3, the model reproduced the observed temperatures within this section of the basin (Figs. 5A and 6). Heat flow at the sediment-basement contact ranged from 49 to 55 mW/m2. This heat flow constraint can be further used to constrain future crustal studies in the region. In the model this value was controlled by the input values for Qo, b, and A, but different distributions and values of heat production and reduced heat flow can also result in the same heat flow at the base of the sedimentary basin. The range in heat production values is consistent with granitoid rock types (Vilà et al., 2010). The model shows increasing surface heat flow from west to east. The thicker sedimentary package occurring on the eastern side of the model results in more heat production and therefore higher surface heat flow.
The surface heat flow calculated along the profile is compared with the heat flow from Frone and Blackwell (2010), who used primarily BHT data and a simplified 1-D stratigraphic model to determine heat flow within the study area (Fig. 6). The region of the highest heat flow values from Frone and Blackwell (2010) correspond with the largest temperature residual range from the 2-D model.
In the first scenario, there is higher heat flow due to higher heat production in the basement or a thicker upper crust beneath the eastern side of the basin. From Equation 1, it is clear that if Qo is constant over large areas, then changes in Q are related to changes in b, A, or a combination of the two. For example, for a 10-km-thick upper crust, a uniform increase in heat production of 1 μW/m3 results in a surface heat flow increase of 10 mW/m2 (Fig. 7). The same increase in heat flow could also be generated from a constant heat production and an increase of the upper crustal thickness from 10 to 20 km (Fig. 7). From preliminary receiver function analysis of the Earthscope Transportable Array over West Virginia, Pennsylvania, Ohio, Virginia, and Maryland, crustal thickness within the Appalachian Basin is 45 ± 5 km (Crotwell and Owens, 2005; Trabant et al., 2012). Therefore, in this study variations in heat flow are interpreted as predominantly changes in basement heat production.
Magnetic anomalies suggest that the eastern boundary of the Rome Trough is part of a large-scale strike-slip fault (King and Zietz, 1978; Steltenpohl et al., 2010) or a Grenville-aged suture between Laurentia and an accreted terrain (Tohver et al., 2004; Mueller et al., 2008). Steltenpohl et al. (2010) suggested a right-lateral offset of ∼220 km that is continuous from New York to Alabama. Both of these interpretations imply that a different basement lithology may be present under eastern West Virginia. Heat flow results for southwestern Pennsylvania from Shope (2012) show elevated heat flow values within the Rome Trough, ∼200 km north northeast of the modeled cross section, that is suggested to be the result of elevated basement heat production. This result is consistent with the interpretation of the magnetic data. This case was tested by allowing the heat production in the basement to vary at the eastern boundary of the Rome Trough by 25%–200%. The results of this modeling show that an increase in heat production of 25%–75% in the eastern basement block still fits the observed data, but does not produce the higher temperatures observed. Higher heat production values can fit the higher observed BHT values; however, the models no longer fit the equilibrium temperature data in wells WV-3 and WV-8 and this does not explain the lower BHT values in the same area. For heat production variations in the basement to explain the higher observed BHT values between these wells, the eastern basement block would need to be divided into subblocks. There is evidence for down to the east normal faulting in the eastern basement block to the south (Ryder et al., 2008); however, the offset of these faults decreases to the north and is thought to terminate ∼25 km north of the modeled cross section (Shumaker and Wilson, 1996). Higher heat production values in smaller subblocks were tested to attempt to explain the elevated BHT values. Results show that an increase in heat production of 50%–70% from the eastern margin of the Rome Trough to the normal fault shown by Shumaker and Wilson (1996) provides a better fit to the elevated BHT values; however, this model also produces higher than observed temperatures in WV-3 and WV-8 and does not explain the low BHT values. Thus first-order changes in heat production in the basement do not provide a simple explanation for the variation in BHT values between WV-3 and WV-8.
Lateral Groundwater Flow
The second scenario is that intrabasin fluid flow could result in the observed temperatures. This case would imply that the thermal anomaly is caused by advective heat transport rather than conduction. A large-scale example of advective heat transport is the heat flow anomaly in South Dakota and Nebraska (Gosnold, 1990). That anomaly is caused by gravitationally driven groundwater flow, laterally transporting heat within the sediments. It has been shown that velocities of 10–8 m/s were needed to produce the observed heat flow anomaly. Darcy velocities for the Appalachian Basin from Gupta and Bair (1997) are used to model the effects of advection on the basin temperatures. Using their maximum velocity value of 3.05 × 10–12 m/s, the maximum temperature change within the basin is ∼1 °C compared to pure conduction. The temperature effect of advection at this velocity is small and the change at any one well is within the error of the measurement. Advection may still cause localized temperature effects; however, basin-scale advection as currently understood cannot reproduce the observed temperatures.
The third scenario of the anomalous BHT values is related to the structures in the basin. The anomalous values are generally located between the Etam and Deer Park–Leadmine anticlines. The high-amplitude anticlines are the result of thin-skinned deformation during the Late Mississippian–Permian Alleghanian orogeny (Wilson and Shumaker, 1992; Kulander and Ryder, 2005). The anticlines were formed by compression and deformation along décollement horizons. The results of this thrusting are high-amplitude thrust fault–cored anticlines (Ryder et al., 2009). Structures like these can affect temperatures in two ways. Isotherms can be either elevated or depressed at the peak of anticlines depending on the thermal conductivity of the deformed formations (Guyod, 1946). Elevated isotherms are within the anticlines in this model (Fig. 5). The equilibrium well WV-8 is located on the crest of the Deer Park–Leadmine anticline and shows elevated temperatures compared to the other equilibrium wells (Fig. 5). A few of the BHT values in this region between WV-3 and WV-8 show the largest deviations from the modeled values.
The second temperature effect that could be caused by the structure is localized fluid flow along thrust faults within the anticlines. Bonini (2007) showed that for thrust fault–related folds, a reduction in lithostatic pressure can cause fluid overpressure and migration along faults. The reduction in lithostatic pressure is caused by surface erosion. The migration of fluids from depth along faults results in elevated temperatures along the faults. The magnitude of the temperature effect in the rock surrounding the fault is dependent upon the temperature of the fluid and the time scale of the flow. Similarly, migration of fluids from the surface downward may result in depressed temperatures. The scale of the temperature effect is dependent upon the fluid volume and velocity. The migration of fluids along faults in the anticlines is likely not uniform along the strike of the structure due to heterogeneities in the subsurface. Advection of fluids along fault pathways is not included in this model. This is an unlikely process in the Appalachian Basin for two reasons. First, slip along a fault plane can result in the formation of phyllosilicate minerals, which in turn decrease the permeability of the fault (Cavailhes et al., 2013). Second, the lack of recorded seismicity in the region suggests that faults have not been reactivated, resulting in sealed faults that are likely to inhibit fluid migration.
The fourth scenario for the observed temperature scatter is related to erosion or sedimentation rates with in the basin. Erosion of the upper surface of a basin results in an uplift of isotherms in the subsurface while sedimentation depresses isotherms, resulting in a higher or lower heat flow, respectively (Jaeger, 1965). This results in a high or low temperature gradient. This gradient will persist if the rate of erosion or deposition is high enough to outpace the dominant heat transport method. The degree to which erosion or sedimentation affects the heat flow is controlled by the thermal diffusivity and the rate of erosion or sedimentation.
From fission track studies (Roden, 1991; Blackmer et al., 1994; Boettcher and Milliken, 1994) and radiogenic isotope studies (Matmon et al., 2003), erosion rates for the Appalachian Basin from 265 to 30 Ma are between 10 and 28 m/m.y. From the Miocene to the present, erosion rates increased to 50–100 m/m.y., possibly due to changes in mantle dynamics (Miller et al., 2013; Blackmer et al., 1994; Boettcher and Milliken, 1994; Galloway et al., 2011). In order to constrain the maximum change in temperatures that could be caused by erosion, a rate of 100 m/m.y. for 20 Ma was used.
In this study a numerical solution was used to estimate the effect of erosion on basin temperatures. The solution instantaneously erodes 100 m/m.y. of rock. The solution shows a decreasing effect with depth, with a maximum temperature increase of 6.8% at the surface and an average temperature increase of 6.7% in the depth range 1–6 km above a steady-state background temperature distribution. This is for a 1-D case, which is preliminarily correct for the study area, so that the temperature increases affect all wells in the study more or less equally (Blackmer et al., 1994). For this portion, and likely all, of the Appalachian Basin, erosion has an effect on subsurface temperatures and heat flow; however, it does not appear to be a viable solution to the variable measured temperatures.
Previous heat flow studies in the Appalachian Basin have had limited available borehole data (Blackwell and Richards, 2004; Joyner, 1960). Use of BHT data can provide a large data set that, when checked against equilibrium data, can be a useful source of temperature data. Incorporating both equilibrium and BHT data with heat flow and heat production provides a more detailed analysis of heat flow and basin temperatures in the northern West Virginia portion of the Appalachian Basin.
This study shows that the thermal regime in this portion of the Appalachian Basin is dominated by conductive heat flow. A uniform thickness model of the crust following the Q-A relationship can explain the heat flow and heat production in the crust below the basin. Results show that for Qo = 33 ± 1 mW/m2 and b = 7.5 ± 0.5 km the heat production in the basement is between 2 and 3 µW/m3 and the heat flow at the base of the basin is between 49 and 55 mW/m2. The heat flow value at the base of the basin had the strongest influence on the modeled temperatures. Therefore, the values for Qo, b, and A could vary outside of the ranges used from Thakur (2011) as long as the heat flow at the top of the basement is within the range of 49–55 mW/m2.
Models using the Q-A relationship for heat flow within the crust show that most of the measured temperatures within the basin can be matched with a simple conductive heat flow model with uniform basement properties. Variation in basement heat production values to 10% does not significantly affect the results of the model. The modeled surface heat flow values between 55 and 63 mW/m2 and are within the expected range for Paleozoic age sedimentary rocks based on the Pollack et al. (1993) analysis.
Heat flow calculations from Frone and Blackwell (2010) using BHT data show that the heat flow east of the Rome Trough is ∼40% higher than within and west of the trough. Numerical models of a 2-D cross section of the basin are consistent with a heat flow increase from west to east; however, the predicted increase is closer to 15% and is more variable. Thus, the heat flow anomaly that was previously reported appears to be overestimated. The values may have been overestimated due to variation in basement heat production, large-scale basin advection, the inclusion of BHT measurements that were affected by local advection caused by fluid migration along thrust faults within anticlines, and erosion. Of these potential causes, complex structure and 3-D effects best explain the variability in the BHT values and the calculated heat flow values.
We thank Associate Editor Anna Crowell, Paul Morgan, and two anonymous reviewers whose questions and comments helped to improve this manuscript. The Roy M. Huffington Department of Earth Sciences and the Southern Methodist University Geothermal Laboratory provided funding for this research.