How to convert urban energy balance into contributions to urban excess temperatures?

Graphical abstract


Resource availability
The software implementing the Urban Heat Island decomposition has been written in R language (www.r-project.com) and the relevant R scripts are reported in this manuscript to allow the reproducibility of the results. ENVI-met is available from www.envi-met.com.

Method details
Attribution techniques are defined as methods " . . . evaluating the relative contributions of multiple causal factors to a change or event . . . " [1]. Typically they are used for decomposing climate change into its causes and determining the magnitude of the anthropogenic contribution. Such techniques can be utilised for different changes or events whereby one of them is the Urban Heat Island (UHI) effect that represents a temperature difference between urban and rural areas. The identification of dominant factors to UHI can improve the understanding of their dynamics and formation.
Often, UHI attribution focuses on the city scale [2,3] or only on the respective contributions of the surface energy balance (SEB) terms. The specific underlying processes responsible for that and their caused temperature change are not considered (e.g. [4]). At smaller scales, like building levels, the urban geometry becomes very complex, which confines such investigations to single objects or street canyons with predefined orientations. The focus on the identification of dominant processes again lacks of a statement to the caused temperature change (e.g. [5]).
In the present publication, we deduce a method for decomposing the Urban Heat Island intensity (DT) at the neighbourhood scale. The procedure is guided by a previously published approach [6]. In addition we share and describe an R [9] code especially developed for the decomposition scheme. The code is part of the work flow shown in Fig. 1 and involves three main steps that are explained in the subparagraphs below.
Step 1: ENVI-met simulation and reading binary files The micrometeorological model ENVI-met (version 3.1, [7]) was applied to simulate all necessary meteorological field variables in the neighbourhood of interest. The presented code can be easily extended and adapted to the needs of the user. Although the formation of an UHI effect requires low wind speeds, it has to be noted that the initial wind speed should be not smaller than 3 m/s. Otherwise the model can crash because of emerging numerical instabilities. The specific humidity at 2500 m (top of model atmosphere) was approximately derived from radio soundings. The choice of the lateral boundary conditions depends on the surroundings of the study area. A comparison with measurements showed that the combination of open (for temperature and specific humidity)forced (for turbulent kinetic energy) boundary conditions was best for our purposes. Anyway, it is highly recommendable to test different combinations in order to achieve the most accurate simulation results. Further details about model specifications and configurations for a specific case study are presented in [8] since these vary according to the local characteristics.

Reading ENVI-met output
The ENVI-met model output is binary and stored within the respective "*.EDT" files. Since the corresponding "*.EDI" files provide information about the data structure, we wrote an R-program extracting all meteorological fields. For any variable, the fields (matrices) of each height level (z = 0 . . . 29) are arranged one below the other starting with z = 0. To handle this structure, our code is divided into three parts: First, the binary format needs to be read into R and converted into a proper format. Second, we loop until we reach the variable of interest (loop index can be derived from "*.EDI" file by counting the position) and then, in a second loop, we select the height level of interest. Third, the resulting matrix is stored as a text file.

R program definitions
The subsequent code exemplifies this procedure for the sky view factor. The first lines (Fig. 2) define the number of grid cells, paths to the ENVI-met output folders and the location where the converted files should be stored as well as further parameters for reading the data.

Extracting required variables
The following lines (Fig. 3) convert the binary files into matrices. The first for-loop refers to the number of ENVI-met output files. The second one selects the sky view factor (here position 25) and the third one overwrites all previous height levels for the variables that are not used. In result, we receive the sky view factor at z = 0 (that is 0.5 m above surface) in a vector format. For later calculations, it is converted into a matrix.
This simple code can be modified and extended according to the variables needed or written as a function with arguments (e.g. "nrow", "ncol"; area size).
Step 2: DT decomposition approach

