Urban storage heat flux variability explored using satellite, meteorological and geodata

The storage heat flux (ΔQS) is the net flow of heat stored within a volume that may include the air, trees, buildings and ground. Given the difficulty of measurement of this important and large flux in urban areas, we explore the use of Earth Observation (EO) data. EO surface temperatures are used with ground-based meteorological forcing, urban morphology, land cover and land use information to estimate spatial variations of ΔQS in urban areas using the Element Surface Temperature Method (ESTM). First, we evaluate ESTM for four “simpler” surfaces. These have good agreement with observed values. ESTM coupled to SUEWS (an urban land surface model) is applied to three European cities (Basel, Heraklion, London), allowing EO data to enhance the exploration of the spatial variability in ΔQS. The impervious surfaces (paved and buildings) contribute most to ΔQS. Building wall area seems to explain variation of ΔQS most consistently. As the paved fraction increases up to 0.4, there is a clear increase in ΔQS. With a larger paved fraction, the fraction of buildings and wall area is lower which reduces the high values of ΔQS.


Introduction
The storage heat flux (ΔQ S ) is the net flow of heat stored within a volume that includes the air, trees, buildings and ground. The ability to absorb, store and release heat depends on the thermal mass and morphology. In urban areas, the net heat stored in the canopy is a relatively large fraction of the net all-wave radiation (Q*) compared to other environments (Nunez and Oke 1977;Grimmond and Oke 1999). In highly urbanised areas, it can account for more than half the daytime net all-wave radiation ) and be two to ten times larger than for simple planar surfaces (e.g. soil). The wellknown nocturnal urban heat island (UHI) is caused by the release of stored heat and enhanced by anthropogenic heat (Q F ). Combined with reduced radiative cooling (or enhanced radiative trapping), the storage heat flux is a major contributor (Oke and Cleugh 1987). Oliphant et al. (2018) demonstrate the importance of building materials such as concrete and asphalt as an essential factor to enhance ΔQ S , as increased surface roughness using light-weight materials neither affect the storage term nor the UHI. As nocturnal cooling is important for recovering from daytime heat stress (Rocklöv et al. 2011;Thorsson et al. 2014), the expected increases in both urban population (UN 2015) and heat wave frequency (Schär et al. 2004) will likely cause increased heat stress and heatrelated morbidity and mortality.
In simple environments, the storage heat flux can be directly measured using heat flux plates buried a few centimetres below the surface with temperature sensors above to determine the flux divergence. However, in complex urban landscapes, this approach is impractical at the local scale. There are a range of methods to assess the storage heat fluxes in urban areas, including OHM, Objective Hysteresis Model Grimmond and Oke 1999); AnOHM, Analytical Objective Hysteresis Model (Sun et al. 2017); RES, Residual, determination of the storage heat flux from the residual of the surface energy balance (Offerle et al. 2005b); CAR, Complete Aspect Ratio (Rigo and Parlow 2007); TEB, Town Energy Balance model (Masson 2000) or other urban land surface models; and ESTM, Element Surface Temperature Method (Offerle et al. 2005a). Some methods (e.g. OHM, CAR) use bulk parameters by material types, whereas other methods (e.g. AnOHM, TEB and ESTM) require the thermal parameters (e.g. heat capacity) for the component materials. AnOHM provides a method to determine OHM parameters.
Studies exploiting Earth Observation (EO) data to derive spatial variations of ΔQ S are very sparse. Rigo and Parlow (2007) make use of the normalised difference vegetation index (NDVI) and net all-wave radiation (Q*) to obtain ΔQ S . Kato and Yamaguchi (2007) exploit the Advanced Spaceborne Thermal Emission and Reflection radiometer (ASTER) sensor system to derive ΔQ S as a residual from the urban energy balance. However, they do not separate the anthropogenic and storage heat flux terms.
In this paper, we use EO data to estimate spatial variations of ΔQ S in urban areas using the ESTM scheme. ESTM (Section 2) accounts for variations in urban morphology, land cover and land use. We evaluate ESTM at four sites with different land covers (grass, deciduous trees, asphalt and an urban canyon) with detailed observations available. We couple ESTM-SUEWS (Section 2) and use this system to address the spatial and temporal variability of ΔQ S in three cities in 2016 : Basel (Switzerland), Heraklion (Greece) and London (UK). As clear skies are required to acquire satellite based surface temperature data, the full temporal range cannot be assessed.

