A method to evaluate the impact of urbanization on ground temperature evolution at a regional scale

The replacement of natural soil and vegetation by arti ﬁ cial surfaces increases temperatures of the surrounding air and subsurface throughout the year, because of indirect solar heating of urban structures, building heat losses and land use change. This phenomenon is called Urban Heat Island and it can be better perceived during the nighttime, when the city releases the heat accumulated during the day. During the daytime, due to relatively small amounts of solar radiation received by the urban surface, especially in high-density cities in arid and semi-arid climates, the Urban Cool Island can be identi ﬁ ed as well. The present work illustrates a mixed probabilistic-deterministic method to estimate ground temperature at shallow depths, starting from information on geology, hydrogeology, climate, but also urban presence, through correlations with global land cover and population density. A dedicated mapping on a regular grid has been produced. Results have been compared with ground and aquifer temperature available in literature, for some representative cities of the Italian Peninsula and the Alpine Zone. Preliminary validations are encouraging and can be taken as a starting point for more comprehensive mapping of ground temperature evolution at a regional scale.


Introduction
The subsurface temperature gradient, or geothermal gradient, depends mainly on endogenous geothermal heat fl ow (Cermak and Rybach, 1979). However, in shallow layers, its contribution is limited, and subsoil temperature follows climate seasonality, which is dampened because of ground thermal insulation that varies with geological and hydrogeological conditions (Kusuda and Achenbach, 1965). A shallow layer ground temperature should also take into account urban ground warming (Ferguson G. and Woodbury A.D., 2004). The replacement of natural soil and vegetation by artificial surfaces increases temperatures of the surrounding air and subsurface throughout the year because of indirect solar heating of urban structures, building heat losses and land use change (Bornstein, 1968). At a district or city level, this phenomenon is called Urban Heat Island -UHI (Landsberg, 1981), and it can be better perceived during the nighttime, when the city releases the heat accumulated during the day.
On the other hand, during the daytime, due to relatively small amounts of solar radiation received by urban surface, especially in high-density cities, the Urban Cool Island -UCI can be identifi ed as well. This effect practically disappears in low-density cities and in continental, oceanic and cold Mediterranean climates, while it cannot be neglected in warm Mediterranean, semi-arid and arid climates (Alonso M.S. et al., 2003;Rasul et al., 2015), in which surrounding environments tend to be warmer than in the urban area (Sauceda, 2014).
These effects are both transferred to the underground, affecting ground temperature in the shallow layers (Oke T.R., 1982). It follows that, generally, cities in colder climates evidence a shallow ground temperature higher than what is naturally expected (even up to 5°C). On the contrary, cities in warmer climates evidence a shallow ground temperature similar to or even lower than what is naturally expected, experiencing a balance between nightly UHI and daily UCI or being UCI prominent (Serageldin A. A. et al., 2015). The deterministic heat transfer models applied to the determination of UHI and UCI and their effects to the ground are limited by unavoidable approximations (Kato & Yamaguchi, 2005). Remote sensing techniques can be used to map surface temperature (Bitelli et al., 2015); some studies were carried out to map subsoil temperature through satellite derived data, but the level of investigation is currently limited to very shallow depths (Wenfeng et al., 2014). Due to the limits of physical methods, a statistical and geostatistical approach is needed to expand the study of ground temperature to Geostatistics provides quantitative descriptions of natural variables distributed in space or in time and space (Chiles and Delfi ner, 2012), hence, it can be useful to predict ground temperature space and time variability in an urban area, through analysis of the available information. The improved quantifi cation of ground temperature by geostatistical techniques, with respect to the application of only deterministic models, can give an advantage to many engineering topics, such as insulation design, heat recovery and shallow geothermal systems.
This work presents the use of statistical and geostatistical tools for evaluating the impact of urbanization on ground temperature evolution at a regional scale. An urban zone, in the Zurich city area, with measurements of borehole temperature profi les, was selected. Temperature profi les were used to calculate the time-varying vertical temperature profi le at the point, knowing climate, geological and geothermal information. Then, a geostatistical estimation was performed, and the spatial variability of temperature in depth was estimated through the Universal Kriging (UK) method (Chiles and Delfi ner, 2012). A spatial correlation among the estimated temperature, depth of bedrock (Pelletier et al., 2016), population density (PDN) (CIESIN, 2016) and global land cover (GLC) (Arino et al., 2010) was calculated in a regular grid for the area of Zurich.
To create a homogenous group of data, a grid of 1 km × 1 km × 1m was identifi ed, and the space-time ground temperature profi le was calculated for each node of the grid -approximation of the combined effect of UHI and UCI on the area. The calculations have been implemented in different locations in Europe and validated comparing average ground temperature at the fi xed depths with measurements and literature data. The fi nal target was the estimation of ground temperature, at each point of the grid and along vertical depth, inclusive of the subsurface urban heat island (SUHI) effect.
The results of the work can be a starting point for the use of geostatistical analysis in estimating ground temperature in complex urban areas, and can be useful for ground temperature assessment on urban settlements.

