Journal of Modern Physics
Vol.5 No.9(2014), Article ID:47101,12 pages DOI:10.4236/jmp.2014.59085

Statistical Analysis of Subsurface Diffusion of Solar Energy with Implications for Urban Heat Stress

M. P. Silverman

Department of Physics, Trinity College, Hartford, USA


Copyright © 2014 by author and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

Received 8 April 2014; revised 2 May 2014; accepted 24 May 2014


Analysis of hourly underground temperature measurements at a medium-size (by population) US city as a function of depth and extending over 5+ years revealed a positive trend exceeding the rate of regional and global warming by an order of magnitude. Measurements at depths greater than ~2 m are unaffected by daily fluctuations and sense only seasonal variability. A comparable trend also emerged from the surface temperature record of the largest US city (New York). Power spectral analysis of deep and shallow subsurface temperature records showed respectively two kinds of power-law behavior: 1) a quasi-continuum of power amplitudes indicative of Brownian noise, superposed (in the shallow record) by 2) a discrete spectrum of diurnal harmonics attributable to the unequal heat flux between daylight and darkness. Spectral amplitudes of the deepest temperature time series (2.4 m) conformed to a log-hyperbolic distribution. Upon removal of seasonal variability from the temperature record, the resulting spectral amplitudes followed a log-exponential distribution. Dynamical analysis showed that relative amplitudes and phases of temperature records at different depths were in excellent accord with a 1-dimensional heat diffusion model.

Keywords:Time Series Analysis, Heat Conduction, Thermal Diffusion, Power Laws, Climate Change, Heat-Island Effect

1. Introduction

Despite a long period of contentious debate, there is now little dissention among scientists who study global climate change that the change is real and likely dominated by human activity [1] . A recent survey of the peerreviewed scientific literature concluded that the number of papers rejecting the consensus on anthropogenic global warming was a “vanishingly small proportion of the published research” [2] .

Evidence for increasing mean global temperature has come from four major sources: 1) the US NASA Goddard Institute for Space Studies (NASA GISS), 2) the US National Oceanic and Atmospheric Administration (NOAA), 3) the UK Climatic Research Unit (CRU), and most recently 4) the US Berkeley Earth Surface Temperature group (BEST) whose acronym reflects what the group believe to be the largest sampling and smallest uncertainties of the major investigations to date [3] .

Reports in both the scientific literature and popular news media have focused nearly exclusively on global aspects of global warming [4] —i.e. on planet-wide catastrophes such as melting of the Antarctic ice sheet, disruption of the thermohaline circulation in the Atlantic, increased frequency and intensity of extreme storm events, and regional perturbations of ecosystems leading to large-scale extinctions. Although global catastrophes are considered within the realm of future possibility, the most serious consequences of climate change to occur soonest will likely be local, not global, and have direct impact on human health [5] , in particular on mortality associated with exposure to ambient temperature [6] .

According to the US Centers for Disease Control and Prevention, heat waves are the most deadly weatherrelated exposure in the US, and account for more deaths annually than hurricanes, tornadoes, floods, and earthquakes combined [7] . Comparable exigencies have also been reported recently for Europe [8] . In measurements of surface land temperatures for the purpose of estimating the background rate of global temperature rise, researchers have systematically excluded measurements in urban areas to avoid the so-called heat-island effect, i.e. the phenomenon that urban areas are ordinarily hotter than rural areas. The effect reflects such urban conditions as fewer trees to provide shade or capture moisture, more asphalt and cement to absorb heat, tall buildings that block heat radiation into space, and other reasons. Thus, the majority of surface temperature studies disregard the principal locations—cities—where impact of temperature change on human living conditions is likely to be the most profound.

Previous attempts have been made to compare urban and rural temperature trends [9] by examining data from surface stations and making statistical adjustments for differences in measurement conditions. Surface temperature measurements, however, are strongly impacted by the variability of daily weather, which increases the uncertainty of low-amplitude trends.