Theory
Pursuant to the formulation of [6], we interpret the decomposition approach for an urban neighbourhood. An urban neighbourhood results from a land use change compared to rural structures. This change is considered as a perturbation to the rural state. In this way, DT represents the UHI intensity in the respective neighbourhood. In order to quantify different heat flux contributions that are caused by the perturbation, we formulate the energy balance of an arbitrary surface in terms of radiation balance, turbulent heat fluxes, heat storage, and the anthropogenic induced heat flux [8].
Together with both sensible heat flux and latent heat flux Fig. 3. Converting variables of interest from binary to numeric matrix format.
this yields the surface energy balance (without advection) where R n is the net radiation balance, S # the incoming short-wave radiation, L # the incoming longwave radiation, Q S the storage heat flux, a the albedo, e the surface emissivity, s the Stefan -Boltzmann constant, r the air density, c p the specific heat of air at constant pressure (1013.25 hPa), T a the air temperature at a reference height, r a the aerodynamic resistance to heat diffusion and b ¼ Q H =Q LE the Bowen ratio. Positive fluxes are directed away from the surface, negative ones are directed towards the surface. For two sites with different surface characteristics, we can assume that at a certain height (reference height) the influence of the surface on the air temperature vanishes. In this case T a is equal for both sites and we linearise the emitted long-wave radiation in Eq. (3) at this point. After rearranging and factorising of T À T a we obtain Isolating T À T a finally yields Eq. (5) translates the surface energy balance into temperature differences [8]. If we consider a particular state as reference ("rural state") then the other one represents a perturbation ("urban state") caused by energy redistribution (right hand side of Eq. (5)). f is the fraction of energy redistribution associated with turbulent processes, which can either damp or amplify the reaction at the reference site to the perturbation. For that reason, Gðf Þ ¼ 1= 1 þ f ð Þ indicates the gain resulting from such turbulent redistribution processes. The gain function G (Fig. 4) tends to amplify (G > 1) for À1 < f < 0 b < 0 ð Þ and to damp (0 < G < 1) for f > 0. A negative gain is physically impossible and f has to be greater than À1. l 0 translates an energy flux to a temperature response (in Eq. (5)) and is called climate sensitivity.
In an urban environment, we consider the difference between urban and rural state. Applying this to Eq. (5) yields The temperature of the urban state T u can be replaced by the temperature of the rural state plus a perturbation caused by the urban structures (T r þ DT). Analogous replacements apply for f u ,R n;u , Q S;u and Q AH;u . D indicates small perturbations. Inserting these replacements into Eq. (6) and calculating the derivatives results in the UHI intensity (DT) of an urban neighbourhood Neglecting higher order terms O D 2 , it follows Eq. (9) is an approximation of the UHI intensity. All parameters on the right hand side of Eq. (9) that are associated with D characterising changes in the energy balance between a rural area and an urban neighbourhood. All other quantities, including l 0 = 1 þ f ð Þ, convert the changes in the respective energy term into a temperature change whose sum result in the UHI intensity for each model cell and time step in the quarter.

R program definitions
The first lines of this code (Fig. 5) define needed libraries, information about the locations where the results should be stored as well as the folder where the extracted ENVI-met variables are located.
Start of time loop and reading of the data Next (Fig. 6), the time loop over all available files starts (loop counter equal to the number of ENVImet output files). After constructing the complete path names of all input variables they have to be converted into a matrix since the chosen visualisation function needs this kind of class type. Finally, to complete the definition and data reading part, some parameters for DT-decomposition are defined.

DT decomposition: radiation
After these leading information, the first term in Eq. (9), characterising the change in radiation balance between the urban and rural state, is calculated (Fig. 7). As the surface emissivity is not available from the ENVI-met output we used the long-wave radiation budget ("lwbudget_sfc_rural") and the incoming long-wave radiation ("lwenvironment_sfc_rural + lwvegetation_sfc_rural") to derive this value for the rural state. The sensible ("H" (urban), "H_rural" (rural)) and latent heat fluxes ("LE" (urban), "LE_rural" (rural)) as well as the Bowen ratios ("beta_urban" and "beta_rural") have to be filtered to avoid critical cases (further discussed in [8]) within the calculation.

DT decomposition: convection and evapotranspiration
The second (change in convection efficiency) and third (change in evapotranspiration) terms in Eq. (9) are determined by the following lines (Fig. 8) DT decomposition: Heat storage and total DT Here, the anthropogenic heat flux (Q AH ) is neglected but can easily be implemented in the code. The effect of the change in heat storage ("Q_SH_term") is calculated as a residual so that, together with the other 3 terms, the UHI intensity ("DeltaT") for a neighbourhood is formed (Fig. 9). The variable "sealed_filter" and the subsequent 4 lines replace all sealed surfaces for the evapotranspiration term with zero values since these are impervious to water. Step 3: visualisation of DT contributions Using the program code described in step 1, all variables required for Eq. (9) are extracted. In step 2 the DT contributions are calculated. In step 3, applying the "image.plot" function, we exemplify the  spatial visualisation for the contribution of the radiative energy (Fig. 10). This can be extended for the other terms of the decomposition procedure or further analysis.