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

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.


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 peer-

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 0.75 cm ± ) below the surface with cables leading to an above-ground data logger and master computer. The manufacturer-specified operating range of the thermistors is ( ) 35, 50 C − +  with measurement error below 0.4 C ±  . 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 47,240 N = 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 ( ) 10 x t , 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  weather. The black trace ( ) 240 x t from the deepest probe, which resembles a smooth sinusoidal function, manifests only seasonal variations in temperature.
The records 10 x and 240 x 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 It is worth noting that the two kinds of averages are structurally different. In effect, the 24-hour average 24 h A 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 365 d

MA
replaces each point by a sum of 365 points, thereby shortening the input series by a length of 365.

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 240 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 Hartford 0.28 0.0032 C y 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] Furthermore, a separate regression analysis (to be published elsewhere [12]) of above-ground temperature measurements ) 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 Hartford τ is anomalous among cities, an analysis was made of the time series (1900-2012) 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.

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 10 x and 240 x to confirm that the moving-average series 365 d 240

MA
x is unaffected by daily fluctuations.
To eliminate a static term leading to a zero-frequency spike, time series d x were transformed to series d y of zero mean From the Fourier amplitudes of d y in which δ is the Kronecker delta function, respectively follow the magnitudes, phases, and power spectral amplitudes For discrete data collected at a sampling rate 1 t ∆ , the cut-off frequency is  in accordance with Shannon's sampling theorem [14]. The harmonic T j corresponding to a particular period T (in units of ∆t ) in a time series of length N is , and the time-variation of that harmonic component is then . 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 10 S , 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 0 16, 000 j ≤ ≤ . The only other statistically significant peak occurs at 4 j = (not shown), corresponding to the annual period y T 2 . A double-log plot of 240 S 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 series 3 ( ) for the diurnal harmonic spectrum ( ) 2 The peak at y T theoretically occurs at 3.74 j = . However, harmonic numbers must be integers. 3 The 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 2 n   for the quasi-continuous spectrum 240 S (Figure 5). Exponent (14) is nearly the same for 10 S upon removal of the diurnal harmonic content.
Further elucidation of the spectral content of 240 S is obtained by examining the power spectra of the 24-hour averaged detrended series 24  MA y as shown in Figure 7. The oscillatory structure seen in the spectrum of 24 h 240 A y does not signify multiple periodicities, but in fact is due entirely to seasonal variability at annual period y T . The sequence of cusps occurs at characteristic of 1D Brownian diffusion. The periodicity at y T 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 ly resembles a log hyperbolic density, which is widely found to characterize particle mass and size distributions  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 0 δ > is a scaling parameter. The upper panel of Figure 10 shows a histogram of 240 log log S , which reveals the hyperbolic form of the exponent in Equation (16)  log MA S in the lower panel of Figure 9 strongly resembles an exponential pdf whose general form is Correspondingly, a plot of ( ) 365 d 240 log log MA S in the lower panel of Figure 10 reveals the linear form of the exponent in Equation (17) with visually fit parameter 0.53 λ = .
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.

Dynamics of Solar Energy Diffusion
A 1D diffusion model for temperature ( ) , T x t , based on the heat-flux equation [19] ( ) ( ) depends on ground thermal conductivity T κ , mass density ρ , and specific heat capacity c . 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 ( ) which reduces to ( ) Equation (22) permits estimation of the diffusivity D by two independent methods: 1) phase shift between corresponding maxima at times ( ) 1 2 , t t of two temperature records of known depths ( ) and 2) amplitude attenuation of the peak temperature max T recorded at two known depths Applied to the plots of Figure 4, Equations (23) and (24) led respectively to a mean diffusivity  (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 365 d 240 MA x .

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 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 FBÑ 3 1 β > > is more predictable than white noise (WN) for which WN~0 β , such as observed in decay of radioactive nuclei [24]. For FBN 2 β > , the future trend follows the past trend, whereas it reverses for FBN 2 β < . Brownian motion with FBN 2 β = 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.