We use the Curie depth derived from spectral analysis of near-surface magnetic anomaly data to constrain the solution of a one-dimensional (1-D) steady-state heat-flow equation. The method, in addition to anchoring the geotherm at deep levels in the crust, yields the ratio of radiogenic heat production to thermal conductivity. We evaluate the utility of this constraint in two granitic regions in the United States where radiogenic heat production, thermal conductivity, and heat flow are well determined. We also examine the results of this method in the context of previous research in the central Red Sea and coastal Saudi Arabia. In a test case using data from New Hampshire, USA, the steady-state approximation applies, and we calculate heat-production values (∼8 μW/m3) that agree with measurements on samples of Conway granite in the region. In the second example, our Wyoming geotherm is hotter by ∼200 °C at the Curie depth, possibly reflecting the 15 km thicker crust (as observed from the dense EarthScope USArray stations in comparison to only two earlier refraction profiles) and the consequent greater radiogenic heat production in the lower crust. In the Colorado Front Range, our radiogenic heat-production value is somewhat higher than predicted from previous heat-flow studies and may reflect heat advected by intrusions. Our results indicate that basal crustal temperatures in northern Colorado may be close to solidus of rocks of felsic composition, a scenario that is consistent with the geological history of the region and has important implications for crustal rheology. Our estimates of depth-integrated heat production from the crustal columns (25 and 55 mW/m2 Wyoming and Colorado, respectively) also agree well with estimates from previous studies.

Solutions of the 1-D steady-state heat-flow equation with different boundary conditions (e.g., heat flow at the surface or the Moho) are commonly used to construct geotherms in the lithosphere in situations where steady-state conductive thermal transfer is assumed (Carslaw and Jaeger, 1986; Fowler, 2005). In addition to heat coming from the mantle, heat produced by decay of radioactive elements, primarily U, Th, and K, in the crust is an important contributor to heat flow at the surface (Clauser, 2011). This is true especially for continental crust where felsic rocks are enriched in the heat-producing radioactive elements (Wollenberg and Smith, 1987). Unless the contribution from radiogenic heat is known from measurements of rocks, it must be estimated, and unless the deep geotherm is constrained petrologically from xenoliths (which may not be representative of present-day temperatures), temperature deep in the crust is largely unconstrained. Where there is no temperature constraint at depth and/or where the assumptions adopted to solve the heat-flow equation are violated, the derived geotherm calculation may misrepresent temperatures deep in the crust.

One parameter needed to estimate the geotherm is radiogenic heat production of the crust. Based on measurements and geochemistry of crustal rocks, the average rate of radiogenic heat production in the upper continental crust is estimated between 2.5 and 2.9 μW/m3. The range for felsic rocks is 0.7–28 with the mean of 4.8 μW/m3 and the range for mafic rocks is 0.04–4.1 with the mean of 0.7 μW/m3 (Wollenberg and Smith, 1987). The continental crust is heterogeneous, however, and in all but a few special geologic terrains (large granitic plutons), surface measurements are not representative of upper crustal radiogenic heat production, and the heat production rate for any specific locality is not known. In this paper, we develop a method to estimate heat production of the crust by applying a constraint of a temperature at depth (specifically, the Curie temperature depth of magnetite inferred from spectral analysis of magnetic data; although other temperature proxies from seismology and electromagnetic methods may be used similarly).

It is well known that the base of magnetization in the lithosphere may be controlled by the phenomenon of Curie temperature of magnetic minerals. Curie temperature, or Curie point, is the temperature above which ferromagnetic materials lose their spontaneous magnetization and become paramagnetic. At depths where the temperature is above the Curie temperature of magnetic minerals, rocks are paramagnetic and do not contribute to magnetic anomalies observed at or above the Earth’s surface. This phenomenon has been used in deriving the Curie point depth in several studies since Bhattacharyya and Leu (1975) developed a method for calculating this depth. The magnetization of the crust is often governed by concentrations of magnetite because generally magnetite is more strongly magnetic than most other minerals. Rocks containing hemo-ilmenite or ilmeno-hematite have been shown to carry a significant amount of natural remanent magnetization (McEnroe and Brown, 2000; Robinson et al., 2002) and will have variable Curie points depending on composition. However, these minerals occur only in certain environments, and, hence, the Curie point of magnetite (580 °C) is generally used to estimate temperatures deeper in the crust from the detected base of magnetization.