Theory of ground temperature evolution
In shallow layers, the fourth-dimensional (x, y, z, t) temperature distribution is the result of many interrelated factors (spatially distributed on a different scale), which are: • seasonal climatic variations; • soil thermal properties; • hydrogeological properties; • urbanization and related UHI; • snow cover; • geothermal heat fl ow. The basic equation for assessing the variation of vertical temperature distribution (T g ) over time is generically a function of the ambient temperature wave, the thermal properties of the ground layers and the geothermal gradient. Since all the variables entering the function are regionalized, the target variable is four dimensional, varying in space and time. Equation 1 summarizes the well-known temperature T distribution in the subsoil (Baggs, 1985): Where T m is the annual surface average temperature at location xy (°C), A is the wave amplitude at location xy(°C), the year is the wave period, z is depth (m), t is time (days), t T0 is the time at minimum temperature (days), a is the equivalent thermal diffusivity on the depth of investigation (m 2 /days) and ∇T is the geothermal gradient at location xy (°C/m), depending on geothermal heat fl ow h f (W/m 2 ) and equivalent thermal conductivity on the depth of investigation l (W/(m•K)) (see Equation 2): (2) When information on these parameters are known or estimated, it is possible to create a hypothetical 4D map of ground temperature T g in an interest area.

The impact of urban heat island
The effect of urban areas on the ground temperature gradient was studied using the measurements of Bayer in the Zurich city area (Bayer et al., 2016). Four datasets of ground temperature values in four locations near Zurich were used -starting at 20 meters below ground level in order to avoid any seasonal infl uence due to surface atmospheric temperature. The datasets of population density (PDN) and global land cover (GLC) were used to obtain a spatial correlation between the reported results and the specifi c localization of the datasets, through probabilistic techniques. A deterministic function was used to compare the temperature curve in a given location, infl uenced by the SUHI and without considering any SUHI effect. Two variables were the objective of the analysis: • The temperature difference ΔT (°C) caused by the SUHI according to PDN and GLC of the location; • The modifi ed temperature gradient of the underground down to the investigation depth (50 m).
Borehole temperature measurements were available in four locations in the area of Zurich, as well as in the suburbs around the lake of Zurich. The four measurements were named according to their geographical positions: Southeast, South, East and Central (Bayer et al., 2016). As shown in the map (see Figure 1), the East (3), South (1) and Southeast (4) boreholes are in peripheral areas while the Central (2) borehole is inside the urban area.
Data on the geology of the area, ambient temperature and amplitude, geothermal heat fl ow, thermal conductivity of sediments and bedrock, thermal diffusivity of sediments and bedrock and thickness of sediments were known for the four selected points with a precision of 1 km ×1 km grid (Tinti et al., 2018). The ground temperature profi les were calculated using Equation 1, without taking into consideration any SUHI effect. To calculate the effect of the SUHI, probabilistic analysis was performed and results were validated by comparing the calculated ground temperature with the measured data in the location of the four boreholes.