This paper reports on subterranean temperature measurements made over a 5+ year period as a function of time and depth in a medium-size US city (Hartford, Connecticut:,) and compares the results with surface temperature measurements made in New York City (NYC:,), which is currently the largest US city by population. The two temperature trends are of comparable magnitude, and both are alarmingly higher than the mean regional temperature, which is close to the rate of global temperature rise.

The stochastic process giving rise to the subterranean temperature series can be understood in its essential details as a one-dimensional (1D) heat diffusion process. Knowledge of the diffusion constant, or diffusivity, which was determined from the data in two independent ways, together with shallow subsurface boundary conditions, permits accurate prediction of the temperature series recorded at greatest depth, thereby substantiating the interpretation of the inferred temperature trend. Power spectral analysis of the time series as a function of depth revealed two different power-law variations, which shed further light on underlying stochastic processes.

2. Experimental Procedure and Observed Time Series

Thermal diffusion of solar energy as a function of time and depth was measured by a series of 6 thermistor sensors positioned in a vertical conduit at depths of 10, 20, 40, 80, 160, and 240 cm (uncertainty) below the surface with cables leading to an above-ground data logger and master computer. The manufacturer-specified operating range of the thermistors is with measurement error below. The six temperatures were recorded every hour on the hour starting at noon on 7 June 2007. Data analyzed in this paper cover a period of approximately 5.4 years (i.e. through 2012) and comprise observations.

A panoramic sample of the resulting time series collected during the first two and a half years and colorcoded for depth is shown in Figure 1. The red trace, designated, from the probe nearest the surface is most sensitive to weather-induced temperature fluctuations, which is the predominant source of noise. With increasing depth d, the corresponding temperature time series show diminishing sensitivity to daily

Figure 1. Panoramic plots (truncated to 2.5 y) of hourly temperatures during the period 2007-2012 by sensors at depths (in cm) of 10 (red), 20 (blue), 40 (green), 80 (magenta), 160 (gold), and 240 (black).

weather. The black trace from the deepest probe, which resembles a smooth sinusoidal function, manifests only seasonal variations in temperature.

The records and are shown in greater detail in Figure 2. Superposed over the hourly records (red), covering a period of 2000 days, are the 24-hour averaged records (black) defined by


for d = 10 and 240 cm, respectively, in which is the floor function1. At a depth daily fluctuations are absent, and there is little difference between and. Solid black traces close to baseline record the 365-day moving average time series defined by

. (2)

It is worth noting that the two kinds of averages are structurally different. In effect, the 24-hour average replaces each suite of 24 points by 1 point, thereby shortening the input series by a factor of 24. In contrast, the 365-day moving average replaces each point by a sum of 365 points, thereby shortening the input series by a length of 365.

3. Abrupt Increase in Rate of Temperature Rise

In combination, transformations (1) and (2) remove nearly all daily and annual variations from the time series. What remains, as shown at larger scale in the plot of in Figure 3, are weakly varying residuals with long-term trend. The maximum-likelihood (ML) line of regression (heavy dash) yields a rate of temperature increase and standard error of

. (3)

To put this number in perspective, one can compare it to (a) the mean temperature increase for the US Northeast reported by the Union of Concerned Scientists [10]


and (b) the global mean annual temperature rise over the past 30 years reported by the US National Research Council [11]

. (5)

Furthermore, a separate regression analysis (to be published elsewhere [12] ) of above-ground temperature measurements (1960-2012) collected at about 6 km outside the city yielded very close to the same rate as (4). Clearly, the much higher rate of temperature rise (3) is a recent phenomenon.

To ascertain whether is anomalous among cities, an analysis was made of the time series (1900-2012)

Figure 2. Comparison of hourly (red), 24-hour time-averaged (thin black), 365-day moving average (heavy black), and mean (thin dashed black) temperature records at depths of 10 cm (upper panel) and 240 cm (lower panel).

Figure 3. Maximum-likelihood (ML) line of regression (heavy dashed black) to the 365-day moving average temperature series (solid black) recorded at a depth of 240 cm.

of surface temperatures of NYC [13] collected at a station in Central Park, Manhattan. The ML slope and standard error of the line of regression to the time series for three different spans of time are


The medium-span temperature trend is consistent with the global mean rate (5), and the short-span is consistent with the Hartford rate (3), but larger because NYC has a higher population and population density.