Early magnetic depth-determination methods assumed sources and ensembles of sources of uniform magnetization in the crust, and the slopes of the annular-averaged spectra from magnetic anomalies were used to estimate depths to magnetic layers (e.g., Spector and Grant, 1970; Okubo et al., 1985; Tanaka et al., 1999). In the 1990s, researchers recognized that the crustal magnetization may be fractal in its distribution (Maus and Dimri, 1994, 1995; Pilkington and Todoeschuck, 1995). A fractal distribution introduces a power-law dependence into magnetic spectra (Pilkington and Todoeschuck, 1993). Without correction for the fractal effect, the depths to magnetic layers and the magnetic bottom can be overestimated. In determining magnetic bottom depth, Ross et al. (2006) and Ravat et al. (2007) introduced the method of forward modeling the spectral peak caused by the bottom of the magnetic layer. However, most magnetic spectra do not form spectral peaks; the lack of a peak is considered to reflect fractal behavior of the magnetic anomaly field and fractal behavior of magnetization. Around this time, methods that consider the fractal nature of crustal magnetization were introduced in magnetic depth-estimation methods (Maus and Dimri, 1996; Maus et al., 1997; Ravat et al., 2007; Bouligand et al., 2009; Bansal et al., 2011). However, until Salem et al. (2014) introduced the de-fractal method, it was not possible to estimate the fractal parameter of the magnetic field needed in applying the correction. Most researchers use a fixed value to determine magnetic bottom depths (e.g., Bouligand et al., 2009; Bansal et al., 2011, 2013; Li et al., 2013); however, the a priori fixed fractal parameter gives correct depths only when the chosen fractal parameter is correct (and if spectral slopes are picked correctly). The method will under- or over-estimate depths if the chosen fractal parameter is higher or lower, respectively. The de-fractal method of Salem et al. (2014) compensates for the fractal parameter of the magnetic anomaly field such that a spectral peak is formed, which then can be modeled using forward-modeling techniques (Ross et al., 2006; Ravat et al., 2007). The de-fractal method also uses the consistency of results from both the centroid method (Okubo et al., 1985; Tanaka et al., 1999; Bansal et al., 2011) and the peak modeling method to identify the fractal correction parameter. The 2-D de-fractal method (where magnetization variation is in the xy-plane and magnetization is constant in depth) is developed and tested in detail by Salem et al. (2014). We modified the procedure of Salem et al. (2014) and made it semi-automatic for the purpose of testing additional feasible interpretations. The details of the procedure and relevant model study results are given in  Appendix A.

In this paper, we use magnetic bottom depths derived from the de-fractal method and the Curie temperature of magnetite as a temperature-depth constraint in the solution of a 1-D steady-state heat-flow equation. We find that the constraint allows us also to estimate the bulk radiogenic heat production. This analysis assumes an exponential distribution of heat production in the crust, and justification for this assumption is given in  Appendix B. We test this method in two areas in the United States with well-constrained surface heat-flow and radiogenic heat-production data; we also reassess our previous results in the Red Sea region and in coastal Saudi Arabia (Salem et al., 2014).

Geotherm Using the Curie-Depth Constraint

We assume a linear relationship of surface heat flow to heat production (Birch et al., 1968; Roy et al., 1968):
where qs is surface heat flow, qr is reduced heat flow, As is surface heat production, and b is the empirical heat-production depth-distribution parameter. The linear relationship is consistent with an exponential depth distribution of heat production (Lachenbruch, 1968):

The linear relation and exponential distribution were introduced for terrains dominated by granitic plutons, and we apply them in the same types of terrains in this paper. The exponential distribution has been challenged, and in  Appendix B, we discuss data that have been used to question its validity. This distribution is used in the present study only for terrains dominated by granitic plutons.

Assuming as boundary conditions,

  •  (i) T = Ts at z = 0 (the surface), where T is temperature, z is depth, and

  • (ii) q = –K dT/dz = –qs at z = 0, where K is thermal conductivity, q is heat flow, and qs is the surface heat flow and neglecting contributions from mass advection, transient cooling or temperature dependence of thermal conductivity, the solution of the 1-D heat-flow equation is:

We add an a posteriori condition to Equation 1: T = Tc at z = zc, where Tc is the Curie temperature and zc is the Curie depth.

We restrict the application of this technique to provinces dominated by granitic plutons from which a mean qs and b are known. Using forumlaEquation 1 becomes

