Estimation of Groundwater Temperatures in Paris, France

Ingolstadt University of Applied Sciences, Institute of New Energy Systems (InES), Ingolstadt 85049, Germany Karlsruhe Institute of Technology (KIT), Institute of Applied Geosciences (AGW), Karlsruhe 76131, Germany University of California San Diego, School of Global Policy and Strategy (GPS), La Jolla, CA 92093, USA Cerema, Direction Centre-Est, 46 rue Saint Théobald, F-38081 L’Isle d’Abeau, France


Introduction
Natural in situ temperatures usually do not substantially vary at a depth of more than 10-20 m. Here, ground and groundwater temperatures are close to the annual mean values in the atmosphere, and commonly only a marginal attenuated influence of the coupled seasonal temperature variation in the atmosphere can be detected [1]. During recent years, attention has been growing to the thermal conditions beneath cities, which were revealed to be completely different when compared to the undisturbed rural surrounding [2][3][4][5][6][7][8][9]. However, the knowledge of the spatial and temporal variability of temperature in urban ground is crucial for proper planning of geothermal energy use. Elevated temperatures offer enhanced opportunities for geothermal heating applications, whereas warmer groundwater is less useful for cooling applications [10][11][12][13][14][15][16][17]. Aside from its role for the geothermal potential, thermal anomalies can influence chemical transport in shallow urban groundwater that often serves as a freshwater resource [18,19]. Finally, anthropogenic accumulation of heat threatens the stability of groundwater ecosystems [20,21].
In most studied cases, temperature in urban ground and groundwater is elevated, which manifests in a large-scale heat carpet underneath a city. This so-called subsurface urban heat island (SUHI) is highly case specific, often with the highest temperatures beneath city centers and strong local spatial variations [7,[22][23][24][25][26][27][28]. Ferguson and Woodbury [29] demonstrated that elevated groundwater temperatures below cities are mainly caused by heat losses from buildings. Accordingly, the investigation of SUHIs enables a determination of the initial period of urbanization in a certain region [30,31]. In addition to building basements, groundwater flow, surface cover, and man-made climate change influence groundwater temperatures [32][33][34][35][36]. Moreover, anthropogenic heat sources like subsurface infrastructure, power plants, landfills, or geothermal installations can have an impact on the spatial distribution of subsurface temperatures [14].
Various studies have reported SUHIs below large cities worldwide, for example, in Moscow [24], London [37], Cologne [7,28,38], Osaka [27,39,40], and Ankara [41]. Within these cities, temperatures are measured in a limited number of boreholes or groundwater wells and are commonly several degrees higher than in surrounding rural soil and pristine groundwater bodies. In the city of Paris, measurements from the 18th century onward are reported in a 28 m deep cellar at the city center [42,43]. The longterm temperature increase of 0.1-0.7 K per decade is attributed to global warming and the developing atmospheric urban heat island (UHI) in Paris. Related studies such as by Pal et al. [44] primarily focus on the atmospheric and surface layer of Paris' UHI. An analysis of surface temperatures shows that the intensity of Paris' UHI, which is the difference between temperatures in the city and undisturbed conditions, can be up to 6 K [45]. Escourrou [46] confirms the presence of a strong atmospheric heat island in Paris and its effect on the local climate (e.g., influence on the rainfall regime and the occurrence of local winds). Cantat [47] underlines that the size and intensity of Paris' UHI depend on various factors such as seasonal variability, climate, or weather conditions. In this work, the focus is set on the shallow subsurface thermal regime in the metropolitan area of Paris. The objective is to examine the SUHI of Paris spatially based on repeated temperature measurements in wells located in and around the city. Due to the low well density within the studied urban area, especially within the city center, we investigate the auxiliary use of satellite-derived land surface temperatures to estimate spatial groundwater temperature distribution. In Method and Data, the measured groundwater temperature data is first described in detail. This information is contrasted to the satellite-based land surface temperatures for the city. Based on this, a new spatial estimation procedure for groundwater temperatures in urban, rural, and mixed environments is explained and applied to the case study of Paris.

Method and Data
2.1. Geographical Setting, Geology, Climate, and Hydrology. Paris is the capital of France and located in the northern part of the country in the region Île-de-France (Figure 1(a)). The focus of investigation here is the entire metropolitan area of Paris, which includes the city center district, the suburbs, and the surrounding rural areas. Paris has a circumference of nearly 55 km and also includes two huge woodlands: Vincennes in the east and Boulogne in the west. The average elevation of the city is 30 m a.s.l. The Seine River with a total length of 13 km inside of Paris crosses the city from northwest to southeast, dividing it into two sections [48].
In general, Paris' city climate is moderate with relatively warm summers and mild winters. The mean annual air temperature is 10.8°C. Temperatures within Paris are slightly higher compared to the suburbs and the surrounding rural areas; the temperature difference typically amounts to 1-2 K. Elevated temperatures within the city center are 2 Geofluids caused by the increasing urbanization (e.g., higher building density) of the capital [47,49]. Paris is eponymous for the Paris Basin (French: Bassin parisien). The Paris Basin is an elliptical geological basin filled with continental, epicontinental, and marine Mesozoic and Tertiary sediments. Postvariscan crustal extension caused subsidence processes during the Permian that formed the Paris Basin alongside other central European basins. It is present in the northeast of France, in the west of Belgium, and in the southeast of England. The center of the basin is located 80 km east of Paris. Aquiferous formations in the Paris Basin are mainly built up by sandstones, limestones, chalk, and carbonates. The aquifer units are intermitted by less permeable claystone and marl layers which results in a complex layered aquifer system [50].

Measured Groundwater Temperatures.
There have been several groundwater measurement campaigns in the metropolitan area of Paris. Data from 1990 to 2015 were measured by local authorities in wells of the Paris metropolitan area and are available in the French national database on groundwater resources [51]. The well coordinates and measuring depths are accessible at InfoTerre [52]. Even though a rich dataset of 3018 GWT measurements exists, the records are heterogeneous, often incomplete over time and information on the specific sampling methods is sparse. The majority of GWT data were measured on site while pumping. Therefore, the given depths represent the mean depth of the screened well interval. During preprocessing, we excluded outliers exceeding a 6 K offset to the median for a single well location for further analysis. Moreover, measured values with an offset above 2.5 times the interquartile range (IQR) were expelled if the IQR was above 1 K. This was necessary to exclude outliers most probably caused by poor handling during sampling and typos. In addition, four outliers that are directly influenced by a waste disposal site north of Paris were eliminated. A total of 74 outliers (2.3% of the available data) were removed from the data pool. The remaining data covers 377 sampling wells in the outer and inner suburbs of Paris. The derived spatial coverage is nonuniform and there are no data points in the city center ( Figure 1).
Diffre et al. [53] provide GWT data that include measurements in the city center (Figure 1(b)) from February 1977. As the corresponding coordinates of the observation wells were not given in that report, they had to be extracted from the online database InfoTerre [52]. Unfortunately, none of the wells measured in 1977 are listed in the ADES [51] database ( Figure 1). The 1977 GWTs focus on the inner suburbs and the city center, while the more recent ones reveal the thermal conditions in the surrounding.

Satellite-Based Estimation of Groundwater Temperatures.
Recently, satellite data was suggested to investigate the thermal ground conditions in cities and SUHIs. Zhan et al. [54] estimated ground temperatures in Beijing, China, from MODIS land surface temperatures. The results were compared with recorded temperatures (measuring depths of 0.05 m, 0.40 m, and 3.20 m). A time delay between the maximum temperatures of the atmosphere and the ground was detected depending on the measuring depth. Hafner and Kidder [55] used Advanced Very High-Resolution Radiometer (AVHRR) satellite data to determine the albedo of the land surface as well as soil thermal and moisture properties around Atlanta, USA. These parameters were utilized for a numerical simulation of Atlanta's atmospheric temperature dynamics. The spatial extension and intensity of atmospheric urban heat islands in 18 Asian megacities were analysed by Tran et al. [56] using satellite-derived land surface temperature (LST) data. This study demonstrated that the spatial heterogeneity of the investigated heat islands mainly depends on vegetation cover and building density, whereas its total extent and intensity are related to the population size.
Benz et al. [38] analysed whether the spatial distribution of interpolated GWT is linked to the spatial distribution of satellite-derived LST. In this context, the four German cities Berlin, Munich, Cologne, and Karlsruhe were analysed. In order to quantify the relation between GWT and LST, the Spearman correlation coefficient was determined and correlations of up to 80% were found. Benz et al. [38] also propose a method for calculating groundwater temperatures with the help of mean annual LST, building density (BD), and basement temperature (BT). This reflects that accelerated ground heat flux in cities is caused not only by land surface changes and modified above-ground temperatures as detected from satellites but also by other heat sources such as heat release from basements, buildings, and underground networks.
Even for undisturbed rural areas, LSTs are generally colder than adjacent GWTs. This offset is caused by a variety of processes in the soil transition zone and varies with local climate, biological activity, and the hydrogeological setting [57,58]. Benz et al. [59] introduce an approach to calculate this offset on a global scale. This approach empirically relates evapotranspiration (ET) and the number of days with snow cover (SD) to the offset between LST and GWT. For a global GWT dataset, they could reduce the root mean squared error (RMSE) of prediction by 0.5 K to 1.4 K and increase the squared Spearman correlation coefficient by 0.5 to 0.95 compared to a prediction relying solely on LST.
Due to the fact that we have a strong spatial variety of our measured GWT in Paris (Figure 1), in the following method, we combine the individual approaches of Benz et al. [38] and Benz et al. [59] to estimate GWT in the region of Paris, both for rural and urban areas at the same time ( Figure 2). This combined approach has been suggested by Benz [60] for German cities, and it has the specific advantage of covering natural and urban processes within one estimate.
Estimated groundwater temperatures (eGWT) are calculated in three steps: eGWT rur = LST + 3 5 ± 0 2 · 10 4 K m 2 · s kg · ET + 6 6 ± 0 3 K · SD, 1 eGWT urb = max LST, The first step consists of estimating GWT in rural areas. Here, the offset between LST and GWT is dominated by the latent heat flux cooling down the surface in the form of evapotranspiration and by snow cover insulating the subsurface during the winter period. An empirical best-fit approximation for estimating rural groundwater temperatures (eGWT rur) from satellite-derived data on a global scale is given by equation (1) [2]. The second step is estimating groundwater temperatures in urban environments (eGWT urb ) by using equation (2), which takes into account building density and basement temperatures [38]. The third step combines rural and urban estimations by replacing LST by eGWT rur in equation (2). The resulting equation (3) is intended to estimate groundwater temperatures in rural, urban, and mixed settings (eGWT c ).
Building densities were calculated at a resolution of 1 km × 1 km from OpenStreetMap [61] building polygons. Basement temperatures were estimated to be 17 5 ± 2 5°C following guidelines by the International Organization for Standardization (ISO) [62] and assuming similar conditions for Central Europe. This value range is also consistent with average basement temperatures obtained from borehole temperature profile inversion in the area of Zurich, Switzerland [30].
Satellite-derived data represent the decadal mean from 2005 to 2014 and include land surface temperatures, evapotranspiration, and the percentage of snow days. Data was retrieved and processed with Google Earth Engine according to the procedure that is explained in detail in Benz et al. [59]. The MODIS Aqua and Terra Land Surface Temperature and Snow Cover Daily Global 500 m products were retrieved from Google Earth Engine, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota [63,64]. In contrast to Benz et al. [59], the updated MODIS data version 6 was used for determining LST values. Snow days were computed with MODIS version 5 due to different snow cover algorithms and product availability in the version 6 data, which could not be used with the empirical variables in equation (1) without biasing eGWT rur . For a detailed discussion of the given raster data, we refer to Benz et al. [59]. Estimated groundwater temperatures were computed by equations (1), (2), and (3) using ArcPy and ArcMap 10.5.1 [65]. Note that extracted raster data on point locations were bilinearly interpolated between the middle of adjacent raster cells. All raster data was exported at a resolution of approximately 1 km × 1 km (Figure 3).
Land surface temperatures (LST) range from 10.4°C in the outer suburbs to 14.9°C within the city center (Figure 3(a)). This dataset exhibits two central LST cold spots (12.7 and 12.5°C), caused by the two large forests in the eastern and western parts of the city center. The annual percentage of snow days (SD) ranges from 0 to 7.3% and as expected SDs are more frequent in the outer suburbs than in the city center. Again, the two forests manifest as anomalies in the spatial SD distribution (Figure 3(b)). In contrast to LST and SD, the data on ET has a lower resolution of 0 25°× 0 25°(approx. 28 km × 28 km). The city center is almost fully covered by one cell of low ET with a value of 10 mg m -2 s -1 . Towards the rural areas, ET increases up to 17 mg m -2 s -1 (Figure 3(c)). Finally, building density (BD) is calculated from building polygons of OpenStreetMap at a resolution of 1 km × 1 km. BD is highest in the city center (up to 60%) and decreases towards the suburbs to nearly zero in areas of low population. 4 Geofluids to evaluate the misfit of the predicted GWT. Calculating RMSE is probably the most common way to evaluate the predictive capabilities of a model. However, it is more sensitive than MAE to large prediction errors or outliers, because errors are squared before being averaged. ME is used to describe the error orientation, as positive and negative errors cancel each other out when calculating ME. It is thus useful to test if observed values are over-or underestimated by the prediction. In addition, a linear least square regression is performed to correlate measured and estimated GWTs by using the Spearman method. Hereby, the obtained Spearman correlation coefficients (r) quantify the link between estimated and measured GWTs, but are less useful as (mis)fit measures.

Spatial and Vertical Distribution of GWT.
For the first step, we characterize the spatial and vertical distribution of GWTs based on the measured data ( Figure 4). For this, the available early measurements from 1977, which focused on the city center, and averages from the years 1990 to 2015 for the greater Paris region were spatially combined. In February 1977, the reported GWT beneath the city center ranged between 14.1°C and 18.3°C. Towards the suburbs, the measured temperatures decreased yielding a temperature range from 12.2°C to 16.6°C. The measurement depths range from 17 m to 100 m below ground level. For the shallower measurements, a slight seasonal bias can be While no clear spatial trend can be deduced for the shown data, the lower limit of GWT clearly declines towards the outskirts. Nevertheless, the lower limit of 14.1°C for the groundwater temperatures in the city center in 1977 is remarkable, being more than 3 K higher than the average air temperature of 10.8°C during that time. On average, they are also warmer than temperatures in the inner and outer suburbs. If we ignore any temporal and or vertical effects, the maximum temperature in the city center (measured in 1977) exceeds rural groundwater temperatures (measured after 1990) by 5 to 6 K. This offset is typical for SUHIs in Western Europe and is on the upper end of observed intensities in previously studied German cities [7].
The vertical resolution of the measured data is depicted in Figure 4. Most sampling points were obtained in shallow depths of a few tens of meters, with some measurements also taken in depths of more than 100 m. However, no clear evidence of the geothermal gradient can be observed. Still, for the comparison between satellite-derived LST and the measured GWT, the lower depth limit for this analysis is set to 80 m in order to zoom in on the regime influenced by land surface effects. Additionally, the shallowest 15 m of the subsurface are disregarded to avoid any influence of seasonal atmospheric temperature variation. The usual seasonal temperature variation in a humid climate is below 1 K for depths of more than 10 m [1].
In the next step, equation (3) is applied to estimate groundwater temperatures using the combined method, eGWT c . This resulting map in Figure 5(b) represents eGWT c at shallow depths with a temperature range of 11.3°C to 16.6°C. Consistent with the available groundwater temperature measurements from different points in time, the city center is characterized by warm temperatures, which are decreasing towards the suburbs.

Comparison of Measured and Estimated GWT.
To assess the quality of predicted groundwater temperatures from satellite imagery, eGWT c , they are compared with well measurements in the field. It is important to note the time difference: data in the city center was only observed in 1977, about three decades earlier than the satellite data considered for calculating eGWT c . However, as GWTs are expected to have increased during that time [30,66,67], these measurements are included in the following comparison viewing them as minimum values. In order to inspect the role of temporal variability and trends, we separately compare measured GWTs from (i) 1977, (ii) 1990-2015, and (iii) 2005-2015 with satellite-derived estimates. The latter GWT subset is extracted as it covers roughly the same time period as the satellite-based product. Since satellite data covers the time period from 2005/01 till 2014/12 and groundwater temperatures are recorded till 2015/11, the data contains a negligibly small amount of 85 unsynchronized single measurements (4.4% of the single measurements between 2005 and 2015). Figure 6 shows the measured GWTs versus LST and eGWT c along with misfit and Spearman correlation indices. The city center wells measured in 1977 (Figure 6(a)) yield misfits (RMSE, MAE) between eGWT c and measured GWTs around 0.5 K higher than the error observed between LST and GWT. Additionally, ME is higher (Figures 6(a) and 6(d)) when comparing measured groundwater data with eGWT c (1.27 K) than with just LST (0.09 K). This means that subsurface temperatures are overestimated by eGWT c when using recent values of LST and old well temperatures in equation (3). This is because LST is expected to have increased by around 1 K from the 1970s to the period between 2005 and 2015 due to climate change and atmospheric urban heating, which is not surprising. At the same time, GWTs in 1977 were as high as 18.3°C, which is warmer than the assumed basement temperature of 17.5°C in equations (2) and (3). Therefore, other heat sources such as heat release from subsurface infrastructures may play a significant role but cannot be reproduced by this method. This is supported by the low sensitivity of the measured GWT to LST as revealed by Figure 6(a).  7 Geofluids being slightly lower compared to the entire dataset from 1990 to 2015. We can therefore conclude that the estimation approach is most applicable for overlapping time periods between subsurface and satellite measurements.
Overall, the combined approach as given in equation (3) delivered groundwater temperature estimates with a RMSE below 1 K. However, even for this time period, extreme temperature anomalies cannot be predicted. It is very likely that these wells are impacted by local heat sources other than buildings. When applying this method to the German cities of Berlin, Cologne, and Karlsruhe, Benz [60] found similar results with RMSEs between eGWT c and GWT of 0.96, 0.86, and 0.82 K, respectively.

Impact of Urban and Rural Estimation Approach on
GWTs. To evaluate the suitability of the different GWT estimation procedures (equations (1), (2), and (3) and LST only), we use the correlation coefficient (r), RMSE, MAE, and ME. Figure 7 shows the correlation to the measured GWTs of the different equations for calculating eGWT urb , eGWT rur , and eGWT c and the sole use of LST.
While RMSE and MAE are lowest for the combined method, indicating an improvement of the estimation accuracy, correlation is highest for eGWT rur . According to these results, this approach which estimates naturally occurring offsets also has an impact on the predicted groundwater temperatures which is higher than that of the urban approach (eGWT urb ). This is a consequence of the location of the observation points that are solely in the inner and outer suburbs of Paris and have BD values below 0.3 (30%). In the city center, BD values are more than doubled, which also doubles the amount of positive temperature correction (equation (2)). At the same time, eGWT rur has a lower impact on the city center because both SD and ET are significantly lower (Figure 3).
All predictions underestimate measured GWTs with a minimum ME of -0.23 K for eGWT c . However, we could improve the negative offset by 0.9 K compared to the misfit between measured GWTs and LST. In contrast, Benz [60] found that the combined method overestimates temperatures in the city centers of Karlsruhe, Cologne, and Berlin. This overestimation can in part be attributed to the insulation of basements, which is not considered in this estimation. Overall, we expect eGWT c to work best in the city centers of larger urban areas, due to higher BD and lower ET and SD values. For the given groundwater data in the outskirts of Paris, eGWT c is only slightly favorable over the use of eGWT rur .
Although the comparison between the results of the individual equations correlates well with measured groundwater, the used approach also has some limitations. Particularly, local thermal disturbances caused by underground anthropogenic heat sources (e.g., geothermal devices, sewer networks, and underground constructions) cannot be detected from above-ground data or their thermal footprints are too small to be resolved by the resolution of 1 km × 1 km of eGWT c . Moreover, the resolution of the ET dataset is very coarse (about 28 km) and cannot detect small-scale or even city-scale features. Applying the combined method 8 Geofluids to Karlsruhe and Cologne in Germany, the ET dataset did not display the reduction of ET commonly associated with urban areas, which added to an overestimation of GWTs in those areas [60]. Unfortunately, no recent GWTs within the Paris city center are available here, and these observations can thus not be tested for such a large-scale city. However, future measurement campaigns could help to eliminate this limitation. Another limitation of the method used is that LST as primary input data already indicates a heat anomaly in the city center, and the data used for ET and SD are also affected by urbanization. Accordingly, it cannot be differentiated between natural and anthropogenic thermal contributions by the use of this method, whereas the processes itself can fairly easily be attributed to one or the other.

Conclusions
For the first time, the large-scale subsurface urban heat island (SUHI) in the groundwater of Paris is spatially investigated. The minimum and maximum groundwater temperatures (GWTs) are 11.4°C and 17.2°C measured between 2005 and 2015 for the shallow subsurface at depths between 15 m and 80 m. These groundwater temperatures represent the inner and outer suburbs of Paris. They should be interpreted as regional trends and are not meant to represent local temperature disturbances caused by anthropogenic heat sources which may be of a much greater magnitude. The minimum and maximum GWTs observed in 1977 below the city center ranges between 14.1°C and 18.3°C. Despite the temporal gap of 28 years between the two measurement campaigns, this indicates a characteristic difference between the center and the suburbs. Based on the examined data, the maximum GWT anomaly reaches 6.9 K in Paris.
In the present study, an approach is applied that estimates groundwater temperatures (eGWT c ) from satellite-derived data and building footprints combining existing methods for estimating rural (eGWT rur ) and urban (eGWT urb ) groundwater temperatures. Estimates are validated with temperatures measured in the wells of Paris. Predicted GWTs have a very good fit (RMSE < 1 K) and only slightly underestimate measured GWT with a mean error of -0.23 K. In addition to that, a coherent spatial estimation of the SUHI of Paris is obtained. Even so, measured extreme GWTs that cannot be reproduced as subsurface heat sources other than in basements are not considered. In an ideal case, the estimation techniques displayed here would supplement existing temperature measurements in wells, which would allow the use of this method as a tool to spatially connect point information obtained from wells.
In future studies, the presented procedure needs to be further validated by estimating GWTs in other rural and urban environments. This is ideally accomplished in regions, where GWT data with a very dense spatiotemporal resolution is available. When applying this approach to other cities, it should be kept in mind that the anthropogenic fingerprint varies among different regions, countries, and climates. The role of this variability will be scrutinized by a comparison of different case studies in the future.

Data Availability
Groundwater temperature data is available in the French national database on groundwater resources [51].

Conflicts of Interest
The authors declare that there are no conflicts of interest.