Inverse modeling of the urban energy system using hourly electricity demand and weather measurements, Part 2: Gray-box model

Accepted Manuscript Title: Inverse modeling of the urban energy system using hourly electricity demand and weather measurements, Part 2: Gray-box model Author: Afshin Afshari Nengbao Liu PII: DOI: Reference:

S0378-7788(17)30197-4 http://dx.doi.org/doi:10.1016/j.enbuild.2017.01.052 ENB 7327

To appear in:

ENB

Received date: Revised date: Accepted date:

22-10-2016 12-12-2016 17-1-2017

Please cite this article as: Afshin Afshari, Nengbao Liu, Inverse modeling of the urban energy system using hourly electricity demand and weather measurements, Part 2: Gray-box model, (2017), http://dx.doi.org/10.1016/j.enbuild.2017.01.052 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights

te

d

M

an

us

cr

ip t

Analysis of the energetic signature of a city Inverse gray-box dynamic model inferred from load and weather data Third-order lumped-parameter thermal network model of the urban system Offers insights into thermal properties of urban structures and microclimate UHI intensity is close to 2 ˚C while the corresponding cooling penalty is 12%

Ac ce p

    

Page 1 of 32

*Manuscript Click here to view linked References

Inverse modeling of the urban energy system using hourly electricity demand and weather measurements, Part 2: Gray-box model Afshin Afsharia,∗, Nengbao Liub Institute of Science and Technology, Abu Dhabi, UAE b University of Chicago, Illinois, USA

ip t

a Masdar

cr

Abstract

an

us

The interpretation of the energetic signature of a city, taken as an integrated urban energy system, has been attempted in this study. Instead of the customary bottom-up approach, we successfully identified an average gray-box dynamic model by inversion of measured hourly system load and weather data. The third-order lumped-parameter thermal network model is derived by applying energy conservation principles to major nodes representing buildings, paved surfaces and urban canopy. This powerful data-driven approach offers deep insights into the thermo-physical properties of the built structures as well as the coupling of buildings and the urban microclimate. In particular, the Urban Heat Island (UHI) intensity has been determined to be close to 2 ◦ C while the corresponding electricity cooling load penalty is 12%.

M

Keywords: Urban canopy model, Urban heat island, Lumped thermal parameter, Thermal network model, Gray-box model, Inverse problem, System identification

d

Nomenclature Specific enthalpy (MWh/Gg)

γ

Surface azimuth angle (radian)

γs

Solar azimuth angle (radian)

ω

Humidity ratio (kg/kg)

θ

Angle between solar zenith and wall orientation (radian)

ζ

Solar zenith angle (radian)

B

Direct beam solar irradiance (W/m2 or MW/km2 ): prescribed

D

Diffuse solar irradiance (W/m2 or MW/km2 ): prescribed

I

Global solar irradiance (W/m2 or MW/km2 ): prescribed

P

Cooling power (electricity) (MW): prescribed

p

Pressure (atm)

Parameters

Q

Thermal energy flow rate into a node (MW)

α

Ac

ce p

te

η

T

Temperature (◦C): prescribed or state variable depending on subscript

U

Wind speed (m/s): prescribed

X

Normalized node temperature (K)

y

Model output (MW)

Z

Height (m)

Constants κ

von Karman constant: 0.40

ρ

Density of air (Gg/km3 or kg/m3 ) :1.2 ×103

σ

Stefan-Boltzmann constant MW/km2 K 4 ): 5.67 × 10−8

L

Latent heat of vaporization of water: 694.5 MWh/Gg or 2265 J/kg

S

Specific heat of dry air: 0.2778 MWh/Gg.K or 1005 J/kg.K