The ratio As/K (Equation 2) is an important parameter calculated in this study. The parameters in Equation 2 are derived as follows: TC is the Curie temperature, which is taken as 580 °C (TC of magnetite); TS is the mean annual surface temperature, which is based on local climate stations and/or latitude and elevation (small errors in TS are insignificant with respect to the uncertainty in TC based on exact magnetic mineralogy); zC is Curie depth, calculated from the spectral analysis of magnetic anomaly data (see  Appendix A); qs = K(dT/dZ)z=0; and b is the depth parameter from the qs-As relation. Substituting As/K and surface (dT/dz)z= 0 or As and qs = K(dT/dZ)z=0 in Equation 1, with the geotherm anchored at T = Tc at z = zc, Tz can be determined for all depths zzc.

If the analysis is restricted to provinces with major silicic plutons, as in this paper, then a reasonable assumption is that the upper crust is silicic in composition. All studies that make calculations of crustal geotherms must make assumptions of crustal thermal conductivities. Fortunately, the greatest variation in thermal conductivity is in quartz-bearing rocks in the uppermost crust. As the crust probably converges to a more intermediate and/or mafic composition with depth, and as crustal thermal conductivity decreases with increasing temperature, reasonable estimates of lower crustal conductivity may be made (Roy et al., 1981; Clauser and Huenges, 1995; Kukkonen et al., 1999). Thus, with a reasonable model of an average K for the crust, As may be calculated. Alternatively, if a reliable regional surface heat qs is known, K may be determined from qs = K(dT/dz).

Application of the Method to Two Tectonically Diverse Regions in the United States with Radiogenic Heat-Production Data

In order to compare the results of our method and observations, we have identified two tectonically diverse regions in the United States where contributions from geology and/or tectonics may be deduced from geological and geophysical data sets. One is in New Hampshire, in the northeastern United States, where highly radiogenic granite batholiths in the upper crust contribute significantly to the geotherm and the temperatures are likely to be in steady state; the second is a boundary between the Archean Wyoming craton and tectonically younger northern Colorado.

New Hampshire

In central New Hampshire, we used the de-fractal method on the long-wavelength–corrected North American magnetic anomaly compilation known as NURE_NAMAM2008 (Ravat et al., 2009). Robustness of long wavelengths is important for the interpretation of deeper magnetic sources including the interpretation of Curie depths. We made three independent estimates of magnetic bottom (ranging from 25 to 32 km). Surface heat flow and its relationship to radiogenic heat production in the eastern United States is well established (Roy et al., 1968); the characteristic depth of heat-producing elements deduced from this is 7.5 km, and the reduced heat flow is 33 mW/m2. In the Southern Methodist University heat-flow and/or thermal properties database (Blackwell and Richards, 2004), three heat-flow determinations in the region of the Conway granite are 79, 90, and 95 mW/m2, and averages of hundreds of heat-production measurements from each of these locations are 7.36, 8.87, and 8.66 μW/m3, respectively (Roy et al., 1968). The Conway granite covers a large area; it was emplaced during Mesozoic separation of North America from Europe (McHone and Butler, 1984) and is extremely rich in thorium (56 ± 6 ppm; Adams et al., 1962). Mean thermal conductivities of samples of rocks corresponding to the three heat-flow locations are 3.27, 3.56, and 3.7 W m–1 K–1 (Roy et al., 1968).

Following the details of the method shown in  Appendix A, we show the de-fractal method’s fractal parameter-corrected and modeled spectra for the New Hampshire example in Figure 1. The method relies on the consistency of magnetic top and bottom estimates from the forward-modeled and the de-fractaled centroid methods (see Salem et al., 2014). We used the Curie-constrained geotherm determination method on the New Hampshire heat-flow data, and the results are summarized in Table 1. We calculated unknown parameters (the ratio As/K) for a range of thermal conductivity values and selected values that closely match the known reduced heat flow in the region (Roy et al., 1968). Identifying K in this manner also determines As. These results show that where heat flow is steady state and the crustal heat production is relatively simple, observed average surface heat production and the calculated surface heat production (exponentially decreasing with depth) agree well.

The example shown in Figure 1 yielded a base of magnetization or Curie-depth estimate of 30 km. The Curie-depth–constrained geotherm and the geotherm computed using the observed radiogenic heat production from the Conway granite are similar and within the error bounds of the range of parameters used to construct them, as shown in Figure 2. More importantly, had one erroneously assumed a typical surface heat-production value of ∼2 μW/m3, which is often done in regions where heat production is not known, it would result in much lower temperatures in the crust and incorrect geophysical and geodynamical inferences (see continuous black geotherm with As of 2.1 μW/m3 and K of 2.5 W m–1 K–1 in Fig. 2).