ESTM
The Elemental Surface Temperature Method (ESTM) (Offerle et al. 2005a) reduces the 3-dimensional urban volume to four 1-dimensional elements (i.e. building roofs, walls, internal mass and ground (road, vegetation, etc.)). The storage heat flux is calculated from element (i) surface temperatures (T i ): where ΔΤ i /Δt is the rate of temperature change over the period for each element i, ρc is the volumetric heat capacity, Δx i is the element thickness and f i is the plan area index of that element.
So, x i f i is simply the total element volume over the plan area, for each element i. The element layers (e.g. wall brick, insulation, wood) average internal temperatures are accounted for, with: where Q is the heat flux through the surface and k is the thermal conductivity. The surface temperature of internal building elements (floors, ceiling and internal walls) is determined from setting the conductive heat transfer out of (in to) the surface equal to the radiative and convective heat losses (gains), as described by Offerle et al. (2005a).
To facilitate ESTM usage, the scheme is incorporated into the Surface Urban Energy and Water Balance Scheme (SUEWS) (Järvi et al. 2011(Järvi et al. , 2014Ward et al. 2016;Järvi et al. 2019). This simulates the urban radiation, energy, water and CO 2 fluxes with each grid characterised by the fractions of seven surface types: paved (e.g. roads, sidewalks), buildings, evergreen trees/shrubs, deciduous trees/shrubs, grass, bare soil and water. At each time step, both the surface water state  and the soil moisture below each surface type (excluding water bodies) are calculated. To force SUEWS, the minimum meteorological data required are downward shortwave radiation, wind speed, (outdoor) air temperature, relative humidity, atmospheric pressure and precipitation (Table 1).
3 Evaluation for the heat storage for simple surfaces

Methods
Given the difficulty of measuring storage heat flux in complex urban areas, we evaluate the performance of ESTM for individual components of the urban environment. The four sites have single land covers: asphalt surface in Säve, near Gothenburg, Sweden, day of year (DOY) 43-106 (Jansson et al. 2006); long grass site in Basel, Switzerland, DOY 197-327 (Parlow et al. 2014); street canyon (Torggatan) in Gothenburg, Sweden, DOY 1-213 (Offerle et al. 2007) and a deciduous forest site, Morgan-Monroe State Forest (MMSF), USA, DOY 60-365 (Oliphant et al. 2004). Material properties for the different sites are given in Table 2.
The meteorological forcing data are compiled from nearby weather observation sites. Evaluation data for ground heat flux are derived from heat flux plates (Säve and Basel). Evaluation data from Torggatan (Gothenburg) where compiled from unshielded fine-wire thermocouples (TC) (Omega, T-type, 0.127 mm) measured surface temperature affixed to the surface facets with a thin layer of adhesive including indoor temperature observation using Tinytag sensors (Offerle et al. 2007). Evaluation data for the deciduous forest (MMSF) are soil, air and biomass storages heat fluxes determined from soil heat flux plates, thermocouples and ventilated thermistors (Oliphant et al. 2004).

Evaluation results for individual surfaces
The ESTM scheme can satisfactorily estimate ΔQ S for the four test sites (Fig. 1). The best performance is for the grass (mean absolute error (MAE) = 5 W m −2 ). The deciduous forest and asphalt MAE are 16 W m −2 and larger for the urban canyon (MAE = 49 W m −2 , Fig. 1). The mean bias error (MBE) is < 1 W m −2 at the three sites and 22.4 W m −2 for the urban canyon.
Of the four areas, two are relatively simple (grass and asphalt) and two are very complex (deciduous forest and urban canyon). In the latter cases, the detailed measurements allow the 3D environment influence on storage heat flux to be assessed. For example, the total ΔQ S of the deciduous forest includes contributions from air, leaves and branches. One explanation for the high accuracy of ESTM for the simple surfaces is that the parameters needed such as thickness, volumetric heat capacity and thermal conductivity can be set with a high accuracy as the material properties for each site can be studied in detail or derived from observed temperature profiles.
Although the more complex deciduous forest and street canyon sites have the larger scatter ( Fig. 1), ESTM is able to capture the variations of ΔQ S in a fundamental descriptor of the city-the street (or urban) canyon (Oke et al. 2017). The uncertainty in observed ΔQ S in more complex environments is greater, because of both the large number of measurements involved and assumptions required (Oliphant et al. 2004, Roberts et al. 2006, Offerle et al. 2007).

Sites and meteorological forcing data
The three H2020 UrbanFluxes project  cities are the focus of this study. The cities range in size from the mega-city of London (UK), to medium-sized central European city of Basel (Switzerland), to the small low latitude Mediterranean city of Heraklion (Greece). For each, the central part of the city and some vegetated areas are included ( Fig. 2) in the (west-east × north-south) model domains ( Fig. 2): Basel, 5.1 km × 4.9 km; London, 21.5 km × 21.4 km; and Heraklion, 13.2 × 6.8 km. As Heraklion is a much smaller city, the domain extends out to the surrounding rural area (Fig. 2).
As continuous forcing data (Table 1) are needed for both SUEWS-ESTM and to permit the net change of storage heat through time, the simulation time step should be 1 h or less. SUEWS-ESTM forcing data may come from observations (e.g. meteorological towers) or larger-scale models (e.g. meso-scale model or re-analysis data). Here, we use data from instruments installed on meteorological towers ( Fig. 2) (Crawford et al. 2017;Feigenwinter et al. 2018;Stagakis et al. 2019).
It is assumed the internal building element temperature is mainly controlled by the internal air temperature (T iair ). This is modelled following Georgitsi (2011), with a sinusoidal variation around a base indoor temperature (T base ) assumed to be at a minimum at 04:00 and a maximum at 16:00 T base is increased (decreased) as outdoor air temperature (T a ) increases (decreases). Time of day (t day ) is expressed in decimal hours. The resulting diurnal range in T iair is typically within a 1-5°C.

Land surface temperature
Landsat 8 and MODIS Terra satellite data, resampled to 100 m resolution (Mitraka et al. 2015), are used to retrieve the land surface temperature (T LST ). For Landsat 8, the thermal infrared sensor (TIRS) bands surface reflectance are used with the ATCOR algorithm (Richter and Schläpfer 2015) assuming a constant surface emissivity (0.98) and mid-latitude atmosphere.
MODIS Terra (1 km × 1 km resolution) TIR bands top-ofatmosphere radiance are downscaled with a spatial-spectral unmixing method (Mitraka et al. 2015). The spectral atmospheric correction uses ATCOR. The surface spectral emissivity is estimated by determining surface cover fractions from the high-resolution visible and near-infrared (VNIR) Landsat data combined spectral libraries (Kotthaus et al. 2014).
The satellite data images are acquired before (morning) and after (evening) the peak surface temperature at times that vary The  Table 3.
As ΔQ S is the net change in heat stored per time from changes in the surface and internal material temperatures, to use ESTM with instantaneous satellite data, a continuous time series is needed. In the morning, a sinusoidal relation between the outdoor air temperature (T a ) and surface temperature (T S ) difference is assumed (Lindberg et al. 2008: where T s -T a is assumed to be 0 K at sunrise (SR) and marks the start of the sine period. The phase is modified and dependent on time of sunrise (t SR ) and time of maximum in T s (t peak ). The surface temperature observations from the four evaluation sites (Section 3) were used to obtain the timing: where t SS is the time of sunset. With no satellite acquisition, the amplitude (a) is calculated as a function of maximum sun elevation angle (α smax ), as described in Lindberg et al. (2016). When satellite data are available, both a retrieved T LST and the satellite overpass time (t) are known. In a second step, a continuous T s is calculated.
The T s decrease in the afternoon and evening after t peak has a more (cf. to the morning) complicated pattern as the cooling rate reaches a maximum and then levels off ). To derive the cooling pattern, the observed surface temperature at the four evaluation sites (Section 3) are analysed (Fig. 3). The common pattern in the surface temperature decreases follows the NOCRA (NOcturnal Cooling RAte) model ) except that the maximum surface temperature cooling rate appears earlier. The different cooling phases (1a, 1b and 2, Fig. 3) are described using sine, cosine and linear fits, respectively. Onomura et al. (2016) provides details.
The Torggatan street canyon data are used to derive the built afternoon and night surface temperature cooling parameters. Phase 1a starts at t peak (i.e. cooling rate is zero) and ends at the time of maximum cooling rate (t maxcool = t SS − 0.08(t SS − t SR )). Phase 1b continues until the start of Phase 2 (t 2,start = t SS + 1.5). Phase 2 ends at sunrise the next day (i.e. cooling rate is zero). The cooling rates during the three phases are: When satellite T LST are available, the temperature rate amplitude (A r ) can be retrieved in a similar way (Eq. 4) to the morning surface temperature model.
The evening surface temperature profile is calculated from the daytime peak surface temperature by integrating evening surface temperature rate over time. If no evening satellite data are available, the morning scheme is used until T s drops below T a . It then stays at T a until t SR the next day. This permits the storage heat flux modelling to continue without satellite data.