(W/m2 K 4

or

Solar absorptance

∗ Corresponding

author Email addresses: [email protected] (Afshin Afshari), [email protected] (Nengbao Liu)

Preprint submitted to Elsevier

December 12, 2016

Page 2 of 32

Linear degradation slope of the chiller COP

f

Infiltration of urban air into the building



Emissivity

h

Horizontal

ψ

Glazing ratio (window to vertical facade)

i

Indoor air

ρ

Solar reflectance

p

Paved surface

A

Area (km2 )

r

Roof

C

Total node thermal capacity (MWh/K)

s

Sky

c

Thermal capacity (MWh/km3 K)

t

Anthropogenic heat due to traffic

F

View factor

u

Urban canyon air

f

Infiltration rate (air changes per hour, ACH)

v

Vertical

h

Time constant (h)

w

Wall

h

Total node thermal conductance (MW/Kkm2 or W/Km2 )

Subscripts - Directional heat flow

l

length of building (m)

R

Aerodynamic resistance (s/m)

w

Width of street (m)

Z0

Roughness layer height (m)

Zd

Zero displacement height (m)

Zh

Average building height (m)

Zr

Reference height (m)

c

Building cooling

1. Introduction

Indoor air node to sky node

iu

Indoor air node to urban air node (via glazing)

iw

Indoor air node to wall node

ps

Paved surface node to sky node

ua

Urban air node to atmospheric air node

up

Urban air node to paved surface node

us

Urban air node to sky node

wp

Wall node to paved surface node

ws

Wall node to sky node

wu

Wall node to urban air node

ce p

1.1. Context and problem statement

is

te

Atmospheric air

d

Subscripts - Nodes a

cr

volume

us

of

an

unit

M

per

ip t

β

Ac

In the UAE, about 60% of the total annual electricity consumption and up to 75% of the summer system peak demand can be attributed to buildings’ Air Conditioning (AC) load [1]. This remarkable fact is explained by the harshness of the climate, the relatively low performance of the buildings’ fabric/systems and the reinforcing interaction between the built environment and the urban microclimate. One of the most significant impacts of the latter is the Urban Heat Island (UHI). The quasi-totality of studies reported in the literature approach the UHI problem in 2 ways. Experimentalists either directly measure it and its consequences in-situ or develop reduced scale test facilities (e.g., wind tunnel). Proponents of numerical modeling build upon certain assumptions about the thermo-physical characteristics of the system (model parameters) based on prior information and then proceed to simulate current and future urban system state/outputin particular the energy consumption of buildingsusing known boundary conditions (system inputs). In the latter case, a limited experimental validation of the model is often conducted. The approach proposed here is radically different. We shall attempt to characterize and interpret the energetic signature of downtown Abu Dhabi, taken as an integrated urban energy system. For that, we define the aggregate urban energy process as a dynamic lumped-parameter input-output system. The ’rural’ meteorological boundary conditions are the input variables while the aggregate city-wide electricity load is the output. For the case at hand, we have a full year of hourly input/output data and would like to identify the model parameters from the available measurements. We use all the available data to directly infer the internal parameters of a gray-box model. To our knowledge this approach does not have a precedent in the field of urban micro-climate analysis. This problem belongs 2

Page 3 of 32

ip t

to the class of ’inverse’ problems. Inverse problems are often ill-posed and their numerical resolution is particularly challenging unless substantial prior information can be incorporated in the resolution procedure [2]. The advantage of the proposed procedure is that very little prior information about the urban system is needed–basically, the structure of the model, itself derived from the energy conservation principle, initial estimates of the model parameters and some light constraints–contrary to most prevailing ’forward’ urban modeling approaches. Once the parameters are identified, they inform on aggregate/average thermo-physical properties of buildings and the urban microclimate. In addition, the model enables the estimation of the UHI intensity and the relative importance of the different fluxes of energy. By varying the model parameters, It is possible to estimate the impact of different urban retrofit interventions such as upgrade of buildings’ envelope/AC, cool surfaces, etc., on aggregate electricity use and UHI intensity.

M

an

us

cr

1.2. Urban heat island Air temperature within urban canopies is, in average, higher than in the surrounding rural areas. This phenomenon is known as the Urban Heat Island (UHI) effect. It is usually more intense under cloudless sky and light wind conditions [3] that prevail in our case study. Several factors contribute to the UHI effect. They include (1) increased absorption of short-wave radiation due to the low solar reflectance of urban materials and urban geometry, (2) increased storage of sensible heat due to increased thermal diffusivity of urban materials, (3) anthropogenic heat production due, mainly, to AC waste heat rejection and motorized transportation, (4) reduced long-wave radiation losses, (5) lower evapotranspiration rates due to lack of vegetation and widespread use of impermeable surfaces, (6) lower sensible heat loss due to urban geometry that slows wind speed [4]. The UHI effect can significantly increase buildings’ energy demand for Air-Conditioning (AC) in warm countries [5]. The increase of building’s AC demand will not only cause more anthropogenic heat, but also intensify the UHI effect because of the lower efficiency of AC systems resulting from higher ambient temperatures. Some studies [6, 7] designate the above interaction as a ’vicious cycle’.

Ac

ce p

te

d

1.3. Simplified urban microclimate modeling Over the years, there has been a continuous refinement of simplified urban canopy modeling approaches 1 . Generally, such models, informed by substantial prior information about urban morphology and thermophysical properties, represent a ’forward’ relationship connecting input/boundary variables to state or output variables. Whenever measurements of state/output variables are available, they are partial and serve for validation or calibration of the model, not for parameter estimation. The most basic forward approach is the Slab model wherein the urban environment is modeled as a plate with a large roughness length and small albedo [8]. Slab models assume that buildings and roads are of the same temperature; buildings’ height and coverage ratio are not explicitly accounted for. Subsequently, more advanced urban canopy models were proposed. They separate the urban energy flows into budgets for roofs, roads and walls while taking the geometry explicitly into account to varying degrees–albeit in the form of an ideal regular and symmetric city [9, 8, 10, 11, 12]. For instance, coupled modeling of buildings and the urban canopy was attempted by Masson for a city composed of regular infinite urban canyons, but the so-called Town Energy Balance (TEB) model focused on the estimation of the urban conopy conditions rather than the buildings’ cooling load [9]. Further development includes multi-layer urban canyon models to resolve the vertical variation of the urban canopy layer [13], vegetation within the canyon [14] and AC waste heat rejection into the canyon [15]. The Table 1 below summarizes the main features of the existing simplified urban canopy energy models. 1.4. Modeling of buildings’ energy demand Usually, detailed characterization of one or several reference/baseline model(s) is essential for the estimation of the energy impact of building retrofits. Most reference models in the literature are derived using one of two approaches: top-down modeling or bottom-up modeling. Top-down models are further categorized as black-box or gray-box. The bottom-up approach relies on ’forward’ modeling of typical buildings. Hence, city models are little more than a collection of individual building representations, assuming rural (airport) climate. This is done using conventional 1 We

deliberately exclude detailed models based on Computational Fluid Dynamics (CFD) from this analysis.

3

Page 4 of 32

ip t cr

Table 1: Comparison of the main simplified urban energy balance model categories

Best [16]

Masson [9], Kusaka [8], Harman & Belcher [12]

Multi-layer

Krayenhoff & Voogt [17], Martilli et al. [13]

ce p

te

Single-layer

d

M

Slab

Description urban area represented as a flat slab, distinguished from rural areas by larger roughness length, smaller albedo, larger thermal capacity, anthropogenic heat release, etc represents a city as a layer of buildings with the overall surface heat exchange being the sum of exchanges on individual surfaces

Advantages

us

Reference

an

Category

Ac

similar to single layer but energy exchanges at multiple levels within the canopy

simplicity

more realistic representations of radiative trapping and turbulent exchange; computationally efficient

allows for varying building heights

Disadvatnages

not a realistic model of the urban canopy; in particular, fails in calculatng nocturnal cooling and mean wind/turbulence

city represented as an average oriented urban street canyon; urban canopy air conditions represented by a single value considers boundary conditions at several atmospheric levels; requires a higher discretization of the mesoscale model near the surface

4

Page 5 of 32

M

an

us

cr

ip t

building energy modelling and simulation software such as EnergyPlus [18]. A detailed physical description of the building typologies of interest is required. It can include the building’s construction material, geometry, geographic location, and the type of Heating, Ventilation and Air-Conditioning (HVAC) system. Such detailed forward models enable the evaluation of the energy performance of technological options and the prediction of the impact of changes in control strategies [19]. However, a large amount of prior information is needed, which makes the process of collecting this data time consuming and costly, especially at the urban scale. What’s more, the necessary calibration process is also a great challenge. Another weakness of this approach lies in its high dependence on the representativeness of the selected building typologies and the accuracy of the collected data [20]. Finally, the urban micro-climate and the associated heat island effect are typically not taken into account–with the possible exception of radiative exchanges [21]. The top-down approach results in an ’inverse’ input-output model–usually with lumped parameters [22] [23] [24]. The model structure has to be determined a priori so as to express the system output as a function of the system input. The model is ’gray-box’ if its structure has been determined in accordance to certain physical principles, such as the energy conservation principle. If, on the other hand, the model structure does not conform to any physical consideration and does not incorporate prior process knowledge, it will be called a ’black-box’ model. The model output can be linear or non-linear with respect to model parameters or with respect to the input. Additionally, the model can be stationary (algebraic equations) or dynamic (differential equations). Once the model structure is set, parameters are derived from the ’inversion’ of the measured inut-output data. In other words, a model identification algorithm is used to minimize an objective function that represents the distance between the measured system output and the modeled output generated by the model driven by the measured input data. In contrast to bottom-up models, the top-down model of urban energy systems usually considers the urban system in its totality (buildings + microclimate). The input/output data for this type of model tend to be general, such as aggregate electricity consumption data for a district or city, Gross Domestic Product (GDP), population growth, meteorological data, etc [25].

ce p

te

d

1.5. Incorporation of prior information in inverse models Black-box inverse models tend to have more parameters in order to achieve acceptable model accuracy (at least during the training stage) [26]. This type of model requires less engineering effort but depends heavily on data availability and quality. It is generally necessary to acquire data over a long period of time with widely varying conditions in order to train the model comprehensively. Otherwise, the model may lack robustness when applied to conditions that differ from those prevailing during the training stage (e.g., a change in control strategy) since the parameters do not necessarily respect the process physics [22]. Time series, linear/non-linear regression, and Artificial Neural Network (ANN) are good examples of black-box models. Gray-box inverse models are used in situations where the detailed system processes are not entirely known or where the entire system is too complicated to model in detail. Nonetheless, certain simple principles, such as the energy conservation principle, are applicable either to the whole system or to lumped-parameter nodes within the system. The physical knowledge about the process reduces the dimensionality of the parameter space. Furthermore, gray-box models tend to be more robust than black-box models [27].

Ac

1.6. Thermal network model of buildings Thermal energy flow in the urban canopy can be represented by a network of thermal resistances and capacitances, analog to an RC electrical network. The thermal network is able to dynamically model the underlying sensible and latent thermal energy flows between buildings and their urban environment. The basic idea of a thermal network model is to represent building elements as thermal elements in an electrical analog circuit and then use electricity circuit analysis techniques, such as Kirchhoff’s current and voltage laws, and nodal/mesh analysis [28, 29], in order to derive the governing differential equations. This approach generally has a simple linear structure (i.e., the output is linear with respect to inputs and state variables), while offering acceptable computational accuracy. The governing linear system of differential equations can, in the special case of a two-node model, be solved analytically [30, 29]. However, larger order systems generally require numerical resolution via standard temporal discretization techniques [31]. The discretized equations can be represented in a state-space formulation and solved efficiently by numerical packages such as Matlab. Alternatively, a network modeling platform such as Modelica or Matlab/Simscape can be used–in which case, certain non-linear modes of heat transfer such as long-wave radiation among surfaces and between surfaces and the sky need not be linearized. 5

Page 6 of 32

Studies of network thermal modeling of buildings and building zones/components abound. The distinction is usually in terms of number of resistors and capacitors: 2 or 3 resistors and 1 capacitor [32, 33], 3 resistors and 2 capacitors [30, 29], 5 resistors and 2 capacitors [34, 35], and so on. Table 2 presents a brief comparison of some of the main RC model types. Table 2: Comparison of the main lumped parameter RC model types

Reference

First order: 2R1C

[32, 33]

Second order: 3R2C

[30, 29]

Enhanced second order: 5R2C

[34, 35]

Description useful for certain types of simplified problem solving but is deficeint in two areas: (1) heat balance of a building space requires surface temperatures of all construction elements (not just one uniform lumped element temperature); (2) inaccurate modelling of the energy split between conduction and surface convection, again due to the uniformity of material temperature inherent in the method (this is especially bothersome in high thermal capacity construction elements) construction element is treated as a second-order lumped capacitance giving two material body temperatures, one of which can be used as a surface temperature prediction; if parameterized correctly, the second-order element is nearly identical to a fully distributed element with equivalent overall thermal resistance and thermal capacitance Underwood [35] shows that this 5-parameter second-order model with 3 series resistances and 2 capacitances can be quite accurate when compared to a range of building element thermal responses generated using a detailed finite-difference numerical model

ce p

te

d

M

an

us

cr

ip t

Types

Ac

K¨ampf et al. [36] started their model with a standard two-node lumped parameter model and then extended the model for a building with an arbitrary number of zones, accounting for the capacitance of internal walls separating the zones. Their study suggests several improvements to the standard handling of external surface radiation, external and internal convective exchanges and the placement of capacitance in multi-layered walls. Despite the great insight that it offers into the physical processes at work in a building, the loss of accuracy due to the simplified and linearized model structure is undeniable [37]. Nonetheless, lumped parameter models are quite appropriate as the basis for decision support tools that are aimed at optimizing energy efficiency or UHI mitigation interventions in situations where the baseline system can only be characterized with input-output time series measurements. 2. Method

The modeling approach adopted in this study consists in first separating the weather-driven air-conditioning load from the aggregate load, then trying to identify a dynamic state-space model to represent the purely weather-driven cooling load. The state-space model is derived from a 3-node lumped parameter thermal network representation of urban energy flows between buildings, urban canopy air, streets and the first level of the atmosphere. State variables are the non-observed temperatures at the different nodes of the model. Given the specificity of the urban area under study, the soil is assumed to be dry and urban vegetation is not considered. The humidity of the urban canopy air is assumed to be the same as its boundary value (first level of the atmosphere). 6

Page 7 of 32

2.1. Modeling stages

M

an

us

cr

ip t

Expanding on [1], a linear regression model of electricity consumption for the city of Abu Dhabi, UAE, was built, which enables efficient and accurate separation of the non weather-dependent cooling load from the hourly aggregate electricity consumption [38]. Different seasonalities (daily and annual), trend (due mainly to population growth), as well as marginal effect of weather variables, such as temperature, relative humidity, wind speed, and solar irradiance are taken into account in constructing this regression model. Afterwards, an urban canopy and building energy model based on a network of thermal resistances and capacitances is derived from first principles, with the previously determined weather-dependent cooling load as its output. The lumped parameter model represents the underlying aggregate energy interactions between buildings and their urban environment. The main advantage of this dynamic gray-box model is that its parameters have physical meaning. Consequently, by varying the parameters of the baseline model representing the current situation, it is possible to determine the impact, at the aggregate city level, of different energy efficiency and UHI mitigation interventions on curbing cooling demand and alleviating UHI intensity. For example, using parameters representing envelope heat conductance or albedo, it is possible to estimate energy savings that can result from city-scale programs aiming at enhancing envelope insulation or ’cool’ surfaces. Another extremely significant advantage of the model is its simplicity and computational efficiency. After constructing the model and formulating it as a space-space input-output system, the parameters were estimated using one year of input-output data (inverse identification). Given the inverse nature of the problem, the numerical estimation procedure is challenging if prior information about the process is not incorporated. The best parameter estimation results are obtained using genetic algorithm and particle swarm techniques. Still, unicity of the solution remains a problem. The problem is addressed by incorporating a priori information: fixing some of the parameters and enforcing min/max limits on the others. Finally, we estimate the UHI intensity, the cooling penalty due to the UHI, and the different energy fluxes in the urban canopy. The methods and conclusions laid out in this study can be applied to a wide range of energy modeling and urban studies.

d

2.2. Input-output data: Boundary conditions

Ac

ce p

te

The accuracy of the baseline model is dependent on the quality of the measured input-output data used to estimate model parameters. This study relies on substation-level hourly electricity data measured by the Supervisory Control and Data Acquisition (SCADA) system of Abu-Dhabi Emirate’s electricity utility as well as robust hourly weather data monitored by Masdar City’s weather station for the calendar year 2010 (January 1st - December 31st). The focus of this study is on buildings’ electrical load and, specifically, the cooling portion of it. Abu-Dhabi’s load is mostly residential and commercial. According to [39], as of 2010, the industrial & agricultural load represented less than 16% of the total load in the Emirate. The Emirate of Abu Dhabi covers a large area with multiple cities. A subset of substations which are all within the municipality of Abu-Dhabi is selected for this study. In aggregate, this subset constitutes a better proxy for the buildings’ load than the Emirate-wide system load. This subset includes all low-voltage substations of 11kV, 22kV and 33kV that directly supply street transformers serving final customers. The aggregate annual load profile of the selected substations for 2010 is shown in Figure 1. Annual electricity consumption of this subset is 6,228 GWh. Demand peaks at 1,096 MW. The base-load component of the load is approximately 40% of the peak. It should be noted that not all of the AC load is weather-dependent. Some AC components such as pumps and fans may be scheduled to operate year-round or during the entire cooling season (from May to October) without any regard for the actual weather conditions. In addition, some of the heat gains inside buildings (lighting, cooking, occupants’ metabolism) are driven by diurnal activity schedule rather than by weather. These result in socalled ’seasonal’ variations with daily, weekly and annual periodicities. These variations will be identified and set aside, together with the baseload and the trend, during the regression stage (see below). Hourly readings of atmospheric dry bulb temperature (T a ), relative humidity (RH), wind speed (Wa ), global horizontal solar irradiance (Ih ), direct normal solar irradiance (In ) and diffuse horizontal solar irradiance (Id ) were obtained from Masdar City’s weather station for 2010. Specific humidity ωa was calculated from relative humidity and dry bulb temperature using standard psychrometric formulas. The overall cross-correlation between some of the measurement streams was calculated and the results are shown in Table 3.

7

Page 8 of 32

ip t

1100

1000

cr

900

700

us

Load (MW)

800

600

an

500

400

0

1000

2000

3000

4000 5000 Time of year (h)

6000

7000

8000

M

300

Ac

ce p

te

d

Figure 1: Hourly load for 2010 of the urban area under study

Table 3: Weather data cross-correlation matrix

Ta ωa Ih Wa

Ta 1.00 0.48 0.48 0.50

ωa 0.48 1.00 -0.13 0.09

Ih 0.48 -0.13 1.00 0.38

Wa 0.50 0.0 9 0.38 1.00

8

Page 9 of 32

As could be expected, with the exception of wind speed, the weather variables are heavily correlated among each other and do not rigorously qualify as statistically independent. Furthermore, taken individually as time series, the weather variables display auto-correlation [40]. 2.3. Extraction of the cooling load from the aggregate load

Ac

ce p

te

d

M

an

us

cr

ip t

In order to define a baseline of cooling electricity consumption, a model accounting for the different variables affecting the electricity consumption is necessary. System-wide electrical load can be divided into two main constituents: the non-weather dependent load and the weather dependent load. The latter is entirely due to cooling of indoor spaces (a heating season is not apparent). The first constituent is mainly calendar-driven. It varies with hourof-day, day-type and season-without displaying any direct correlation with instantaneous variation of the weather variables. This non-weather dependent constituent of the load also includes the trend (growth). The goal of the pre-processing stage is to identify, via (mostly) linear regression, the non-weather dependent portion of the load. The remainder, the weather dependent portion of the load is then used in the next stage to estimate the parameters of the lumped-parameter network model. Trend In the last decade, Abu Dhabi has seen significant population growth. The energy consumption has increased accordingly. To account for this trend, two methods were tested: an additive linear trend and a multiplicative trend. The additive linear trend is based on the assumption that the rate of population change, new construction and new appliance acquisition result in a constant hour-over-hour increase in the electricity consumption throughout the year. This method considers that at every hour the load is increased by a constant amount from the previous hour, independently of other factors. In this case, the growth term is a linear parameter and can be estimated at the same time as the other model parameters after the data normalization process. The multiplicative trend is based on the assumption that the rate of growth affects all parameters in the model (in other words, all constituents of the load) and cannot be accounted for by a constant hour-over-hour increase of the load. The multiplicative coefficient multiplies all other coefficients, creating a model that is non-linear in the parameters. The multiplicative trend assumption is more consistent with the reality and produces more accurate results. It is therefore adopted despite the major nonlinearity (in the parameters) that it introduces. Seasonality Electricity demand is known to exhibit calendar driven fluctuations. These so-called ’seasonal’ fluctuations, not explainable by the weather variability, are caused by annual patterns, as well as higher frequency, weekly and daily fluctuations. The seasonality terms must have zero mean over their period of seasonality. In other words, in average, they should not introduce an offset to the baseload (the intercept of the regression equation). To account for the main seasonality terms-annual and daily/workday-Fourier series decomposition proved its superiority. Fourier series decomposition is based on the fact that a periodic function can be expressed by the sum of trigonometric functions, under certain conditions [40]. Considering a zero-mean seasonality function, a continuous Fourier series decomposition can be expressed as a summation of, possibly infinite, harmonics of sine and cosine waves. Over a typical year, the load follows a sinusoidal behavior, which justifies the use of a Fourier series representation of the annual seasonality term. The best annual seasonality representation for the lowest number of coefficients is achieved for a fourth-order Fourier series representation (8 parameters to estimate). Given the relative high frequency of the data (hourly), in addition to the annual sinusoidal variation, we can also detect intra-day sinusoidal variations. The regular workweek for most UAE office workers ranges from Sunday to Thursday, while Friday and Saturday constitute the weekend. In addition, certain major holidays affect the aggregate load profile and present similarities with the Friday profile, therefore public holidays are assimilated to Fridays. The month of Ramadan, a period of 29 or 30 days, also presents a specifically distinct diurnal profile. Therefore, a total of 4 day-types are to be distinguished: Workday (Sunday to Thursday), Friday/Holiday, Saturday and Ramadan. The best Workday daily seasonality representation for the lowest number of coefficients is achieved with a sixth-order Fourier series representation (12 parameters to estimate). 24 dummy variables (0/1) were used to account for special day types: Friday/Holiday, Saturday, Ramadan. Since we have a total of 3 special-day profiles (in addition to the normal Workday profile), there is a total of 72 hourly offset parameters to estimate.

9

Page 10 of 32

ip t cr us

an

Figure 2: Lumped parameter model expressing the energy balance between building and outdoor environment

2.4. Urban energy balance

ce p

te

d

M

Figure 2 displays the lumped parameter model expressing the energy interactions between the indoor and the outdoor environments. The atmospheric conditions are assumed to be fully prescribed since their determinants are monitored by our rural weather station (details below). The building’s indoor space is represented by a single thermal zone which is assumed to be at constant temperature (perfect AC control2 ). In the case presented in Figure 2, the building envelope, the urban canyon air and the outdoor paved surfaces are, each, represented by one lumped-parameter node. The model is therefore of third order since it has three free nodes3 . These three nodes exchange energy among themselves and with other boundary nodes that are either constant (the indoor temperature) or prescribed (external temperature, sky temperature). Solar radiation is equivalent to a current generator in an electrical circuit, and is also prescribed. In addition, wind speed and humidity variables–both prescribed–are used in the model, as will be shown below. Although this model is mainly focused on the sensible heat transfer flows, the outdoor and indoor humidity ratios are incorporated in the calculation of the infiltration/ventilation heat gains. The outdoor humidity is prescribed by the weather station measurements while the indoor humidity is assumed to be constant. The different modes of energy exchange are:

Ac

• Atmosphere and urban canyon: In presence of UHI effect, the temperature inside the urban canyon is different from (generally higher than) that in the first atmospheric level immediately above the canopy level. Energy exchange exists on the urban canyon canopy boundary with the atmosphere, and the exchange velocity is the main driver of this energy exchange process. • Urban canyon and paved surfaces: Our model assumes that paved surfaces represent all horizontal urban surfaces that do not belong to buildings. Indeed, in our domain of study, urban vegetation and water surfaces are negligible. Paved surfaces are in direct contact with the urban canyon and exchange energy with it via convection. In addition, the traffic-related anthropogenic heat (heat generated by motorized vehicles) will be released in the urban canyon. 2 On

rare occasions, when the outdoor temperature drops below the setpoint, this may result in a need for heating–or negative cooling possible implementation would prescribe separate nodes for walls and roof. However, our experience shows that the thermal state of the wall and roof nodes are not distinct enough to be estimated separately in an inverse procedure. 3 Another

10

Page 11 of 32

• Urban canyon and buildings: The main mode of energy exchange between urban canopy and buildings is convection at buildings’ surfaces (walls, windows, and roofs). Also, some of the waste heat from the AC system is released into the urban canyon. Finally, infiltration/ventilation is an important mode of energy exchange (both sensible and latent) between buildings and the urban canopy. • Solar (shortwave) radiation: Solar radiation is a special heat source as it heats all urban surfaces and even the indoor environment directly, through windows.

ip t

• Longwave radiation: Longwave radiation exists among different surfaces, such as road, wall, roof and sky. Heat loss to the sky can be high especially at night.

cr

2.5. Urban thermal network model

te

d

M

an

us

Heat storage capacity and heat conductance are two essential concepts underlying the thermal network representation of the interactions between buildings and urban climate. The structure of buildings, which includes walls, roof, floors and so on can store heat and also can transfer heat from outside to inside and vise versa. A convenient representation of the thermal network is an electrical circuit analog where temperature plays the role of potential difference and heat flow plays the role of current. Figure 2 shows the model used in this study. Air inside the building, building envelope, urban air which is inside the urban canopy, urban paved surfaces, sky and atmospheric air which is above the urban canopy are modeled as lumps. In this model, indoor air temperature, which is assumed to be constant throughout the year, ambient air temperature above the urban canopy, and sky temperature are prescribed boundary variables. All unknown temperature nodes are associated with capacitors. Solar radiation is modeled as a heat flux incident on the paved surfaces, building walls/roof, and the indoor space (through windows). Furthermore, longwave radiation between urban surfaces and the sky, as well as among surfaces, is modeled. In parallel to the sensible thermal network model, a relationship between enthalpy of air inside and outside buildings is derived in order to account for heat gains due to ventilation and infiltration. The advantage of the proposed thermal network model is that the parameters have physical significance. This means that the model can be easily adjusted to account for retrofits and UHI mitigation interventions. For example, if internal insulation is added to the building envelope, the model can be modified by adjusting the parameter corresponding to the indoor envelope conductance/resistance while keeping other parameters unchanged. Assumptions underlying our thermal network model include:

ce p

• Longwave radiation heat exchange between wall, roof and road are considered–albeit linearized. • Longwave radiant heat exchanges between sky and wall, road, roof, indoor environment are taken into account– albeit linearized. • Constant indoor zone temperature, resistances/conductances and capacitances.

Ac

• The Coefficient Of Performance (COP) of the AC system is dependent on the outdoor temperature and decreases as the temperature rises. The baseline COP is fixed while the degradation coefficient β (slope) is estimated. • Air inside the canyon is well mixed and homogeneous. • Wind velocity inside the canyon is homogeneous. 3. Problem formulation 3.1. Energy balance As the RC model is derived from energy conservation principles, for each dynamic node (connected to a capacitance), the rate of change of its internal energy is equal to the algebraic sum of the heat fluxes reaching it. Therefore, in our case, we obtain four differential equations expressing the sensible heat transfer process and another one which expresses the relationship between the indoor and the outdoor (urban canopy) enthalpies. 11

Page 12 of 32

For each node, a differential equation represents the heat transfer process. The general formulation of the heat transfer for a given node can be written: Cj

X dT j X = h jk (T k − T j ) + Qj dt j k

(1)

us

cr

ip t

where C j is the capacitance of node j and h jk is the heat transfer coefficient (or thermal conductance) between node j and node k. h jk is the inverse of r jk , the thermal resistance between node j and node k. T j is the temperature of node j while the T k are the temperatures of adjacent node(s) k. Q j is the thermal energy directly injected into node j (positive for heating, negative for cooling). Nodes include envelope (wall/roof), urban canyon air, paved surfaces and indoor air. We assume that unpaved outdoor surfaces have negligible area. Furthermore, the outdoor greenery and water surfaces are not accounted for in the model. These are acceptable assumptions, given the case at hand. The thermal state of the indoor zone is assumed to be perfectly controlled so that the indoor temperature and humidity are constant: T i = T i ∗ , ωi = ωi ∗ , ηi = ηi ∗ . Hence, the left side of the differential Equation 1 will be zero in the case of the indoor temperature T i . Furthermore, a relationship between indoor air enthalpy and outdoor air enthalpy is needed in order to model heat gains due to ventilation and infiltration: (2)

an

Q f = ρ f [L(ηa − ηi ∗ ) + S (T u − T a )]

M

where L is the latent heat of water in J/kg, S is the specific heat of air in J/Kkg, ρ is the density of air in kg/m3 and f is the infiltration/ventilation rate in m3 /s. In this study, we assume that ωi is perfectly controlled and equal to ωi ∗ = 0.008g/kg. The controlled indoor temperature, T i ∗ , on the other hand is unknown and will be estimated as part of the parameter identification procedure. Furthermore, the urban canopy air’s specific humidity ωu is assumed to be equal to that of the atmosphere ωa . This is very nearly true as demonstrated by our simultaneous measurements of urban and rural humidity. We can now write the main state equations:

d

0 = −hiw (T i − T w ) − his (T i − T s ) − hiu (T i − T u ) + αi Ai Iv − eP sen + Qi Cw dTdtw = hiw (T i − T w ) − hwu (T w − T u ) − hwp (T w − T p ) − hws (T w − T s ) + αw Aw Iv + αr Ar Ih u Cu dT dt = hwu (T w − T u ) + hiu (T i − T u ) − hua (T u − T a ) − hup (T u − T p ) + a(1 + e)P + Qt dT p C p dt = hup (T u − T p ) + hwp (T w − T p ) − h ps (T p − T s ) + hep (T e − T p ) + α p A p Ih

te

              

(3)

Ac

ce p

where subscripts i, w, u, p, s, a of conductance and capacitance parameters correspond to indoor air, wall/roof (envelope), urban canyon, paved surfaces (road), sky and atmosphere. αi , αw , α p are heat gain coefficient–solar absorptance (one minus albedo) for opaque surfaces, solar heat gain coefficient (SHGC) for glazed surfaces–of windows, wall/roof and paved surfaces. Ai , Aw , A p are the total areas of windows, wall/roof and paved surfaces. Ih and Iv are the global horizontal irradiance (GHI) and the global vertical irradiance (GVI). qa is the anthropogenic heat gain rate, e is the Coefficient of Performance (COP) of the AC systems, a is the fraction of waste heat going into the urban canyon, f is the ventilation and infiltration air flow rate, ηu , ηi , ηa are the enthalpy values of inside-canopy air, indoor air, and atmosphere, S a is the sensible heat of air. Internal gain variable Qi is set to zero hereafter, since, being mainly schedule driven, its impact on the cooling load does not depend on weather and is extracted during the pre-processing stage. The cooling power needed to compensate for the internal gain is deleted during the pre-processing stage (regression), together with all the other constant or schedule/seasonality driven or trend related components of the load. The resulting ’cooling power’ P is the portion of the cooling electricity load that is purely weather driven4 . P is the sum of the sensible cooling power P sen and the latent cooling power Plat . P sen is the electricity consumption used to offset the external sensible gains, Plat is the electricity consumption used to offset sensible and latent gains due to ventilation and infiltration: eP˙ lat = Q f . From Equation 2 we can derive the following final equation, which relates the state (node temperatures) to the output of the system P: 4 Strictly speaking, this means that the urban canopy air temperature calculated by our model slightly underestimate the true UHI intensity because it does not account for the portion of the chiller waste heat that results from offsetting non-weather dependent gains

12

Page 13 of 32

P sen = P − ρ f [L(ηa − ηi ∗ ) + S (T u − T a )]

(4)

Replacing for P sen in the energy balance equation of T i and solving for y = P, gives us the output equation of the state-space system. We then replace P in the differential equation expressing the energy balance of T u to obtain the final expression of the state-space system.

ip t

3.2. State-space representation

an

us

cr

The above equations have an equivalent state-space representation, which can be solved efficiently using numerical simulation tools, such as Matlab. Before proceeding, we reduce the number of state variables (temperature nodes) from 4 to 3 by combining the roof and wall nodes. We retain only the node T w which now represents the state of whole opaque envelope (walls + roof). This decision was taken because it was recognized, after many attempts, that the data did not enable us to separate the parameters pertaining to the walls from those of the roof–despite the fact that each node receives a different solar irradiance, respectively GVI and GHI. To simplify the equations, we define a new normalized state vector: X = T − T i ∗ where T is the vector of all 3 system nodes. In addition, we define Xa = T a − T i ∗ , X s = T s − T i ∗ and ∆η = ηa − ηi ∗ . The general state-space formulation then becomes: ( dX(t) dt = A · X(t) + B · U(t) (5) y = C · X(t) + D · U(t)

  hwp /Cw  Xw    hup /Cu   Xu    −h p /C p X p } |{z}

d

hwu /Cw (g − hu + ae hiu )/Cu hup /C p {z A

te

  −hw /Cw dX  = (hwu + ae hiw )/Cu  dt hwp /C p |

Ac

ce p

  0  + (hua − ae S f )/Cu  0 |

h y = hiw /e |

M

where X(t) is a vector of state variables that correspond to the normalized temperatures of each node associated with a capacitance. U(t) is a vector of inputs (prescribed boundary conditions). A, B, C and D are coefficient matrices. Expanding the previous system of equations, we obtain:

(hiu + S f )/e {z C

X

hws /Cw ae his /Cu h ps /C p

0 ae /Cu 0 {z

αw Aw /Cw ae αi Ai /Cu 0

   Xa    X  0   s    f ∆η 1/Cu    (6)   I  0  v  }  Ih  Qa |{z}

αr Ar /Cw 0 α p A p /C p

B

U

  i Xw  h 0  Xu  + −S · f /e } Xp  | |{z} X

his /e

1/e {z D

αi Ai /e

0

   Xa   X   s    i  Xe  0  f ∆η }  Iv     Ih    Qa |{z}

(7)

U

where y = P is the (weather dependent) cooling power, ae = a·(1+1/e), h p = h ps +hup +hwp , hw = hwu +hws +hwp +hiw , hu = hwu + hup + hua + hiu . Let us define the following time constants: τw = Cw /hw , τu = Cu /hu and τ p = C p /h p .

13

Page 14 of 32

Equations 6 and 7 are in continuous time-domain. We choose the semi-explicit temporal discretization method due to its balance of stability and accuracy. We obtain: (8)

X(n + 1) = (2I − dt · A(n + 1))−1 [(2I + dt · A(n)) · X(n) + dt · (B(n) · U(n) + B(n + 1) · U(n + 1)]

(9)

ip t

X(n + 1) − X(n) 1 = · (A(n + 1) · X(n + 1) + B(n + 1) · U(n + 1) + A(n) · X(n) + B(n) · U(n)) dt 2

where dt is the sampling interval.

cr

4. Implementation 4.1. Morphological characterization of the site

M

an

us

The model is a simplified representation of a portion of downtown Abu Dhabi (UAE). The geometric and morphological parameters were derived from GIS data for the site. Our representation of the urban canopy relies on only 3 parameters: average building height Zh , plan area index λP and frontal area index λF . h is area-weighted average of all buildings in the studied region. λP is the ratio of the plan area of buildings over the total site area. Finally, λF is the ratio of average vertical area of roughness elements (buildings) opposing the urban wind flow over the total site area. In the general case it is dependent on the actual wind flow direction and even on height above the ground. In most studies it is approximated as 1/4th of the vertical faade area divided by the site area. Instead of attempting to model the full geometric complexity of the original site, we analyze an equivalent regularized arrangement of buildings and streets that is based on the above 3 morphological parameters. This approach is based on the assumption that the site geometry is regular enough to be adequately represented by an array of regularly positioned identical building with a square footprint l · l interlaced by streets having width w. l and w can be derived from the primary 3 parameters mentioned above using straightforward geometric properties: λP λF

te

d

l=h

! 1 w=l √ −1 λP

(10) (11)

ce p

For the downtown area considered in this study, Zh = 34.7m, λP = 0.270 and λF = 0.280. Therefore, the regularized dimensions are:   Z = 34.7m    h l = 33.5m (12)     w = 30.9m

Ac

Roughness height and zero-plane displacement height are estimated using the following simple formulas [41]: Z0 = 0.15Zh = 5.20m

(13)

Zd = 0.67Zh = 23.1m

(14)

14

Page 15 of 32

ip t cr

Figure 3: Solar angles

ce p

te

d

M

an

us

4.2. Incident solar radiation 4.2.1. Global vertical irradiance Several nodes in our RC model receive solar radiation: wall/roof, paved surfaces, windows. The solar radiation incident on horizontal surfaces is different from that incident on vertical surfaces (walls). While measurements of Global Horizontal Irradiance are usually readily available from weather monitoring stations, the Global Vertical Irradiance (GVI) is not measured since it would be different for each surface orientation. In order to derive the GVI, we need beam radiation on the vertical surface as well as diffuse radiation (including reflection from other surfaces). The initial diffuse radiation from the sky is composed of three parts: isotropic part, circumsolar diffuse part, and horizon brightening part [42]. Different models have been proposed to properly calculate the total solar irradiation on a tilted surface, and their differences largely lie in the way they treat the three components of the diffuse radiation. The isotropic diffuse model, derived by Liu and Jordan [43] includes three components for the radiation on the tilted surface: beam, isotropic diffuse, and solar radiation diffusely reflected from the ground. The circumsolar diffuse part and the horizon brightening part are set to zero as all diffuse is assumed to be isotropic. The isotropic model is easy to understand and calculate. More complete models have been developed in order to take into account the two missing constituents of the diffuse radiation in the isotropic model. Hay and Davies [44] only consider one missing constituent, the circumsolar part, and estimate its fraction. Reindl et al. [45] refine the model of Hay and Davies by adding the horizontal brightening term. Perez et al. [46, 47] conduct a more detailed analysis of three diffuse components. In this study, the isotropic sky model is adopted due to its simplicity and sufficient accuracy given the other sources of uncertainty present. Hence, the diffuse solar radiation on a vertical surface (tilt angle β = π/2) is given by [42]:  1 + cosβ 1 − cosβ 1  + Ih (1 − αg ) = Dh + (1 − αg )Ih (15) 2 2 2 where αg is the solar absorptance of the ground is taken to be the same as α p . The beam irradiance incident on a vertical surface can be calculated by using the following trigonometric relationship [48]:

Ac

Dv = Dh

  Bv = cosθz cosβ + sinθz sinβcos(γ s − γ) B

(16)

where B is the direct normal (beam) irradiance, θz is the zenith angle, γ s is the solar azimuth angle and γ is the surface azimuth angle. These angles are shown in Figure 3. In the domain of study most building walls are oriented in the following directions: gamma = [3π/4, −3π/4, −pi/4, pi/4]. We can calculate the average diurnal profile of the global solar irradiance incident on each wall as well as the mean GVI on all four walls for the full year. The average annual GVI peaks at about 300Wm−2 . Using a cylindrical representation of the building as suggested in [49] results in very similar results. 4.2.2. Shading Equation 16 refers to a situation where the vertical wall is not shaded by neighboring structures. This assumption can result in significant errors if applied to a dense urban environment. Unlike roofs (assumed to be unshaded in this 15

Page 16 of 32

study), the solar radiation incident on walls and paved surfaces must be adjusted for shadowing. Let us consider the case of an idealized regular two-dimensional urban geometry with infinitely long canyons and buildings. The building cross-section has width l and height Zh corresponding to average characteristics of the real urban domain. Following [9], let θ be the angle between the sun direction (zenith angle) and the canyon axis and define θ0 as the critical canyon orientation for which the road is no longer sunlit. Then:
w , 1) Zh tan(ζ)

(17)

ip t

θ0 = sin−1 min(

Bw = Bh


w 1 θ0
1 ( − ) + tan(ζ) (1 − cos(θ0 ) Zh 2 π π

cr

where ζ is the zenith angle. Integrating over all canyon orientations, we can derive the following expressions for wall and paved surfaces direct beam radiation: (18)


1
2θ0 2 Zh − tan ζ + tan(ζ) (1 − cos(θ0 ) (19) π πw π where Bh = Ih − Dh is the direct beam irradiance on an unobstructed horizontal surface (e.g., roof). Note that this is not the same thing as B (DNI). Formulas for diffuse irradiance on walls and paved surfaces, accounting for infinite reflections, are also provided in [9]. Reflections are assumed to be isotropic (diffuse). Let us first derive the initial diffuse irradiance incident on walls and paved surfaces: Dw = Dh · F w (20) D p = Dh · F p

M

an

us

B p = Bh

Now we express the first (diffuse) reflection from walls and paved surfaces: Rw = (1 − αw )(Bw + Dw ) R p = (1 − α p )(B p + D p )

(21)

d

Finally, the total irradiance on walls and paved surfaces can be derived:

ce p

Where Mw and Mw are given by:

te

Iw = (Bw + Dw ) + (1 − 2Fw )Mw + Fw M p I p = (B p + D p ) + (1 − F p )Mw

Mw = Mp =

Rw +Fw (1−αw )Rp 1−(1−2Fw )(1−αw )+(1−F p)Fw (1−α p )(1−αw ) R p +(1−F p )ρ p (Rw +Fw (1−αw )R p ) 1−(1−2Fw )(1−αw )+(1−F p )Fw (1−α p )(1−αw )

(22)

(23)

Ac

Figure 4 displays the average annual diurnal profile of the global irradiance incident on walls and on paved surfaces with and without obstructions for ’typical’ values of l = 33.5m, Zh = 34.7m and w = 30.9m obtained from GIS data for the domain of study. The corresponding wall and pavement sky view factors are Fws = 0.381, Frs = 1 and F ps = 0.276 (see derivation of Fws and F ps below). 4.3. Longwave radiation

Urban structures often exchanges significant amounts of energy with the sky via long wavelength (’longwave’) radiative transfer. Thus, a suitable sky temperature estimate is needed in order to model the radiative heat transfer between those surfaces and the sky.

16

Page 17 of 32

unshaed (rural) shaded (urban) 800

600

200

0

2

4

6

8

10

12 Hour of day

14

16

20

22

24

unshaed (rural) shaded (urban)

us

250 200 150 100

an

Global Vertical Irradiance (W/m2)

300

50

2

4

6

8

10

12 Hour of day

14

16

18

20

22

24

M

0

18

ip t

400

cr

Global Horizontal Irradiance (W/m2)

1000

Figure 4: Average annual diurnal profile of obstructed and unobstructed GVI/GHI

te

d

4.3.1. Sky temperature Different sky temperature models have been proposed in the literature. The one reported to be the most accurate[50] is the TRNSYS model, which accounts both for clear and cloudy sky conditions. The TRNSYS model is as follows [50]:

ce p

ε1 = 0.711 + 0.005T d p + 7.3 × 10−5 T d2p + 0.013cos(2π ε2 = ε1 + (1 − ε1 ) fcloud εcloud Ts =

T a ε1/4 2

t ) + 12 × 10−5 (pa − p0 ) 24

(24) (25) (26)

Ac

where ε1 is the emissivity of clear sky, ε2 is the emissivity of cloudy sky, T d p is the dew point temperature at ambient conditions (◦C), t is the hour of the day, pa is the atmospheric pressure in atm, p0 is the atmospheric pressure at sea level, fcloud is a value between 0 and 1 indicating the degree of cloud, εcloud is 0.8, T s is sky temperature K), and T a is the ambient temperature (K). In this study, pa is set to 1 as Abu Dhabi is only slightly above sea level, and fcloud is set to 0 as it is rarely cloudy throughout the year. 4.3.2. Linearized radiative heat transfer coefficient between wall, road and sky The linearized longwave radiation heat transfer coefficient between walls, road and sky is given in Equation 27 hi j = 4(1 − εcan )εi ε j σFi j T i j

3

(27)

where εcan is the urban canyon air emissivity, ε is the emissivity of surface i, j ={wall, road, sky(ε sky = 1)}; σ = 5.67 × 10−8 Wm−2 K −4 is the Stefan-Boltzmann constant; Fi j is the view factor between surface i and j; and T i j is the average temperature of surfaces i and j in K. 17

Page 18 of 32

In our case, to simplify the estimation process, we assume, based on observations, that the average envelope temperature T w is 30 ◦ C, while the average road temperature T p is 35 ◦ C. The average annual value of the sky temperature T s is about 15 ◦ C. In the literature, the values of εw and ε p are reported to be 0.90 (brick) and 0.95 [51]. The air emissivity (εu ) is calculated from the following expression as a function of the humidity content and the size of the space [52]: εu = 0.683(1 − exp(−1.17χ1/2 )) (28)

(29)

us

  F ps = [(Zh /w)2 + 1]1/2 − Zh /w    w F ws = 0.5 Zh (1 − F ps )     F = 0.5(1 − F ) pw ps

cr

ip t

(pa + bpb ) and pw and pa are the partial pressures of water and air (in atmsopheres), respectively. where χ = pw Le T300 urb Le is the mean beam length calculated as 3.6 times the ratio between the transversal area of the urban canyon and its 1/2 perimeter and b = 5( 300 . Tu ) The sky view factors for road and wall as well as the wall-road view factor can be computed using an idealized two-dimensional representation of the canyon [9]:

an

where F ps and Fws are the sky view factors for paved surfaces and walls, Zh is the height of buildings, w is the width of road. In our case, we estimate based on the GIS data of downtown Abu Dhabi, that the plan area weighted average of the building height is Zh = 34.7m while the average street width is w = 30.9m. Replacing in equation 29, we obtain: F ps = 0.381, Fws = 0.276 and F pw = 0.362. 4.4. Output variable: Hourly cooling load

ce p

te

d

M

The output variable used in our model is the hourly weather dependent cooling load extracted by regressing the total system load against independent variables including hourly weather data and calendar data [38]. Figure 5 plots the actual system load versus the one estimated via regression for the year 2010. Clearly, weather has a significant impact on the load since the summer peak is more than double the winter peak. The non-weather dependent portion of the load, composed of daily/weekly/annual seasonality and trend (growth) can be separated from the total load and discarded. The remainder is the weather dependent cooling load5 . As shown on Figure 5, it constitutes 43% of the total annual load and 62% of the total peak demand. The cooling load due to infiltration/ventilation (i.e., driven by enthalpy) is also plotted to show that it is far from negligible. The model achieves an RMSE of 16.85MW (1.54% of peak) and an adjusted R2 of 0.9932. The intercept (baseload) is 486MW and the multiplicative growth rate is 2.64%. All parameters have t-statistics with absolute values exceeding 20 [38]. 4.5. Input variables: Hourly weather

Ac

The weather related variables available to our model are: dry-bulb temperature, relative humidity, wind speed/direction, diffuse horizontal irradiance (DHI), direct normal irradiance (DNI) and global horizontal irradiance (GHI). Figure 6 shows the annual variation of temperature and humidity. The shape of the temperature curve resembles that of the cooling load, indicating a strong correlation. The peak temperature in summer can exceed 45 ◦ C and seldom goes below 10 ◦ C. Relative humidity oscillates around 60% throughout the year, and is usually higher at night. In order to maintain the indoor thermal comfort, AC systems that provide fresh air must often both cool down and dehumidify the air before delivering it to the indoor space. To achieve the desired humidity level, the air must be cooled down to the corresponding dewpoint temperature (around 12 − 13 ◦ C). A procedure for accurately estimating the infiltration/ventilation load Q f is provided in [53]. Here, we simply assume that Q f is equal to the product of the air flow rate by the indoor/outdoor enthalpy differential. The average annual wind speed is 2.82 m/s although it fluctuates greatly. Wind affects cooling load in two opposite ways, one is to enhance convective heat transfer, which decreases the cooling load, the other is to aggravate infiltration, which increases the cooling load. The cooling load does not show discernible sensitivity to wind direction. 5 As mentioned previously, certain schedule driven constituents of the cooling load would also be discarded since the algorithm has no way of distinguishing them from non-cooling loads. In other words, the weather dependent cooling load estimate resulting from the regression probably slightly under-estimates the actual buildings’ cooling load.

18

Page 19 of 32

Regression model (Intercept = 486.14 MW , Peak = 1096.3 MW , Total = 6227.8 GWh) 1200 Measured load Modeled load

1100 1000 900 800

600 500 400 0

1000

2000

3000

4000

5000

6000

7000

8000

cr

300

ip t

700

Cooling contribution to annual consumption & peak demand: 43% & 62%

us

700 Baseload Total cooling load Infiltration load

600 500

an

400 300 200

0

0

1000

2000

3000

M

100

4000

5000

6000

7000

8000

te

d

Figure 5: Cooling Load

50

90

ce p

45

100

80

40

70

35

60

50

Ac

30

40

25

30

20

20 15 10

10

0 0

2000

4000

6000

8000

0

2000

4000

6000

8000

Figure 6: Temperature and relative humidity

19

Page 20 of 32

5. Initialization of the model 5.1. Morphology

Table 4: Configuration of Three Building Types

15 3.5 52.5 30 50 248

Retail 3 3.5 10.5 100 300 424

cr

Number of Floors Height per Floor (m) Height (m) Width (m) Length (m) Energy Intensity per year(kWhm−2 )

Office and Institutional 15 3.5 52.5 40 40 336

us

Residential

ip t

The domain of study is downtown Abu Dhabi. In this area, there are mainly three building types [54]: residential buildings, office/institutional buildings and retail buildings. The configuration of the average building in each of these categories is shown in Table 4.

M

an

From [55], we find that there are approximately 208, 000 residential units in Abu Dhabi. Assuming that the average area of each unit is 100m2 , then the total floor area of residential building is 20.8km2 . The total floor areas for office, institutional and retail buildings are, respectively, 2.9km2 , 5.85km2 and 1.78km2 [56]. Table 5 summarizes the total floor area of each building type and their respective fraction. Table 5: Reported Floor Area of Each Building Type

Residential Buildings 20.8 0.66

d

Reported Total Floor Area (Km2 ) Fraction

Offices and Institutional Buildings 8.75 0.28

Retails 1.78 0.06

ce p

te

The floor area above is for the Emirate of Abu Dhabi, while our study is restricted to the downtown area of Abu Dhabi city. Assuming that the downtown area has the same fraction of total floor area for each building type, then the total floor area of all buildings can be estimated using the ratios of aggregate energy consumptions. According to the data at our disposal, the total electricity consumption for Abu Dhabi downtown was 6228GWh in 2010, and the energy intensity for each building type can be found in Table 4. Then the estimated total floor area for all buildings in Abu Dhabi downtown is given in the equation below:

Ac

6228 = 22.0km2 (30) 0.66 · 248 + 0.28 · 336 + 0.06 · 424 Based on the assumption that Abu Dhabi downtown has the same fraction of each building type as reported in reports[55, 56], then we can calculate the total floor area for each building type, and further with the building configuration shown in table 4, we can estimate the number of buildings for each building type. Table 6 shows the number of buildings for each building type. Table 6: Estimated Number of Buildings for Each Building Type

Residential Buildings Estimated Total Floor Area (km2 ) Floor Area per Building (m2 ) Estimated Number of Buildings

14.6 22500 650

Offices and Institutional Buildings 6.16 24000 256

Retail 1.25 90000 14

20

Page 21 of 32

From the estimated number of buildings, we can further derive the total roof area, wall area and window area for each building type based on the building configuration shown in Table 4. Table 7 shows the estimated total area of roof, wall and window. Table 7: Estimated total roof, wall, window area

Estimated Total Roof Area (km2 ) Estimated Total Wall Area (km2 ) Estimated Total Window Area (km2 )

Retail

Sum

0.42 0.10 0.01

1.80 4.78 2.95

cr

0.98 3.41 2.05

Office and Institutional 0.41 1.26 0.89

ip t

Residential

us

We can further infer the total area of the ground (assumed to be entirely composed of paved surfaces) using the plan area index λP = 0.270. The total site area is 1.80/0.27 = 6.67km2 . The total floor area is 22km2 as mentioned before. For the model, we need to specify a ’typical’ building in an equivalent regularized urban area. The average building footprint in our domain is 666m2 . We can therefore assume that there are approximately 2702 ≈ 1.80/(666 · 10−6 ) ’typical’ buildings. We further assume that this building type has a glazing to (opaque) wall ratio of 0.5.

M

an

5.1.1. Initial values of the heat transfer coefficients Heat transfer coefficient between indoor zone and external wall The heat transfer coefficient between the indoor node and the external wall node combines two heat transfer modes in series: indoor convection and conduction through wall. Assuming, for simplicity, that the wall node is in the middle section of the wall, hiw , the heat transfer coefficient between the indoor zone and the wall node in Wm−2 K −1 can be expressed as in Equation 31 1 (31) hiw = 1 0.5 hconv,in + hww

Ac

ce p

te

d

where hconv,in is the indoor convective + radiative heat transfer coefficient of the wall and hww is the thermal conductance of the envelope exclusive of air film convective coefficients. In the literature, typical value of hconv,in and hww (for Abu Dhabi) are reported as 8.29Wm−2 K −1 [57] and 1.77Wm−2 K −1 [58], resulting in hiw = 2.48Wm−2 K −1 . Heat transfer coefficient between indoor zone and urban canyon The only direct link between the indoor node and the urban canyon air node is through the glazed surfaces. The corresponding heat transfer is the combination of three processes, in series: indoor convection, conduction through the window, and convection outside. Since the window U-value reported in 4 already includes so-called ’film’ conductances on both sides, we can simply set hiu equal to the windows’ U-value Uwin reported to be 3.88Wm−2 K −1 in average in downtown Abu Dhabi [58]. Heat transfer coefficient between wall and urban canyon The heat transfer between the wall node and the urban canyon air node combines two heat transfer modes in series, conduction through wall and outdoor convection: hwu =

0.5 hww

1 1 + hconv,out

(32)

hconv,out is assumed to be about 8.39Wm−2 K −1 which results from the application of Jurges’ formula [59] to the average urban canopy wind speed [60], resulting in hwu = 2.97Wm−2 K −1 . Heat transfer coefficient between paved surfaces and urban canyon The heat transfer between the road node and the urban canyon air node combines two heat transfer modes in series, conduction through the road mass and outdoor convection. A thermal conductivity of 1.00Wm−1 K −1 was reported for paved streets in [61]. Assuming a thickness of 0.5m, the resulting thermal conductance is h pp = 2.00Wm−2 K −1 . We can reasonably assume that the thermal mass of the road is much closer to its surface than to the bottom of the layer. If, for instance, we position the node at 0.05m below the surface (10 % of the depth), we can write: 21

Page 22 of 32

hup =

1 0.1 h pp

+

(33)

1 hconv,road

ip t

A value of 15Wm−2 K −1 has been reported in the literature for hconv,road , the heat transfer coefficient between the urban canyon air and the paved surfaces [3], resulting in hup = 8.57Wm−2 K −1 . Heat transfer coefficient between urban canyon and ambient air The thermal conductance coefficient between the urban canyon and the first level of the atmosphere is given as: hua = ρS Ch Ua = gua Ua

(34)

an

us

cr

where gua = ρS Ch , Ch is the bulk heat transfer coefficient, ρ is the density of air (1.2kg/m3 ), S is the sensible heat of air (1005J/kgK) and Ua is the wind speed at a reference height, assumed to be the height at which the first level of the atmosphere starts. According to [62], this height is approximately three times the average canopy height: Za = 3Zh . The bulk heat transfer coefficient Ch is related to Cd , the drag coefficient for momentum, according to: Ch = Cd /R. R was estimated by [63] to be 0.74. Following [64] and assuming neutral conditions (i.e., no buoyancy, which means that the Monin-Obukhov stability function [65] is negligible), we can estimate the drag coefficient, Cd , at height Z as in Equation 35: 2      κ   Cd =      ln ( Z−Zdu ) 

(35)

Z0u

te

d

M

where κ is the von Karman constant (κ = 0.40), Z0u is the urban roughness length for momentum and Zdu is the urban zero displacement height. As shown previously, for the downtown area considered in this study, Z0u = 5.20m and Zdu = 23.1m. We can now determine the drag coefficients for momentum and heat at reference height Za . According to Equation 35, Cd at height Z = Za is 0.0212 which means that Ch is approximately 0.0287. Replacing this value of Ch in Equation 34 gives: gua = 34.6J/Km3 . In order to estimate an initial value for hua we need the reference wind speed, i.e., the wind speed at the reference height. As suggested by Erell, the wind speed at reference height Za above the urban canopy can be estimated based on the assumption that the geostrophic wind speed is the same above neighboring rural and urban areas [66]:

ce p

Ua =

Z0u Z0r

!0.09

du ) log( ZaZ−Z 0u

log( ZZ0rm )

Um

(36)

Ac

Z0r is the roughness length of the rural site where the rural meteorological station is located (determined, via extensive wind speed gradient measurements, to be approximately 0.059m in our case), Z0u is the roughness length of the rural site, Zm is the height at which the wind speed is measured (in our case, Zm = 10m) and Um is the measured wind speed. Given that the average annual wind speed measured at our rural monitoring station at a height of 10m is about 2.82ms−1 , we can determine the average annual atmospheric wind speed at the reference height to be Ua = 2.26ms−1 . Hence, according to Equation 34, hua = 78.2WK −1 m−2 (per square meter of contact area between the top of the urban canopy and the first level of the atmosphere). Anthropogenic heat rate due to traffic As shown on 7, the average diurnal profile of traffic-related anthropogenic heat rate (sensible + latent) for the city of Abu Dhabi presents two distinct peaks, at 8 am and 9 pm. In both cases, the peak value is approximately 12 − 13W/m2 per m2 of total site area [67]. 6. Estimation results In order to estimate the model parameters, we minimized an objective function which coincides with the Root Mean Square Error (RMSE) of the model output. In other words, the cooling load calculated by the model was 22

Page 23 of 32

ip t cr us an

Figure 7: Normalized diurnal profile of anthropogenic heat gain in Abu Dhabi

te

d

M

compared to that resulting from the pre-processed system load data. The Particle Swarm and Genetic Algorithm optimization functions of Matlab were tested and the latter was retained. The calculation time was brought to reasonable proportions by encoding the RMSE calculation function in C. This function calculates the total RMSE over a full calendar year of 8760 hours (2010). Estimating all model parameters would result in an ill-posed problem characterized by the non-unicity of the solution. The analysis of the previous section enables us to derive ’best guess’ initial values as well as min/max thresholds for some of the model parameters. This is crucial in order to achieve convergence in this severely ill-posed problem.

ce p

6.1. Monte Carlo simulation Initially, we approach the problem by attempting to estimate all model parameters–with the exception of e0 , a and Qt which are unidentifiable and are set, respectively to 2.5, 0.50 and 12.6W/m2 . Multiple runs of the estimation algorithm reveal a lack of unicity of the solution which is due to the ill-posedness of the inverse problem. We conduct a Monte Carlo analysis by running, several hundred times, the parameter estimation algorithm until convergence. We record the UHI intensity ∆T u−r (◦C) and the UHI cooling penalty ∆Qu−r (%) resulting from each run. The histograms in figure 8 are proxies for the probability density distributions of these variables.

Ac

6.2. Fixed parameters In order to attempt to restore the unicity of the solution, some of the parameters are fixed during the estimation process. These parameters are selected because they fall under one or several of the categories described below: • Forcing them to zero does not create a bias in the model (for instance, the stack and temperature effect parameters f1 and f2 ) • They can be determined from prior information with an acceptable level of accuracy (for instance, window U-Value, surface albedos) • Their estimated value is stable and varies little during the Monte Carlo runs6 (for instance, the road time constant) 6 In

other word, the spread of the distribution/histogram of the parameter is ’small’.

23

Page 24 of 32

60

40

0

0

1

2

3 ∆Tu−r (°C)

4

5

6

us

60

40

20

10

15

20 25 Cooling load penalty (MW)

30

35

40

M

5

an

Number of runs per bin

80

0

ip t

20

cr

Number of runs per bin

80

Figure 8: Results of Monte Carlo Simulation

d

• They cannot be estimated from the data (AC waste heat fraction a, AC nominal COP e, traffic related anthropogenic heat rate qa )

te

These parameters and their values are given in Table 8.

Table 8: Fixed model parameters

Value 3.88 32 0 0 0.70 0.90 0.50 2.5 12.6

Unit W/Km2 h ACH/K ACH(m/s)−1

Ac

ce p

Parameter hiu τp f1 f2 αw αp a e0 Qa

Comment Glazing U-Value [58]

50% rooftop, 50% window/split Derived from experimental data cf. [60]

W/m2

6.3. Estimated parameters After running the estimation for about 1000 iterations, convergence was achieved (RMS Echange < 10−9 ) 9. The estimated values of the unknown model parameters are reported in Table 9. The RMSE is 16.99MW and the adjusted R2 is 0.9880. The estimated parameter values and UHI intensity are quite stable when the estimation procedure is repeated, but not strictly unique. Therefore the model is best suited to comparative analysis studies, for instance, for evaluating the relative impact on the urban system of changing some of the parameters. In order to be able to compare the urban and rural cases, we need to determine the cooling load of the same building as in the urban model but situated in a rural setting. A modified thermal network model, without the 24

Page 25 of 32

25

hiw kwu his hwp hps hus hws

20 15 10 5 0

0

100

200

300

400

500

600

0.3

0.1 0

0

100

200

300

400

500

800

900

600

700

800

900

an gua kup tw tu Ti

40 30 20

0

100

200

300

M

10 0

900

700

600

50

800

us

f0 ai gg

0.2

700

cr

0.4

ip t

30

400

500 Iterations

ce p

te

d

Figure 9: Parameter convergence after approximately 1000 iterations

Ac

Table 9: Estimated model parameters

Parameter hiw hwu gua hup hus his hws h ps hwp τw τu f αi β Ti

Value 4.17 6.47 36.4 20.4 24.7 0.02 2.33 2.97 8.42 36.4 11.7 0.185 0.0287 0.325 22.2 ◦C

Unit W/Km2 W/Km2 J/Km3 W/Km2 W/Km2 W/Km2 W/Km2 W/Km2 W/Km2 h h ACH

25

Page 26 of 32

urban air and road nodes, is therefore constructed. It is then run with the same building parameters as the converged urban model. Using the building cooling load determined in this way, we can calculate the so-called ’UHI cooling penalty’: The incremental cooling energy that is due to the urban microclimate. The properties of the estimated model are presented in Table 10. The UHI penalty for the cooling load (thermal) is close to 15% while the cooling electricity penalty is 12%7 .

Unit MW

◦

cr

GWh GWh

C

us

Value 16.99 0.9880 44370 38610 14.9% 12.0% 1.78

ip t

Table 10: Model properties

Parameter RMSE adjusted R2 Annual cooling load - urban Annual cooling load - rural Cooling load (thermal) penalty Cooling load (electricity) penalty UHI intensity

an

Tables 11, 12, 13 and 14 present the annual energy balance at each node. All constituents of the energy balance are expressed as percentages of the annual cooling load (note that the latter is different for the rural and urban cases). A positive value denotes net gain while a negative value denotes net loss. Table 11: Annual heat balance at the indoor node (heat gains expressed as a % of the cooling load)

rural 55% 21% 21% 3% 100%

M

te

d

envelope outdoor (glazing) outdoor (infiltration) solar (shortwave) sum

urban 54% 24% 20% 2% 100%

ce p

Table 12: Annual heat balance at the opaque envelope node (heat gains expressed as a % of the cooling load)

rural -54% -50% -36 % 140% 0

urban -54% -26% -24% 1% 103% 0

Ac

indoor outdoor sky (longwave) road (longwave) solar (shortwave) sum

Figure 10 displays the average diurnal variation of the UHI intensity ∆T u−r and the +/-σ (one standard deviation) confidence interval. Actual/modeled cooling load and the RMSE of the model are presented in Figure 11. 6.4. Discussion The estimated parameters inform on the average properties of the urban thermal processes in Abu Dhabi. In particular, the building related parameters characterize an average Abu Dhabi building. One of the main parameters of interest is the overall U-Value of the envelope hww including air film on both sides. We can determine the value of this parameter to be: 7 Note

that the determination of one of these values from the other is not a straightforward matter since the COP in our model is highly variable.

26

Page 27 of 32

Table 13: Annual heat balance at the road node (heat gains expressed as a % of the cooling load)

ip t

-78% -17% -1% 96% 0

cr

urban air sky (longwave) envelope (longwave) solar (shortwave) sum

Table 14: Annual heat balance at the urban air node (heat gains expressed as a % of the cooling load)

us

26% -24% -20% 78% -18% -124% 72% 10% 0

M

an

envelope indoor (glazing) indoor (infiltration) road atmosphere sky (longwave) anthropogenic (AC waste) anthropogenic (traffic) sum

te

d

15

ce p

10

∆Tu−r (°C)

5

Ac

0

−5

−10

2

4

6

8

10

12 Hour of day

14

16

18

20

22

24

Figure 10: Average diurnal profile of UHI intensity

27

Page 28 of 32

700 Measured urban Model urban Model rural

600

400

ip t

Cooling electricity (MW)

500

300

200

cr

100

−100

1000

2000

3000

4000

5000

1000

2000

3000

4000 Hour of year

50 0 −50

5000

6000

7000

8000

7000

8000

M

−100

6000

an

Residual (MW)

100

us

0

d

Figure 11: Cooling load (actual/model) and RMSE

1 hiw

1 +

te

hww =

1 hwu

= 2.54W/Km2

(37)

Ac

ce p

The average building time constant of 36.4h is quite plausible. So is the conductive + convective (mostly convective) heat transfer coefficient of the paved street surfaces estimated to be 20.4W/Km2 . The estimated infiltration/ventilation rate of 0.185ACH is somewhat lower than expected but this can be explained by the fact that most cooling is done by recirculation of indoor air and mechanical ventilation to bring in fresh air in a controlled manner is infrequent. The estimated rate of degradation of the nominal COP with outdoor temperature (0.325) is similar to what we observed during a pilot project where 10 buildings in downtown Abu Dhabi were monitored over the cooling season of 2014 (unpublished). The linearized longwave coefficients his , hws , h ps and h pw are, for the most part, well within the expected ranges. his is negligible. hws and h ps are between 2 and 3. On the other hand, hwp , at 8.42W/Km2 , is several times higher than the initial value calculated previously. The estimated value of gua is very close to the initial value determined previously. The estimated indoor temperature of 22.2◦C is consistent with the local preference for the AC thermostat setpoint (generally lower than the standard 24◦C). The UHI intensity is of the same order as that determined by a forward urban energy model, derived for the same city [60]. The cooling load penalty sue to the UHI effect is significant. It causes a 12% rise of the system-wide cooling electricity. In comparison to the total electricity consumption (cooling and non-cooling), the rise is about 5%. One fact that has often been overlooked in simplified urban energy models is the relatively significant thermal inertia of the urban canopy air node which influences the dynamic behavior of the urban system. The biggest surprise however is the low value of the solar gain coefficient for the indoor node. Given the generally high Solar Heat Gain Coefficient (SHGC) of the glazing in Abu Dhabi, a much higher value was expected. Table 11 shows that while about half of the building’s heat gains are from the opaque envelope, infiltration explains close to 20% of the total heat gains. The solar gains through the glazing are low. The highest gain on the building envelope node is of solar origin. Although the latter is smaller in the urban setting than in the rural setting, the 28

Page 29 of 32

ip t cr us

an

Figure 12: Annual energy flows in the thermal network as a percentage of the cooling load (urban case)

te

7. Conclusion and future work

d

M

longwave losses and convective losses are also severely curtailed oin the urban case. The different energy flows for the urban case are shown on Figure 12. As shown on Table 14, while the contribution of the traffic related anthropogenic heat to the thermal balance of the urban canopy air node is small, that of the AC waste heat is more than 7 times larger, making it equivalent to the convective gain from the paved surfaces. According to our model, the main mode of cooling for the urban air is the longwave radiative exchange with the sky.

Ac

ce p

This study demonstrates the possibility to identify, via an inverse resolution procedure, a simple urban energy model using hourly weather and system load data. The identified model can then be used to estimate the UHI intensity and the UHI cooling load penalty. In addition, because the model is based on the expression of a dynamic energy balance of the urban canopy, its lumped parameters have physical meaning and can be used to understand the properties of the aggregate urban system. The city of Abu Dhabi was taken as a case study, although this method is general enough to be applied to any other urban area. The best approach is to first run a Monte Carlo simulation in order to get an idea of the spread of the parameters and detect unstable ones. The estimation procedure can be regularized by fixing the unstable parameters and setting min/max thresholds on the others. Most identified parameters are within expected range, with the exception of the solar gain through the glazing which is very small. The average diurnal profile of UHI intensity in Abu Dhabi reaches a maximum of about 7 degrees at 6 am and a minimum of about -5 degrees at 1 pm. The average annual UHI intensity is probably between 1.5 and 2 degrees. The UHI cooling penalty (thermal) is not a linear function of the UHI intensity and sharply increases as the UHI becomes more severe. In our case, the increase stands at about 15% in comparison to the rural cooling load. If we consider the cooling electricity penalty (electricity used the AC equipment), the increase is 12%. Overall, therefore, mitigation of the UHI effect in Abu Dhabi’s can bring about a 5% decrease in the aggregate load of the grid. The impact on peak load is even more significant since most of the summer peak demand is due to air conditioning. Using this figure, it becomes possible to conduct a life cycle cost analysis of candidate interventions. Future work will include validation of the model using on site measurements of urban temperature (ongoing measurement campaign). We will also focus on improving the model and stabilizing the estimation procedure. For instance, alternative thermal network model structures or regularization methods can be tested. The aerodynamic parameterization of the model will be improved and the sensitivity of the model to different formulations will be tested. The surface convection parameters can be variable functions of the wind speed in the urban canopy rather than 29

Page 30 of 32

constant. The physical foundation of the model makes it possible to avoid physically implausible parameters. This ensures that the model not only fits the historical data but is robust in a prediction (scenario generation) application. This predictive capability will be used to conduct life cycle assessment of city scale UHI mitigation and energy efficiency scenarios. Finally, we will attempt to apply the procedure to other cities where the UHI is an issue. Acknowledgment

ip t

This research was supported by grant number 12MAZA5 of the Executive Affairs Authority (EAA) of Abu Dhabi, United Arab Emirates, as part of the Comprehensive Cooling Program.

cr

8. Bibliography

Ac

ce p

te

d

M

an

us

[1] L. Friedrich, P. Armstrong, A. Afshari, Mid-term forecasting of urban electricity load to isolate air-conditioning impact, Energy and Buildings 80 (2015) 72–80. [2] A. Tikhonov, V. Arsenin, Solutions of Ill-Posed Problems, Wiley, New York, 1977. [3] B. Bueno, L. Norford, G. Pigeon, R. Britter, A resistance-capacitance network model for the analysis of the interactions between the energy performance of buildings and the urban climate, Building and Environment 54 (2012) 116–125. [4] L. Gartland, Heat islands: understanding and mitigating heat in urban areas, Routledge, 2010. [5] M. Mohan, Y. Kikegawa, B. Gurjar, S. Bhati, A. Kandya, K. Ogawa, Urban heat island assessment for a tropical urban airshed in india, Earth Environmental Sciences. [6] Y. Ashie, V. Thanh Ca, T. Asaeda, Building canopy model for the analysis of urban climate, Journal of wind engineering and industrial aerodynamics 81 (1) (1999) 237–248. [7] Y. Kikegawa, Y. Genchi, H. Yoshikado, H. Kondo, Development of a numerical simulation system toward comprehensive assessments of urban warming countermeasures including their impacts upon the urban buildings’ energy-demands, Applied Energy 76 (4) (2003) 449–466. [8] H. Kusaka, H. Kondo, Y. Kikegawa, F. Kimura, A simple single-layer urban canopy model for atmospheric models: Comparison with multilayer and slab models, Boundary-Layer Meteorology 101 (3) (2001) 329–358. [9] V. Masson, A physically-based scheme for the urban energy budget in atmospheric models, Boundary-layer meteorology 94 (3) (2000) 357–397. [10] M. Thatcher, P. Hurley, Simulating australian urban climate in a mesoscale atmospheric numerical model, Boundary-Layer Meteorology 142 (2012) 149–175. [11] M. Kanda, T. Kawai, M. Kanega, R. Moriwaki, K. Narita, A. Hagishima, A simple energy balance model for regular building arrays, Boundary-Layer Meteorology 116 (3) (2005) 423–443. [12] I. Harman, S. Belcher, The surface energy balance and boundary layer over urban street canyons, Quarterly Journal of the Royal Meteorological Society 132 (621) (2006) 2749–2768. [13] A. Martilli, A. Clappier, M. Rotach, An urban surface exchange parameterisation for mesoscale models, Boundary-Layer Meteorology 104 (2) (2002) 261–304. [14] S.-H. Lee, S.-U. Park, A vegetated urban canopy model for meteorological and environmental modelling, Boundary-layer meteorology 126 (1) (2008) 73–102. [15] Y. Ohashi, Y. Genchi, H. Kondo, Y. Kikegawa, H. Yoshikado, Y. Hirano, Influence of air-conditioning waste heat on air temperature in tokyo during summer: numerical experiments using an urban canopy model coupled with a building energy model, Journal of Applied Meteorology and climatology 46 (1) (2007) 66–81. [16] M. Best, Representing urban areas within operational numerical weather prediction models, Bound.-Layer Meteorol. 114 (2005) 91–109. [17] E. Krayenhoff, J. Voogt, A microscale three-dimensional urban energy balance model for studying surface temperatures, Boundary-Layer Meteorology 123 (2007) 433–461. [18] D. Crawley, L. Lawrie, F. Winkelmann, W. Buhl, Y. Huang, Pedersen, S. C.O., R.K., R. Liesen, D. Fisher, Witte, J. M.J., Glazer, Energyplus: creating a new-generation building energy simulation program, Energy and Buildings 33 (2001) 319–331. [19] L. Swan, V. Ugursal, Modeling of end-use energy consumption in the residential sector: A review of modeling techniques, Renewable and Sustainable Energy Reviews 13 (8) (2009) 1819–1835. [20] L. Shorrock, J. Dunster, The physically-based model brehomes and its use in deriving scenarios for the energy use and carbon dioxide emissions of the uk housing stock, Energy Policy 25 (12) (1997) 1027–1037. [21] C. Reinhart, C. Davilla, Urban building energy modelling - a review of nascent field, Building and Environment 97 (2016) 196–202. [22] J. Braun, N. Chaturvedi, An inverse gray-box model for transient building load prediction, HVAC&R Research 8 (1) (2002) 73–99. [23] S. Wang, X. Xu, Parameter estimation of internal thermal mass of building dynamic models using genetic algorithm, Energy Conversion and Management 47 (13) (2006) 1927–1941. [24] W. Colea, K. Powella, E. Haleb, T. Edgar, Reduced-order residential home modeling for model predictive control, Energy and Buildings 74 (2014) 69–77. [25] A. Grandjean, J. Adnot, G. Binet, A review and an analysis of the residential electric load curve models, Renewable and Sustainable Energy Reviews 16 (9) (2012) 6539–6565. [26] U. Forssell, P. Lindskog, Combining semi-physical and neural network modeling: an example of its usefulness, in: Submitted to the 11th IFAC Symposium on System Identification (SYSID’97) to be held in Fukuoka, Japan, Citeseer, 1997. [27] H. Nielsen, H. Madsen, Modelling the heat consumption in district heating systems using a grey-box approach, Energy and Buildings 38 (1) (2006) 63–71.

30

Page 31 of 32

Ac

ce p

te

d

M

an

us

cr

ip t

[28] A. Tindale, Third-order lumped-parameter simulation method, Building Services Engineering Research and Technology 14 (3) (1993) 87–97. [29] T. Nielsen, Simple tool to evaluate energy demand and indoor environment in the early stages of building design, Solar Energy 78 (1) (2005) 73–83. [30] J. Crabb, N. Murdoch, J. Penman, A simplified thermal response model, Building Services Engineering Research and Technology 8 (1) (1987) 13–19. [31] C. Carter, Computational methods for passive solar simulation, Solar energy 45 (6) (1990) 379–384. [32] C. Lombard, E. Mathews, Efficient, steady state solution of a time variable rc network, for building thermal analysis, Building and Environment 27 (3) (1992) 279–287. [33] E. Mathews, P. Richards, A tool for predicting hourly air temperatures and sensibles energy loads in buildings at sketch design stage, Energy and buildings 14 (1) (1989) 61–80. [34] F. Lorenz, G. Masy, M´ethode d’´evaluation de l’´economie d’´energie apport´ee par l’intermittence de chauffage dans les bˆatiments, Traitement par differences finies d’un model a deux constantes de temps, Report No. GM820130-01. Faculte des Sciences Appliquees, University de Liege, Liege, Belgium. [35] C. Underwood, An improved lumped parameter method for building thermal modelling, Energy and Buildings 79 (2014) 191–201. [36] J. K¨ampf, D. Robinson, A simplified thermal model to support analysis of urban resource flows, Energy and buildings 39 (4) (2007) 445–453. [37] M. Gouda, S. Danaher, C. Underwood, Building thermal model reduction using nonlinear constrained optimization, Building and Environment 37 (12) (2002) 1255–1265. [38] A. Afshari, L. Friedrich, Urban energy analysis using aggregate electricity demand and weather measurements, part 1: Forward model, Energy and Buildings [submitted]. [39] Statistical Yearbook of Abu Dhabi, Statistics Centre, Abu Dhabi, 2013. [40] W. Wei, Time Series Analysis : Univariate and Multivariate Methods (2nd edition), Pearson, 2005. [41] R. Britter, S. Hanna, Flow and dispersion in urban areas, Annu. Rev. Fluid Mech 35 (2003) 469–496. [42] A. Duffie, A. Beckman, Solar engineering of thermal processes, 4th edition, Wiley, 2013. [43] B. Liu, R. Jordan, The long-term average performance of flat-plate solar-energy collectors: with design data for the us, its outlying possessions and canada, Solar Energy 7 (2) (1963) 53–74. [44] J. Hay, J. Davies, Calculation of the solar radiation incident on an inclined surface, in: Proc. of First Canadian Solar Radiation Data Workshop (Eds: J.E. Hay and T.K. Won), Ministry of Supply and Services Canada, Vol. 59, 1980. [45] D. Reindl, W. Beckman, J. Duffie, Evaluation of hourly tilted surface radiation models, Solar Energy 45 (1) (1990) 9–17. [46] R. Perez, R. Seals, P. Ineichen, R. Stewart, D. Menicucci, A new simplified version of the perez diffuse irradiance model for tilted surfaces, Solar energy 39 (3) (1987) 221–231. [47] R. Perez, R. Stewart, R. Seals, T. Guertin, The development and verification of the perez diffuse radiation model, Tech. rep., Sandia National Labs., Albuquerque, NM (USA); State Univ. of New York, Albany (USA). Atmospheric Sciences Research Center (1988). [48] J. Duffie, W. Beckman, Solar engineering of thermal processes, John Wiley & Sons, 2013. [49] M. Ali, M. Mokhtar, M. Chiesa, P. Armstrong, A cooling change-point model of community-aggregate electrical load, Energy and Buildings 43(1) (2011) 28–37. [50] G. Parnis, Building thermal modelling using electric circuit simulation, Master’s thesis, University of New South Wales, Australia (2012). [51] T. Oke, Boundary layer climates, Routledge (2nd edition), 1978. [52] J. R. Howell, R. Siegel, M. Meng, Thermal radiation heat transfer, CRC press New York, NY, 2011. [53] B. Bueno, G. Pigeon, L. Norford, K. Zibouche, C. Marchadier, Development and evaluation of a building energy model integrated in the teb scheme, Geosci. Model Dev. 5 (2012) 433–448. [54] UPC, Pearls design system new buildings energy and water benchmarking study, Tech. rep., Urban Planning Council (2010). [55] J. L. LaSalle, Abu dhabi real estate market overview, Tech. rep., Jones Lang Lasalle (2013). [56] R. International, Comprehensive cooling plan (ccp), Tech. rep., RTI International (2011). [57] ASHRAE Handbook of fundamentals, ASHRAE, 2013. [58] Abu dhabi municipalitys energy efficiency program - demand-side management in existing buildings, Unpublished internal report, Abu Dhabi Municipality (November 2011). [59] W. Jrges, The heat transfer at a flat wall (der wrmebergang an einer ebenen wand), Beihefte zum Gesundh.-Ing. 1. [60] A. Afshari, A new model of the urban cooling demand and heat island application to vertical greenery systems (vgs), Energy and Buildings [submitted]. [61] B. Bueno, J. Hidalgo, G. Pigeon, L. Norford, M. V., Calculation of air temperatures above the urban canopy layer from measurements at a rural operational weather station, Journal of Applied Meteorology and Climatology 52 (2013) 472–483. [62] I. Panagiotou, M.-A. Neophytou, D. Hamlyn, R. Britter, City breathability as quantified by the exchange velocity and its spatial variation in real inhomogeneous urban geometries: an example from central london area, Science of The Total Environment 442 (2013) 466–477. [63] J. Businger, J. Wyngaard, Y. Izumi, E. Bradley, Flux-profile relationships in the atmospheric surface layer, Journal of American Meteorological Society 28 (1971) 181–189. [64] J. Louis, A parametric model of vertical eddy fluxes in the atmosphere, Boundary-Layer Meteorology 17 (1979) 187–202. [65] A. Monin, A. Obukhov, Basic laws of turbulent mixing in the surface layer of the atmosphere, Trud. Geophiz. Inst., Akad. Nauk. S.S.S.R. 151 (1954) 163–187. [66] E. Erell, D. W. T. Pearlmutter, Urban Microclimate: Designing the Spaces Between Buildings, Routledge, 2010. [67] A. Afshari, F. Schuch, P. Marpu, Estimation of the anthropogenic heat release due to traffic using btex measurements a case study in abu dhabi, Urban Climate [submitted].

31

Page 32 of 32