The derived high surface heat-production estimates give an approximate depth to the solidus of ∼97–110 km (with mantle thermal conductivity value of 3.4 W m–1 K–1), which is identical to the estimates of seismic lithosphere-asthenosphere boundary (LAB) from P to S (Ps) phase conversions at 97–110 km for this region (Rychert et al., 2005). However, the P to S (Ps) phase conversions from the nearby stations from which the interface was interpreted are the least clear and weakest among the data they had analyzed; and Rychert et al. (2007) suggest that it could be shallower than their 2005 study estimate. Regional anisotropy of seismic waves in the mantle and steep-dipping structure (Menke and Levin, 2002) could also introduce errors in the determinations. Nonetheless, the Curie depths derived in the study are consistent with seismic LAB determinations. Seismic LAB is controlled by factors such as composition, anisotropy, and water content, which are not identical to the factor that controls the thermal LAB (i.e., the solidus at a particular composition).

Northern Colorado to Wyoming Craton

The Wyoming block is an Archean cratonic core to which younger terranes were accreted (Whitmeyer and Karlstrom, 2007). Decker et al. (1988) have collected, compiled, and analyzed heat-flow, radiogenic heat generation, and conductivity measurements in Colorado and southern Wyoming. There is much variability over short distances in heat-flow observations in Colorado (64–140 mW/m2; Decker et al., 1988, their table 3), and Decker et al. (1988) showed that the variability was related to upper crustal radiogenic heat-production variations (reaching values of up to 8–9 μW/m3). In contrast, the Wyoming craton consists of once middle and lower crustal rocks lying presently in the upper part of the crust and thus explaining representative low radiogenic heat generation (1–2 μW/m3) and low heat flow (<50 mW/m2), except where there are rocks of granitic composition (Decker et al., 1988). Decker et al. (1988) synthesized their thermal data with interpretations of upper crustal granitic compositions from gravity lows (Brinkworth et al., 1968; Tweto and Case, 1972; Isaacson and Smithson, 1976; Johnson et al., 1984; Decker et al., 1988) to remove the upper crustal heat-production variability and called the resulting heat flow the residual heat flow. Using finite difference computations, they constructed crustal models of thermal parameters. Their models suggested near-solidus conditions for felsic compositions (∼800 °C) near the bottom of the thickened crust in Colorado (Prodehl and Pakiser, 1980). The average upper crustal radiogenic heat generation for their models was 3.8 and 1–2.4 μW/m3 for Colorado and Wyoming, respectively. In Decker et al.’s (1988) model, which is based on two seismic-refraction profiles between 105°W and 107°W longitudes (Jackson and Pakiser, 1965; Prodehl and Pakiser, 1980; see also Keller et al., 1998), a 10 km crustal thinning occurs from Colorado into Wyoming. However, joint analyses of the EarthScope seismic receiver function and gravity anomaly data show that the crustal thickness in both the regions is uniformly 45–55 km, and the average crustal Vp/Vs varies between 1.7 and 1.8, indicative of the bulk crustal composition from granites to diorites (Lowry and Pérez-Gussinyé, 2011). Similar high crustal thickness estimates in the Wyoming craton have also been derived in the recent work by Schmandt et al. (2015; digital results also presented in IRIS DMC, 2016). How these differences in Moho depth will influence the thermal model is difficult to judge because it depends on the distribution of heat-producing elements in the crust and whether the radiogenic elements extend all the way to the Moho or not. However, crossing the solidi alone in some regions would complicate the situation and modeling. To circumvent this issue, we compare the depth-integrated radiogenic heat production (which is the contribution of crustal radiogenic heat production to surface heat flow) between the model of Decker et al. (1988) and our results. Any advective heat or near-solidus temperatures in the affected parts of the crust will also reduce thermal conductivity (below 2 W m–1 K–1 at 800–1000 °C, Clauser and Huenges, 1995) and could reduce thermal conductivity for deep-seated rocks in the crust in comparison to the steady-state scenario modeled by Decker et al. (1988).