Surface parameters from geospatial data
For each city, SUEWS-ESTM is run with 100 m × 100 m resolution. The input parameters for the models (Table 1) are prepared using Urban Multi-scale Environmental Predictor (UMEP) .
High-resolution (e.g. 1 m) geospatial datasets, derived from EO data using advanced machine learning techniques and detailed spectral mixture models (Mitraka et al. 2016;Marconcini et al. 2017), are used to derive both land cover fractions and other morphological parameters (e.g. wall height, wall area and frontal area index). The digital surface models (DSM) either include both ground and building heights or only building heights above ground. In the former case, digital elevation models (DEM) of ground heights are used to obtain relative heights of object. For Heraklion, the DSM are derived from very high-resolution optical stereo imagery and for Basel and London, airborne LiDAR observations are exploited (Marconcini et al. 2017;Lindberg and Grimmond 2011).
Urban areas are often described using a street canyon (Nunez and Oke 1977) with a mean building height (z H ) and street width (W). The real 3-dimensional urban morphology is simplified into a 1-dimensional infinitively long street canyon with roof, wall and ground facets. To ensure conservation of heat and momentum, the 3D to 1D transformation (Lindberg et al. 2015) used here is the Martilli (2009) approach. The fractions of the three canyon facets are set to be the same as the real morphology, so that: where f wall is the fraction of the wall area relative to the total horizontal area. For details, see Martilli (2009) or Lindberg et al. (2015. The urban form parameters are derived from highresolution DSMs (Table 1). To derive f wall , a 4-directional 3 × 3 kernel majority filter on the DSM is applied. Differences between the original DSM and the raster produced from the filtering are identified. A threshold is set for a wall height (e.g. ≥ 3 m) allowing wall pixels to be identified. f roof is derived from high-resolution ground and building DSM in conjunction with a ground only DEM.
The fraction of internal building surface elements (f ibld ) depends on fractions of wall (f wall ) and roof (f roof ), mean building height (z H ) and the number of rooms per floor (n room ). An idealised indoor building geometry is assumed with two rows of equally sized rooms separated by a corridor on each floor. From this geometry, f ibld is: where z floor is the floor height (3.1 m used). In the last term, − 1 is used to exclude the outer roof. With a small number of rooms per floor, f ibld increases rapidly but as the number grows so does the wall fraction. Beyond 10 rooms per floor, the change of the contribution of internal building surface to the total urban surface area is small. The morphometric parameters can be derived using vector data (e.g. polygon building footprint data) also. Although vector data allow situations such as two attached buildings with different roof heights to be better resolved, these conditions are proportionally extremely rare. Furthermore, a direct conversion of linear vector walls will result in an overestimation of wall areas (Lindberg et al. 2015). For these reasons, a raster dataset is used in this study.
For the thermal parameters (Table 2), both land cover (e.g. buildings, paved, bare soil) and land use (e.g. residential, industrial, agricultural. areas) are considered. ESTM treats evergreen trees/shrubs, deciduous trees/shrubs, grass, bare soil and water as having constant thermal values across the city with variations in phenology and soil moisture not considered in this study.
Paved and building land cover classes are sub-divided into three and five land use classes, respectively. The Urban Atlas (EAA 2017) is used to separate roof types (e.g. suburban and city centre may have ceramic tiles and concrete roofs, respectively) and wall characteristics (e.g. fraction of glazing, insulated or not). Manual ground inspections, and comparison with Google Satellite View (ground and roofs) and Google Street View photography (walls) provide external information (Google 2016). The element layer attributes (Table 2) are based on typical construction practices.

Spatial storage heat fluxes in three cities
Storage heat fluxes are calculated for the three cities using all available satellite images (> 50% clear sky) for 2016 in Basel (206), Heraklion (300), and London (142). As both morning and evening satellite data are available on 2016 July 19 for all sites, we select this day to show the storage heat flux maps (Fig. 4). In the morning, storage heat fluxes have large (positive values) in the dense building areas indicating warming of the surfaces. In London, and to a limited extent in Basel, tall buildings with a big volume for heat storage have large fluxes. This is apparent in the eastern part of the City of London and further east in the Canary Wharf business districts (further east) where buildings are 200 m and taller (cf. Figs. 2 and 4). In Basel, there are a few scattered tall buildings (generally < 80 m) and in Heraklion the building height rarely exceeds 30 m (Fig. 2). Extensive vegetated areas, especially where trees are present (e.g. parks, Fig. 2), stand out with low ΔQ S . The road network, most discernible in London, is where intermediate (~150 W m −2 ) size ΔQ S values occur. Water bodies (e.g. Rhine and Thames) are not well represented by ESTM. Generally, in the evening, the areas that stored most heat during the day release (negative values) the most heat (Fig. 4).
The range of ΔQ S varies substantially between the cities on this date (2016 July 19), with Basel having both extremes from~− 340 to 400 W m −2 cf. Heraklion, − 190 to 200 W m −2 and London, − 200 to 300 W m −2 . The larger range in Basel, compared to London with its taller buildings with higher thermal mass, is caused by big differences between air and surface temperature (18-32°C) on 2016 July 19 (Fig. 5). The storage heat flux depends on the surface temperature, which varies with the incoming shortwave radiation and the resulting outdoor air temperature (Eqs. 1 and 2) (amongst other things). The air temperature is quite different between the three cities and the timing of the satellite overpasses relative to the air temperature change through the day (Fig. 5).
The magnitude of the storage heat flux is in principle dependent on the thermal mass (e.g. fractions of buildings, paved and vegetated areas, height and density of buildings, types of material) and the morphology of the urban setting (i.e. sky view factor). These relations are investigated for four key parameters, (i) mean building height (z H ), (ii) wall area, (iii) building fraction and (iv) paved fraction. In Fig. 6, all summer (June, July and August) morning satellite acquisition storage heat fluxes (Table 3), retrieved from the ESTM model, are presented. The fluxes are normalised by the measured incoming shortwave (K ↓ ) radiation and modelled incoming longwave radiation (L ↓ ) for each satellite overpass (ΔQ S / (K ↓ + L ↓ )). London has considerably higher ΔQ S /(K ↓ + L ↓ ) than the other two urban areas. The overall pattern between the different measures of surface characteristics and ΔQ S / (K ↓ + L ↓ ) is similar for the three study areas. Building fraction, z H and wall area have a linear pattern. There is a peak in ΔQ S / (K ↓ + L ↓ ) at around 0.4 in paved fraction across all three cities. This is consistent with Loridan and Grimmond (2012) analysis of eddy covariance and surface energy balance closure data for multiple sites around the world. Wall area is the surface characteristic which shows the least scattered ΔQ S /(K ↓ + L ↓ ). This is also evident for all three study areas.
Neither building fraction nor z H provide the complete 3D information of the urban area. For example, a large fraction of buildings may include a few extensive buildings (e.g. warehouses) with small areas of walls (i.e. material that will store and release heat). As building walls with large thermal mass can significantly contribute to the storage heat flux (Offerle et al. 2005a), this has the best summer daytime ΔQ S relation. This is evident for all three study areas (Fig. 6). When the paved fraction is high, the fraction of buildings and wall area is low and hence ΔQ S /(K ↓ + L ↓ ) decreases from the maxima of around 0.4 (paved fraction). Thus, buildings have a larger effect on ΔQ S compared to paved areas. Although Basel has the highest ΔQ S values (Fig. 4), London has higher overall ΔQ S when all available morning satellite overpasses are examined for 2016, thus, exemplifying the importance of meteorological conditions on ΔQ S . As expected, increased vegetation fractions (trees, grass) are linked to a decrease in ΔQ S (not shown) for all three study areas. Meteorological station (yellow dot). Cloud-masked areas (white). Note that the scales are different between maps

Concluding remarks
The SUEWS-ESTM scheme (available since version 2017c) is used to model ΔQ S in urban areas using EO data combined with ground-based meteorological forcing data and surface morphology, land cover and land use information.
Exploiting EO data to derive ΔQ S is challenging but the method presented has promise and allows the spatial variability of ΔQ S to be explored. The impervious surfaces (paved and buildings) contributes most to ΔQ S . Building wall area seems to explain variation of ΔQ S most consistently. Up to about 0.4 paved fraction, the increase is associated with a clear increase in ΔQ S /(K ↓ + L ↓ ); beyond this, ΔQ S /(K ↓ + L ↓ ) decreases. As areas with larger paved fraction, the fraction of buildings and wall area decreases, reducing the thermal mass required for high values of ΔQ S . The three cities have similar patterns between surface characteristics and ΔQ S /(K ↓ + L ↓ ). However, areas with higher urban density (e.g. central London) have larger fluxes as the greater building volume contributes to the ΔQ S term.
There are several challenges to estimating ΔQ S . Some issues are intrinsic to using EO satellites for T LST : the bias to clear sky conditions, and the momentary but infrequent nature of their sampling. The latter is critical given ΔQ S is a measure of the change in energy stored (or released) within the urban volume. We have resolved this by constructing a continuous T s dataset starting from the Lindberg et al.'s (2008) methodology. The original linear relation between maximum solar elevation and maximum (T a -T s ) for clear days is combined with diurnal sinusoidal variations in T s and clearness index (i.e. weather conditions) to adjust T s . Here, T LST is used to derive the T a and T s difference. Thus, as T a controls the change in both (T a and T s ), this may cause ΔQ S discrepancies, especially if T a variability is not accounted for. Improvements in surface temperature for different facets and their relation to different cooling/heating rates are being explored (Morrison et al. 2018(Morrison et al. , 2020. Other challenges are information received by the satellite sensor, i.e. what surfaces are seen from the sensor used to derive T LST . This is a well-known issue (Voogt 2008;Voogt and Oke 1997;Morrison et al. 2018Morrison et al. , 2020 and not considered in this study. Furthermore, the downscaling procedure can introduce biases in T LST (Mitraka et al. 2015). In addition, the accuracy and up-to-date status of the spatial information should also be considered. Although urban areas might seem relatively static, central London is undergoing constant urban densification . These factors will impact the estimated ΔQ S if the data used are not current. In the application here, material properties such as albedo, emissivity, volumetric heat capacity and thermal conductivity ( Table 2) do not vary with phenology and hydrology or other factors through the year. Yet, soil moisture will vary the soil thermal properties and LAI changes of vegetated surfaces modulate the intra-annual surface albedo. However, these effects are generally small due to the small contribution to ΔQ S from these land covers compared with built-up surfaces.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.