4. Power Spectra and Autocorrelation

Power spectra of the subterranean temperature time series provide detailed information about the stochastic processes by which solar energy propagates through the ground. Of particular interest is the comparison of the spectra of time series and to confirm that the moving-average series is unaffected by daily fluctuations.

To eliminate a static term leading to a zero-frequency spike, time series were transformed to series of zero mean

. (7)

From the Fourier amplitudes of

, (8)

in which and is the Kronecker delta function, respectively follow the magnitudesphases, and power spectral amplitudes



. (11)

For discrete data collected at a sampling rate, the cut-off frequency is, beyond which data are aliased, i.e. folded into a lower frequency range. The fundamental is the inverse of the duration of the time series. Each frequency is labeled by a harmonic number j. The largest harmonic is in accordance with Shannon’s sampling theorem [14] .

The harmonic corresponding to a particular period T (in units of) in a time series of length N is, and the time-variation of that harmonic component is then

. (12)

Figure 4 shows the six waves of harmonic corresponding to annual period

for a series truncated at. These waves reveal the phase shifts and amplitude differences in the large-scale variation of the time series in Figure 1 without accompanying noise.

Of particular interest is the comparison of power spectra obtained from the 10 cm and 240 cm temperature sensors. A double-log plot of, shown in the upper panel of Figure 5, reveals a striking pattern of discrete peaks superposed over a quasi-continuum of harmonics over the range. The only other statistically significant peak occurs at (not shown), corresponding to the annual period 2. A double-log plot of in the lower panel shows the quasi-continuum of harmonics, but no evidence of discrete peaks. The sequence of discrete peaks, shown separately as a double-log plot in Figure 6, corresponds precisely to the harmonic series3 of the diurnal period. All three double-log plots reveal a power-law dependence on frequency over the observed range. Regression analysis leads to ML exponents


for the diurnal harmonic spectrum (Figure 6) and


Figure 4. Components at the fundamental frequency in the Fourier analysis of temperature time series (a); (b); (c); (d); (e). Phase shifts and amplitudes relative to allow determination of the diffusivity D and agree with results of a one-dimensional (1D) diffusion model.

Figure 5. Double-log plots of power spectra (red points) (top panel) and (bottom panel) superposed by ML lines of regression (solid black). Both plots manifest a quasi-continuum of Brownian noise. also reveals a discrete series of diurnal harmonics.

for the quasi-continuous spectrum (Figure 5). Exponent (14) is nearly the same for upon removal of the diurnal harmonic content.

Further elucidation of the spectral content of is obtained by examining the power spectra of the 24-hour averaged detrended series and the 365-day moving average of the detrended series as shown in Figure 7. The oscillatory structure seen in the spectrum of does not signify multiple periodicities, but in fact is due entirely to seasonal variability at annual period. The sequence of cusps occurs at

Figure 6. Double-log plot of (open circles) for the diurnal harmonics, superposed by ML line of regression (dashed).

Figure 7. Double-log plots of (black dots) and (red dots) with ML line of regression (dashed black).

frequencies close to where ideally (i.e. in absence of noise) the Fourier transform of a pure sinusoid of period vanishes, and the log approaches. This interpretation is confirmed by the autocorrelation function (red trace) of in Figure 8, which is fit exactly by the theoretical autocorrelation function

(dashed black trace). Upon removal of seasonal variability in the series, the sequences of cusps vanishes from the power spectrum in Figure 7 (red dots), which is then fit very closely by a ML line of regression leading to power law exponent


characteristic of 1D Brownian diffusion. The periodicity at largely vanishes from the autocorrelation in Figure 8 (solid black trace), which then likewise (in accord with the Wiener-Khintchine theorem) manifests a long-range decay characteristic of a Brownian process.

Consider next the statistical distribution of the power spectral amplitudes—or a function of such amplitudes—which is useful in revealing empirically the probability density function (pdf) of a stochastic process. The objective of such an examination is to arrive at a recognizable form of pdf. In Figure 9 are shown histograms of (upper panel) and (lower panel). In appearance the histogram of strongly resembles a log hyperbolic density, which is widely found to characterize particle mass and size distributions