We calculated a large number of Curie depths at 125 km spacing in Colorado and the Wyoming craton (Fig. 3). Shallow Curie depths in this map generally correlate well with the areas of high surface heat flow (except where our 500 km spectral windows cannot capture the variation in the narrow zone in the Snake River Plain). The import of this regional correlation is that the spectral magnetic bottom estimates have no input from the heat flow or any other geophysical or geologic observations. The Curie depths in Colorado range between 22 and 32 km, with an average of ∼26.5 km. Estimates in the Wyoming craton range between 24 and 52 km, with an average of 34 km. We calculated geotherms using the Curie-depth constraints at the end points of the crustal thermal parameters model profile of Decker et al. (1988), corresponding to the Laramie Range with Curie depth of ∼30–33 km and surface residual heat flow of 45 mW/m2, and the Colorado Front Range with Curie depth ∼24–25 km, surface residual heat flow 85 mW/m2. Using a radiogenic depth parameter, b, of 10 km, we first estimated the ratio As/K from Equation 2 and then varied the thermal conductivity to match the reduced heat-flow values of Decker et al. (1988) for the two provinces (30 and 20 mW/m2, respectively, for Colorado and Wyoming). The results of these calculations are presented in Table 2.

Decker et al. (1988) used a three-layer crust model with different thermal parameters in each layer and a crustal boundary close to the Colorado-Wyoming state line (Fig. 4); in the figure, we show also the comparison between their and our geotherms. They did not have the benefit of a Curie-depth estimate in constructing their thermal model. For the Laramie Range region in Wyoming, our analysis yields temperatures in the lower crust ∼200 °C hotter at the Curie depth than the Decker et al. (1988) models. The discrepancy between the two approaches arises primarily from two differences: First, Decker et al.’s (1988) thermal models have high middle and lower crustal radiogenic heat-production values and make the equivalent surface heat production in the exponentially decreasing formulation higher. The second difference is that the Wyoming craton crust is thicker than assumed by Decker et al. (1988) as previously discussed. Thus, it is possible the amount of radioactive elements may be greater in both these regions than estimated by Decker et al. (1988). Our (exponentially decreasing) surface radiogenic heat production in the Colorado Front Range is higher than the constant heat production of the upper layer in Decker et al. (1988), but our high values are compensating for high heat production of the middle crustal layer in the Front Range of Colorado used by Decker et al. (1988) (1 μW/m3 higher than the corresponding layer in Wyoming). In the Wyoming craton, our magnetic bottom estimates (∼34 km) are shallower than the Moho (except near the northeast edges of the craton), and therefore this magnetic bottom may correspond to the Curie depth. The mantle is often believed to be non-magnetic, and, if the Curie isotherm falls below the Moho, the bottom of magnetic sources is assumed to be the Moho (Wasilewski et al., 1979). The depth-integrated heat-production values (i.e., the contribution of crustal radiogenic heat production to surface heat flow) in Decker et al.’s (1988) model and our results are compared in the last two columns of Table 2, and they agree very well, showing that the Curie-depth constraint can yield reasonable bulk estimates for the crustal thermal parameters. The high surface radiogenic heat production values derived with the steady-state assumption could also be a reflection of magmatic heat in the lower crust in this case.

In Salem et al. (2014), we used a similar method of Curie-depth–constrained geotherm modeling but with uniform radiogenic heat production (without derivations, albeit the same steps as in our derivation above) for deriving thermal parameters in the central Red Sea rift and the adjacent eastern and western margins. Our results from that study suggest that, in the central portion of the oceanic rift, where the Curie depth was only 6.4 km below the water depth (at the crust-mantle boundary at the window size of 100 km of that analysis), the temperature-depth profile was essentially linear. We do not have oceanic heat production measurements in this region, but the composition of the ocean crust at the mid-oceanic ridge can be assumed to be tholeiitic basalt and gabbro and the heat production somewhat uniformly low (<0.1 μW/m3, Clauser, 2011). With the uniform heat-producing layer model and surface heat flow of 160 mW/m2, surface temperature of 22 °C, and the Curie depth of 6.4 km below the ocean bottom (the parameters in Salem et al., 2014), we could only increase the average thermal conductivity in the Red Sea rift example to 1.83 W m–1 K–1. This value is at the very low end of the range for silica-poor volcanic rocks (Clauser and Huenges, 1995) and possibly represents the effect of water-filled pore space in the uppermost section and reduction of thermal conductivity by higher temperatures in the lowermost mid-ocean ridge crust. The radiogenic heat production from this calculation was 0.14 μW/m3, and the mantle heat flow was 159 mW/m2. Radiogenic heat production for the layer is essentially negligible, consistent with the oceanic crust, leading to an approximately linear geotherm. When assuming different thermal conductivity (which means changing also the mantle heat flow), the derived ratio of As/K changes.