Probabilistic analysis
To perform the probabilistic analysis, the calculated temperature (T g ) in depth from a deterministic function was used for variography. Moreover, the sample variograms of temperature T g in the vertical direction was calculated for the four borehole temperature measurements. Theoretical variogram models were fi tted on sample variograms. Table 1 shows, on the left column, the sample variograms and models from the measured temperature data, reported by Bayer (Bayer et al., 2016) -infl uenced by the SUHI-; on the right column, the table shows the sample variograms and models from the average annual temperature data calculated by applying Equation 1 -without the infl uence of the SUHI.
The fi tted theoretical variogram models have two structures: the Spherical Model, (Sph, see Equation 3) and the Power Model (Pow, see Equation 4). They conform to variogram models already used for ground temperature analysis in an Alpine environment (Kasmaee et al., 2016).
(3) (4)  Where g (h) is the variogram in the vertical direction, h is the lag, C is the sill and a is the range of the variogram. The selected parameters for the variogram models are presented in Table 2, as well as the value of PDN (inhabitants/km 2 ) and GLC (-) for the selected grid. GLC is expressed in codes. Code 190 means "Artifi cial surfaces and associated areas (Urban areas >50%)". All the codes of GLC are presented in Table 3.
The UK method was performed to estimate two target zones: I) For the four boreholes, where data in vertical direction was not available (using the variogram models, infl uenced by the SUHI); II) The average annual temperature below ground surface for all target points of the grid (using the variogram models, without infl uence of the SUHI). To calculate the impact of the SUHI on temperature in depth, the estimated temperature curves were used (see Table 4): • The temperature curve by the measured data (mixed with the estimated data from depth of 0 to 20 m) that are infl uenced by the SUHI: T B (z); • The temperature curve obtained from deterministic function (see Equation 1

Regression analysis
The variogram models and the temperature curves were then studied and compared with each other in order to interpret the temperature gradient deviation caused by the SUHI, and to establish a correlation between the in-dicators (PDN and GLC) and the temperature difference ΔT in different depths.
The sample variograms in all locations and in both conditions (with/without the SUHI) have the same nested structures (spherical model + power model) but with different parameters. Moreover, for all four boreholes, the angle between the temperature curves (with/without the SUHI) is approximately the same. Hence, the ΔT caused by the SUHI could be calculated.
For each selected location (Southeast, South, East and Central), the initial necessary data is: • T 0 -temperature at ground level not inclusive of any SUHI effect; • T N -temperature at investigation depth where the SUHI contribution becomes negligible; • PDN -population density in the grid node; • GLC -global land cover in the grid node. Firstly, angle C between the ground level and the temperature trend was calculated (see Equation 5).

(5)
Secondly, angle A, formed by the intersection of the calculated temperature curve in the natural state -T g (z) -and the estimated temperature trend, infl uenced by the SUHI -T B (z) -, was found to be approximately constant in all four borehole locations and correspondent to 2.5°, but with an intersection located at a different depth. Once known, A and C for each point of the grid, the angle B was easily calculated (see Equation 6).

(6)
The information acquired on the four borehole measurements allowed performing an estimation of the SUHI infl uence in other points of the grid, where only standard calculation of ground temperature, using Equation 1, was possible. In order to do so, it was necessary to use information from Indicators PDN and GLC. In particu-lar, the task was to estimate the value T B (z) at ground level (z = 0 m). The indicator GLC was divided in two groups: urban area (value 190 of Table 3) and rural/forest area (all other values). The indicator PDN was used without changes. A peripheral area was defi ned as a location which, although having a PDN value higher than 0, does not indicate a consistent urbanization. The new Indicator "Population infl uence" D was then calculated. It could assume all values of PDN, but following the rule of Equation 7: Thanks to D, it was possible to take into account the fact that UHI in peripheral areas, with low or very low urban density, is very limited, tending to disappear. Therefore, for the four measurement points (Southeast, South, East and Central), a linear regression was calculated, relating the indicator D and the DT at ground level (DT = T b (0) -T 0 = a), representative of the temperature increase from the natural state due to the thermal impact of urbanization (see Figure 2-left).
The urban infl uence to the underground temperature cannot grow indefi nitely and is expected to smooth after a certain point. Due to the combined effect of UHI and UCI, it is expected that, ground temperature will tend to grow everywhere according to D, but with a limit indicatively set at 20°C, the comfort temperature inside buildings. Therefore, a correspondent threshold of D should be added, after which temperature increase tends to taper off. The same concept can be applied at very low values of D, when the presence of the SUHI is negligible, with DT tending to 0. These thresholds were found by calculating a logarithmic regression (see Figure 2-right) and fi nding the intercepts between the two functions (see Figure 3). Parameters of linear and logarithmic regression are reported in Table 5.
As shown in Figure 3, the lower and upper thresholds were found at D = 700 and D = 3600, since they repre- sent the intersection between the linear and the logarithmic behaviour. Having D in each node of the grid, a correspondent DT can be found, using the linear or logarith-   8, also taking into account UCI contribution. In Southern Europe, the SUHI effect is smoothed since annual average temperatures are higher than 20°C, with subsequent potential UCI. The topic needs additional investigation, however, for the purposes of the present work, a smoothing factor calculated as the difference between annual average temperature and comfort temperature 20°C was subtracted from the T B results.
Once angle B and T B are known at ground level, the temperature in any depth T B (z) can be obtained for each coordinate of the grid (see Equation 9).

(9)
This new temperature value substitutes the theoretical temperature value T g (z), based on the calculation of the geothermal gradient (see Equation 2), until it reaches the interception between two curves, where T B (z) = T N . From that point down, the SUHI effect disappears and standard geothermal gradient should be used to estimate ground temperature.
The result was a matrix of ground temperature values, where each vector includes a node of the grid in one specifi c period of the year, modifi ed from its original theoretical trend by taken into account the SUHI effect.

Results
According to the available information, maps of ground temperature could be potentially created for all of Europe. For the purpose of the present work, the Italian peninsula, including the Alpine territory was selected. Figure 4 shows a map of ground temperature, at a depth of neutral zone, where climate seasonal infl uence disappears, estimation average (see Figure 4-left) and estimation standard deviation (see Figure 4-

right).
As shown in the map, ground temperature at shallow depths is mainly infl uenced by climate, with temperature increasing from north to south and with altitude. Despite this, some "hot spots" are visible, due to geothermal anomalies (as an example: Euganean Hills ancient volcanic zone, near Padua) but also to the UHI effect transferred to the underground. Using the method described in Section 2, consistent temperature differences compared with the surrounding areas have been calculated down to the neutral zone below major cities like Milan and Rome. Due to the uncertainties generated by the model, and by the lack of initial information in many nodes of the grids, a map of estimation standard deviation has been connected to the map of ground temperature. Ground temperature calculations at depths of 20, 35 and 50 m have been selected to evidence the presence / absence of the SUHI effect.
Five cities were selected: Zurich, Milan, Padua, Rome and Palermo. Zurich was selected, being the location of the four borehole measurements. Milan, Rome and Palermo were selected, representing three different climatic zones of Italy, North, Central and South. Finally, Padua was selected being the Italian target city of the GEOTeCH Project, fi nancing this research, with an installed monitoring system of underground temperature (www.geotech-project.eu). An area of 50 km x 50 km around each city centre was chosen to extrapolate the temperature values. Relative frequencies of temperature results at different depths have been calculated and reported bellow for fi ve target cities in Figures 5-9.

Discussion
Relative frequencies shown in the previous paragraph show a general temperature increase from north to south, as expected. The temperature behaviour is not uniform and depends on geological variability, but also on the SUHI effect.
Regarding geological variability, it has an impact in cities like Zurich, Palermo or Padua, while the infl uence is limited for Milan and Rome because of the bedrock depth. In particular, Milan is located in an alluvial plain, and this makes the ground temperature spatially stable at all depths.
Regarding the quality of estimation of ground temperature, frequency values can be compared with standard measurement values for the cities. Regarding Zurich, measurements from Table 4 should be compared with frequencies in Figure 5. The central borehole has a measured temperature from 14.5°C at 20 m to 13.5°C at 50 m. The frequency values are compatible with this measurement, since 10% of data falls in this interval, decreasing in percentage from 20 to 50 m. Since Zurich is a relatively small urban area, with respect to its surroundings, 10% of data presumably refers to the downtown area. Southeast and South boreholes, located in the countryside, range around 12°C. In this case, the frequency of ground temperature correspondent to this value tends to increase from 20 to 50 m. Moreover, the highest frequency corresponds to values below 10°C. This is due to geological and altitude variability. Around the Zurich city area, mountains are vastly present (without borehole measurements) and according to the model, ground temperature is much lower than in the plains (where measurements have been taken). It was calculated that ground temperature around 12°C represents 20% at 20 m depth and around 60% at 50 m depth. Finally, the East borehole, with a ground temperature higher than 13°C from 20 down to 50 m deep, represents a peripheral area, always comprised in the 10% affected by the SUHI effect, even if the contribution is limited.
Regarding Milan, the entire 50 km x 50 km area falls in the urbanized zone of Milan and its surroundings. The ground temperature is very stable. Frequencies are in line with the standard aquifer temperature of Milan (14.5°C), located between 20 and 40 m deep. It seems that the model does not fi t well with the measurements, since in the central part of the city, the temperature should be much higher (17-18°C) than the one estimated (constant around 14°C, caused by UHI) (De Caro, 2018). This is probably due to the hydrogeological conditions typical of underground of Milan, with a very shallow aquifer, with high water fl ow. The model created needs adjustments for this particular situation.
Regarding Padua, the area 50 km x 50 km covers many villages located in the surroundings, as well as the Euganean Hills, an ancient volcanic zone. Ground temperatures are in line with what was measured in the framework of the GEOTeCH Project in the peripheral area (14.5°C), while the higher temperature refers both to zones in proximity of Euganean Hills and the urban area of Padua. The temperature tends to increase with depth, due to the geothermal heat fl ow infl uence occasionally reaching 20°C and more at 50 m depth. The model does not take into account localized hot water rise in fractured zones.
Regarding Rome, higher frequency falls around 17 °C at 20 m depth, which is 2°C higher than annual ambient average. Down to the depth, peaks higher than 23°C are reached at 50 m depth, in a few zones. The contribution of UHI seems evident for the case of Rome, and temperature values are in line with typical ones measured in the aquifer, with values between 18°C and 22°C (Capelli et al., 2008).
Finally, regarding Palermo, higher variability with respect to Rome and Milan can be evidenced, with frequencies ranging from 15 °C to 17.5°C and from 17.5°C to 20°C with approximately the same percentages (40%) at 35 m depth. From 20 to 50 m, temperatures tend to increase, with frequencies in the range from 15 °C to 17.5°C, decreasing from 40% to 30% and with frequencies in the range from 17.5 °C to 20 °C increasing from The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2018, pp. 1-12, DOI: 10.17794/rgn.2018.5.1 30% to 40%. Ground temperatures are in line with the average ambient temperature of Palermo, being around 18°C. In the case of Southern Italy, contributions of UHI are expected to be counteracted by the contribution of UCI, with temperatures not exceeding 20°C, except for a few zones.

Conclusion
This paper presented a mixed probabilistic-deterministic method to estimate ground temperature at a shallow depth, starting from information on geology, hydrogeology, climate, but also urban presence (the UHI effect), through correlations with global land cover and population density. A 4D (x,y,z,t) database of ground temperature was created. Ground temperature estimations below major cities, in the form of frequency distributions on a regular grid, from North to South in the Alpine Area and the Italian Peninsula were compared with ground and aquifer temperature available in literature. The comparison is encouraging, even with some discrepancies with particular reference to cases of high groundwater fl ow. Therefore, still a lot of work is supposed to be done, to improve the precision of estimations. The impact of an accurate mapping of ground temperature affected by the SUHI is potentially high. Firstly, it can support policy makers and stakeholders, bringing a general benefi t for the market of ground source heat pump systems. Secondly, an improvement in guaranteeing the sustainability of the system is expected in the design phase. Thirdly, the understanding of the SUHI effect can also support spreading unconventional ground heat exchangers, which take benefi ts more on heat losses from urbanization than on "real" geothermal energy from the Earth's crust. The presented approach fi ts very well with the topic of circular economy, being the objective to transform a (heat) waste in a (energy) resource.