Figure 8. Autocorrelation functions of (solid red) and (solid black); theoretical autocorrelation (dashed black).

Figure 9. Upper Panel: histogram of conforms to a hyperbolic distribution. Lower Panel: histogram of conforms to an exponential distribution. Elements of both histograms are partitioned among 100 classes covering a harmonic range.

in fragmentation processes [15] [16] , as well as stochastic processes involving diffusion of matter, energy or the movement of stock prices [17] . The pdf takes the general form [18]


in which and are the slopes of the asymptotes of the hyperbolic exponent, is a location parameter, and is a scaling parameter. The upper panel of Figure 10 shows a histogram of, which reveals the hyperbolic form of the exponent in Equation (16) with visually fit asymptotic parameters,

. The histogram of in the lower panel of Figure 9 strongly resembles an exponential pdf whose general form is

. (17)

Correspondingly, a plot of in the lower panel of Figure 10 reveals the linear form of the exponent in Equation (17) with visually fit parameter.

Although further investigation of the dynamical origin of these distributions is in progress, for the present purposes it suffices that spectral analysis and autocorrelation have revealed all statistically significant periodicities in the temperature time series.

5. Dynamics of Solar Energy Diffusion

A 1D diffusion model for temperature, based on the heat-flux equation [19]

, (18)

allows adequate prediction of time series from the non-stochastic component of. The diffusivity


depends on ground thermal conductivity, mass density, and specific heat capacity. The positive direction of flow (x axis) is into the ground with origin at the surface.

Solving Equation (18) for experimental conditions of this paper leads to


in terms of a theoretically known or empirically determined function of angular frequency. When initial amplitudes and phases are known for a discrete set of fundamental frequenciessubstitution of into (20) leads to the general form

, (21)

which reduces to


in the case of a detrended time series of single fundamental.

Equation (22) permits estimation of the diffusivity D by two independent methods: 1) phase shift between corresponding maxima at times of two temperature records of known depths

, (23)

and 2) amplitude attenuation of the peak temperature recorded at two known depths

. (24)

Applied to the plots of Figure 4, Equations (23) and (24) led respectively to a mean diffusivity

Figure 10. Upper Panel: histogram of with visual fit to hyperbolic asymptotes yielding slope parameters,. Lower Panel: histogram of with visually fit line of regression yielding slope parameter.

. (25)

Given the value in (25) and the amplitudes and phases


from Fourier analysis of the mean-adjusted 3-year temperature record, including contributions at periods

, , and (first diurnal harmonic), one can predict from Equation (21), as shown in Figure 11.

The empirical (solid blue) and predicted (dashed black) time series agree very closely in amplitude and phase. The evidence supports the hypothesis that the observed subterranean temperature time series are the result of predominantly one-dimensional solar heat diffusion and that only this solar energy diffusion contributes to the trendline of.

6. Conclusions

In contrast to numerous previous studies of global temperature trends by measurements at rural surface stations, this paper reports the measurement and statistical analysis of subterranean temperature within a medium-size city to determine a specifically urban temperature trend. Subterranean measurements at depths of 2 m or more are virtually unaffected by daily temperature fluctuations and are sensitive only to seasonal variations. The mean rate of temperature rise in a medium-size North American city (Hartford CT) between 2007 and 2012, after

Figure 11. Comparison of detrended temperature series (red) and truncated Fourier series (solid black) at fundamental periods and. Comparison of detrended temperature series (blue) and solution (dashed black) to the diffusion Equation (18) with boundary conditions based on.

transformations eliminating diurnal and seasonal variability, was found to be more than ten times the background rate of global warming. This finding is consistent with surface temperatures recorded in central New York City, the second largest city in North America by population. There is no reason to think that these findings may be atypical for cities in industrialized countries. Given recent occurrences of extreme heat waves in the US and Europe, a precautionary conclusion, substantiated by climate modeling, is that, left unmitigated, the high rate of urban temperature rise relative to global warming will lead to more frequent and intense urban heat stress [20] [21] .