We do not have enough high-quality, long-wavelength magnetic data available in the coastal regions of either margin of the Red Sea to use 500 km windows as we did in the United States. Using the Curie depths obtained over the Saudi Arabian margin in the Salem et al. (2014) study (∼17 km), the estimated surface heat production using the method of exponential heat production of this study (∼1.7 μW/m3) is consistent with the high end of measurements (1.46 ± 0.18 μW/m3) by Gettings and Showail (1982), although the corresponding thermal conductivities obtained through both these analyses (1.3 W m–1 K–1) is very low (assumed b = 10 km and qr = 33 mW/m2). Even increasing the Curie depth by 5 km (within the margin of error) and maintaining mantle heat flow at around 33 mW/m2 does not increase thermal conductivity over 1.65 W m–1 K–1. The interpretation of seismic refraction data 2°–3° south of this location does not suggest a great thickness of sediments (Healy et al., 1982; Mooney et al., 1985; Milkereit and Fluh,1985; Prodehl, 1985; Gettings et al., 1986), and only if the surface heat flow were higher than the 50 mW/m2 assumed by Salem et al. (2014), or if the Curie depth were grossly underestimated due to the small window used in that study, could the conductivity approach 2 W m–1 K–1. It is possible that the surface heat flow could be higher because the value of 50 mW/m2 is based on a singular observation in the region (Rolandone et al., 2013), significantly distant from the magnetic anomaly window. It is important to note that these considerations are imposed by the lack of heat-flow data and uncertainties in the Curie depth due to a small window size rather than the method of constraining temperatures at depth for deriving the geotherm. Unfortunately, there are no heat-flow and thermal parameters available near the eastern margin to compare with our results on that margin.

We have developed a new method to incorporate a temperature-depth constraint in the solution of 1-D heat-flow equation for the geotherm. In this study, we use results of spectral analysis of magnetic anomaly data to derive the depth to the bottom of a magnetic layer, which can be associated with the Curie temperature of magnetite and, hence, yields a temperature-depth constraint in the crust. Used in conjunction with the Curie depth and surface heat flow, the method yields a geotherm anchored at the Curie depth and also a ratio of radiogenic heat production to thermal conductivity (As/K). Using the estimate of the mantle heat flow as a criterion (from the intercept of the regional surface heat-flow versus heat-production relationship), one can further constrain the range of As/K and K. If the resulting value of K is reasonable, then this method also can be used to deduce As (which is an important parameter for modeling the heat flow and often not measured in the field).

With the method outlined above, we obtained the Curie-depth–constrained geotherm and estimates of surface heat production (∼8 μW/m3) in New Hampshire, which are high but comparable to observed surface heat production values over the Conway granite in the region. This result validates the method for applications where steady-state heat-flow conditions apply. We also analyzed the region across the Archean Wyoming craton to Colorado where there is a transition from a low to a high heat-flow regime. The Curie-depth map also shows a variation across these domains from the Wyoming craton (Curie depths 24–52 km, with an average of 34 km) to northern Colorado (Curie depths of 22–32 km, with an average of 26.5 km). We compared our results of the Curie-depth–constrained geotherm modeling with the finite-difference models of thermal parameters by Decker et al. (1988). Our depth-integrated crustal heat-production values compared very well with the values of Decker et al. (1988) (see Table 2).

We also placed into perspective previous Curie-depth results of Salem et al. (2014) from the Red Sea (a profile from the western to the eastern margin). In that study, which was our first application of the Curie-depth–constraint technique (albeit with the uniform radiogenic heat production assumption), we obtained reasonable parameters for the heat production of the oceanic crust (close to zero) and thermal conductivities near the low end of the spectrum of the basaltic crust in high-temperature conditions. Inland from the Saudi Arabian margin of the Red Sea, we obtained heat-production values with our method, and these are consistent with the range of a few heat-production observations in the region.

In tectonothermally younger regions, when the estimated mantle heat flow is matched in these computations, the method appears to yield higher radiogenic heat generation values and lower thermal conductivity values than observations or estimates in the region (e.g., our Colorado and the central Red Sea examples). Despite this limitation, it appears that the method will prove useful in constraining geotherms and deriving heat-production estimates where steady-state situations apply. Temperature-depth constraints are also useful in deriving geotherms in transient temperature regimes with time-varying temperature-depth formulation (e.g., Wendlandt and Morgan, 1982; Morgan, 1983; Harrison et al., 1986; Lachenbruch et al., 1994), where steady-state assumptions would be inappropriate. We intend to carry out such computations in Sierra Nevada region in the near future.

We are very grateful to Dr. Ahmed Salem for our interactions during the development of the de-fractal method. We thank the editors of Geosphere, reviewers (C. Bouligand, C. Buecker, and C. Gerbi), and Ms. Caitlyn Hartig for their comments, which improved the manuscript. Ms. Leah Newman’s cross-checking of some of our spectral depth results is greatly appreciated. This research was funded by National Science Foundation grants EAR-1246921 to DR and EAR-1246977 to ARL.


In the general procedure of Salem et al. (2014), the fractal parameter of the observed magnetic anomaly field is increased until there is consistency between depth estimates from the centroid method and the spectral peak modeling method. The procedure also requires determination of amplitude offset between the observed and modeled spectra for the peak modeling. The amplitude offset may be determined either manually or automatically by the least-squares fitting of the spectra. We have modified the procedure to make some of it interactive with some parts automatic such that in combination it results in more repeatable performance. The automation also allows us to examine a large number of solutions, in comparison to the manual approach, from which the best-fitting solutions can be selected manually and upper and lower bounds can be estimated from other equally well-fitting solutions. In the modified selection, we first fit the two amplitude (or power) spectra with the method of least squares and then readjust the fit by least-squares fitting the logarithm of amplitude (or power) over the range of wavenumbers selected on the basis of variations in the observed spectrum. The former method fits the peak region the best, and the latter fits better the overall curve of the amplitude spectrum over the selected wavenumber range. This process made the fit of the curves more comparable among different α. At the point of this writing, no mathematical misfit norm appeared to lead to the same selection as we might visually make because misfit functions minimize the difference between the observed and modeled curves, and as a result, the curves cross over in different segments (see Fig. A1a and the description below). On the other hand, for selecting correct depth to the bottom, we want to ensure that particular segments of the curves fit well without crossing over. We also tried several mathematical misfit norms (L2, L1, and L-infinity), but none could avoid the cross-over problem.

Based on model studies, repeatable and accurate estimates of the depth to the bottom are achieved routinely at a fractal parameter of the field (α) in the neighborhood of the one that leads to a minimum in the misfit. In Figure A1a, the spectra at the minimum root-mean-square (RMS) misfit between the observed and modeled spectral curves (α = 2.1) from a model study are shown (where fractal parameter of magnetization β is 3; and so the fractal parameter of the field (α) should be close to 2; Maus and Dimri, 1994, for near-surface magnetic field data); one can observe cross-overs in the fit in the wavenumber range from 0 to 0.1—the key part of the curve for the depth to the bottom. Figure A1b shows the RMS misfit curve as a function of α. In the example shown, the fractal parameter at the RMS tightly hugs the left face of the peak; however, routinely, the best solutions are within ±0.2 of αminimum misfit. Another change in the software that improved the repeatability of our results is, instead of picking the centroid slope graphically (which led to different number of wavenumbers selected in different runs), determining the steepest centroid slope using a fixed number of consecutive wavenumber points (e.g., 2 or 3 points depending on the waviness of the curve right of the peak). Making sure that in the accepted solutions the slope of the de-fractaled centroid spectrum in the selected range is straight (and also not turned over due to de-fractalization) is also important.

Figure A2 shows the result of the model study for four models with different depths to the bottom of magnetization and the average and standard deviation of the results obtained with nine offset windows. The match is very good for a 500 km × 500 km window but results in underestimation of basal depths obtained from smaller window sizes (not shown to avoid clutter). In the figure, we also make comparisons with the approach of Li et al. (2013). They used automatically determined least-square slopes with preselected wavenumber ranges with the fractal parameter–corrected centroid method (Bansal et al., 2011) to compute depths to top and bottom with small windows. Based on the results shown in Figure A2 and our experience with smaller windows, neither preselected wavenumber ranges nor small windows give consistent results.

During this study, we also developed a 3-D de-fractal approach where the fractal parameter of 3-D magnetization was iteratively estimated for each fractal parameter of the field being computed; however, the fractal parameters of the field and the magnetization may be dependent on each other, and consequently the misfit curve between observed and modeled spectra could not be used effectively in the selection of the solution. In model studies, the 2-D de-fractal method performed well on 3-D fractal magnetization distributions, and we attribute this primarily to the larger influence of layers near the top on the magnetic field than the layers at greater depth. Thus for all practical purposes, the magnetic field from a 3-D fractal crustal magnetization behaves as if it was generated by a 2-D fractal magnetization distribution.