Spectral analysis of the subterranean temperature time series at different depths revealed two distinct power laws suggestive of self-affine fractal behavior over the observed frequency range [22] . The exponent ~2 of the quasi-continuum spectrum is essentially independent of depth and corresponds to fractal Brownian noise (FBN) [23] . In general, FBN with is more predictable than white noise (WN) for which, such as observed in decay of radioactive nuclei [24] . For, the future trend follows the past trend, whereas it reverses for. Brownian motion with lies at the threshold between persistence and anti-persistence because increments are uncorrelated.

The diurnal harmonic spectrum in the shallow subsurface temperature record can be shown [12] to arise from a non-sinusoidal heat flux due to unequal durations of daylight and darkness and the difference in rates of daytime heat absorption and nocturnal cooling. The large power-law exponent ~5.8 characterizes a relatively smooth function with behavior predictable by linear extrapolation. The development of an empirical model that accurately reproduces the diurnal power spectrum for purposes of understanding surface energy fluxes is in progress and will be reported when completed.


The author thanks Drs Jon Gourley and Christoph Geiss for several helpful conversations and the staff who maintain the Trinity College Weather Station and Temperature Well. Students who participated at short intervals in the project include Sarthak Khanal, Matthew Cohen, and Robert Tella.


  1. Tollefson, J. (2011) Nature News.
  2. Cook, J., et al. (2013) Environmental Research Letters, 8, Article ID: 024024.
  3. Rohde, R., et al. (2013) Geoinformatics & Geostatics: An Overview, 1, 1.
  4. IPCC (2013) Climate Change 2013: The Physical Science Basis.
  5. McMichael, A.J., Woodruff, R.E. and Hales, S. (2006) The Lancet, 367, 859-869.
  6. Hajat, S., Vardoulakis, S., Heaviside, C. and Eggen, B. (2014) Journal of Epidemiology & Community Health, 68, 641-648.
  7. Centers for Disease Control and Prevention: Climate and Health Program.
  8. Robine, J.-M., et al. (2008) Comptes Rendus Biologies, 331, 171-178.
  9. Stone Jr., B. (2007) International Journal of Climatology, 27, 1801-1807.
  10. Union of Concerned Scientists (2006) Climate Change in the US Northeast. UCS Publications, Cambridge, MA, 10.
  11. Committee on America’s Climate Choices (2011) America’s Climate Choices. National Academy Press, Washington DC, 15.
  12. Silverman, M.P. (2014) A Certain Uncertainty: Nature’s Random Ways. Cambridge University Press (to be published), Cambridge.
  13. NYC Temperature Data.
  14. Shannon, C. (1949) Proceedings of the IRE, 37, 10-21.
  15. Silverman, M.P., Strange, W., Bower, J. and Ikejimba, L. (2012) Physica Scripta, 85, Article ID: 065403.
  16. Fieller, N.R.J., Flenley, E.C. and Olbricht, W. (1992) Applied Statistics, 41, 127-146.
  17. Bibby, B.M. and Sorensen, M. (1997) Finance and Stochastics, 1, 25-41.
  18. Barndorff-Nielsen, O. (1977) Proceedings of the Royal Society of London. Series A, 353, 401-419.
  19. Evett, S.R., et al. (2012) Advances in Water Resources, 50, 41-54.
  20. Wu, J., et al. (2014) Environmental Health Perspectives, 122, 10-16.
  21. Meehl, G. and Tebaldi, C. (2004) Science, 305, 994-997.
  22. Turcotte, D.L. (1997) Fractals and Chaos in Geology and Geophysics. Cambridge University Press, Cambridge, 148-149.
  23. Hergarten, S. (2002) Self-Organized Criticality in Earth Systems. Springer, Heidelberg, 48-66.
  24. Silverman, M.P. and Strange, W. (2009) Europhysics Letters, 87, Article ID: 32001.


1 is the greatest integer less than or equal to n.

2The peak at theoretically occurs at. However, harmonic numbers must be integers.

3The word “harmonic” has two meanings, both of which are used in this paper. In physics, it is a multiple of a fundamental frequency. In mathematics, it is also a series of terms proportional to .