The original papers introducing the linear heat-flow (Q)–heat-production (A) relation only claimed its application in major “granite” plutons and used these data to define “heat-flow provinces” (e.g., Birch et al., 1968; Roy et al., 1968). Lachenbruch (1968, 1970) introduced the proposal that heat production decreased exponentially with depth in these terrains in order to explain how the Q-A relation could survive differential erosion. He discussed how the heat-producing elements could be integrated in a vertical profile in a granitic pluton during melting. Lachenbruch also proposed a linear decrease in A with depth. This relation does not satisfy the condition of differential erosion, but it was within the error limits of the small data set of A versus depth available to Lachenbruch (1968). Other workers have extended the relation to mixed plutonic and/or metamorphic terranes or just metamorphic terranes. There is a general (statistical) Q-A linear relation in these terrains but of much poorer quality than the relation when applied to provinces with major “granite” plutons. There is no geochemical reason why the relation should be applicable to metamorphic terranes because the surface rocks in these terranes have no mechanism to be a representative “sample” of the upper crust.

Jaupart and Maraschal (2007) have claimed that there are no data to support the exponential model. Of the crustal sections that they cite, only two are granitic sections—the Vredefort structure and the Sierra Nevada (we exclude sections including metamorphic crust because they were not in the original definition of Q-A provinces and are not considered in this paper). If the data published for these sections are examined, the results are clearly indicative of weathered samples. Unweathered samples should have Th/U ratios of ∼3.8–3.9 (Taylor and McLennan, 1985; Rocholl and Jochum, 1993) and K/U ratios of ∼1–1.3 ± 0.3 × 104 (Taylor and McLennan, 1985; Rocholl and Jochum, 1993; Arevalo et al., 2009). Data published for the Vredefort structure (Nicolaysen et al., 1981) have Th/U and K/U ratios ranging from 4.7 to 24.4 (mean 15.5 ± 9.6, n = 14) and 1.0–14.9 × 104 (mean 6.3 ± 5.0 × 104, n = 14), respectively. Data published for the Sierra Nevada plutons (Brady et al., 2006) have Th/U and K/U ratios ranging from 0.1 to 91.7 (mean 6.6 ± 6.5, n = 57) and 0.3–12.5 × 104 (mean 2.1 ± 2.5 × 104, n = 57), respectively. These ranges of ratios are much greater and their means are far different from the bulk earth ratios expected for unweathered samples than ratio from borehole samples, e.g., from Decker et al. (1988): Th/U mean 4.4 ± 2.6 (n = 38) and K/U mean (1.0 ± 0.6) × 104 (n = 38).

Another indication that the samples from which the results presented by Brady et al. (2006) were probably weathered is that the mean heat production from these samples is 1.5 ± 1.2 μW/m3 (n = 57), whereas the mean heat production from borehole samples reported by Saltus and Lachenbruch (1991) is 2.1 ± 1.1 μW/m3 (n = 28). Although these means are not significantly different, the lower mean and larger standard deviation with the surface samples, although the sample size is larger, is what would be expected with weathered samples.

The analysis of the data presented above that have been used to “prove” that the exponential distribution is not valid demonstrates that samples used for the “proof” were probably weathered. Therefore the results from the surface samples from the Vredefort structure and the Sierra Nevada cannot be used to test the validity of the exponential distribution of heat production through the crust in granitic plutons.

The heat-flow data for the three test areas used in this paper are primarily from silicic plutons—the U.S. New England plutons, granitic plutons in northern Colorado and southern Wyoming in the United States, and plutons in the Saudi Arabian Shield adjacent to the Red Sea. Only New England is a well-defined heat-flow province because the other two areas are transitional heat-flow zones. Our analysis, as any model calculating crust geotherms, required an assumption of the distribution heat production with depth. As we have demonstrated above, there are no data invalidating the exponential distribution in major silicic plutons, and this distribution satisfies the documented observation of the linear Q-A relation with differential erosion in granitic pluton heat-flow provinces with high-quality heat-flow data. We do not claim that the distribution is valid in metamorphic terrains. One goal of the present study is to find and test temperature tie points deep in the crust that may be used to constrain crustal thermal parameters.