Spatio-temporal analysis of heating and cooling degree-days over Iran

In this research, temporal trends and change points in heating degree-days (HDD), cooling degree-days (CDD), and their simultaneous combination (HDD + CDD) were analysed for a 60-year period (1960–2019) in Iran. The results show that less than 20% of the study stations had significant trends (either upward or downward) in HDD time series, while more than 80% of the stations had significant increasing trends in CDD and HDD + CDD time series. Abrupt changes in HDD time series mostly occurred in the early 1980s, but those in CDD time series were mostly observed in the 1990s. The cooling energy demand in Iran has dramatically increased as CDD values have raised up from 690 ºC-days to 1010 ºC-days in the last 60 years. HDD, however, almost remained constant in the same period. The results suggest that if global warming continues with the current pace, cooling energy demand in the residential sector will considerably increase in the future, calling for a change in residential energy consumption policies. Monotonic trends in HDD, CDD and HDD + CDD time series were analysed in Iran. The series were also analysed for possible change points using the non-parametric method. Energy demand for cooling houses increased, while that for heating houses remained constant. Up to three upward abrupt change points were observed in the annual CDD series. Monotonic trends in HDD, CDD and HDD + CDD time series were analysed in Iran. The series were also analysed for possible change points using the non-parametric method. Energy demand for cooling houses increased, while that for heating houses remained constant. Up to three upward abrupt change points were observed in the annual CDD series.


Introduction
Increasing emissions of greenhouse gases after the Industrial Revolution caused warming of Earth's climate system (IPCC 2014). British Petroleum (2019) reported that world population may grow up to 50% in the next forty years and thereby energy demand may increase by 80%. About 40% of the world's energy consumption is allocated to the residential sector (United Nations 2016), leading to 36% of the world's CO 2 emissions (Delgarm et al. 2016;Naderi et al. 2020). By continuation of the current CO 2 emissions rate, global temperature increase by 3.2-4 ºC in comparison to pre-industrial level until 2050 is not unexpected (IPCC 2014). Some studies in the United States showed that 1 ºC increase in air temperature leads to 3-15% decrease in heating demand and 5-20% increase in cooling demand (EPA 2016; Obringer et al. 2020). To stop global warming's increasing trend and control climate changes, we need to signi cantly and sustainably reduce greenhouse gas emissions. One effective method for controlling greenhouse gas emissions is to enhance buildings' energy performance (Ferrara et al. 2014). An undeniable fact regarding greenhouse gas emissions is that this important task can only be done by cooperation of all countries.
Weather directly impacts buildings' heating and cooling energy demand (Atalla et al. 2018). Energy consumed for buildings' heating and cooling is in uenced by climate change. Nowadays, heating and cooling energy supply has become a critical problem in many countries. Also global warming and climate change are considered as a threat for human population by many scientists and they have been analysed via different viewpoints. Degree-day method is considered a very important index in evaluating climate changes (You et al. 2014) and one of the simplest yet useful methods for calculating the energy consumption of a building. This simple method is also utilized in other scienti c elds such as agriculture, architecture, entomology, etc. The degree-day method has been used for calculating energy consumption speci cally for building's heating and cooling since 1932 (Dufton 1934) and it is generally considered a trusted method for this purpose (Atalla et al. 2018). Adjusting living environment's temperature status with a reasonable consumption of energy along with effective environmental management can reduce expenses and maximize bene ts from resources for the sake of society (Christenson et al. 2006).
HDD and CDD can be studied using different temperature thresholds. The temperature threshold is de ned as the comfortable temperature for human so that there is no heating or cooling demand. The thresholds may differ depending on economic status, energy supply and human physiological needs. For example, heating and cooling temperature thresholds are respectively 18 ºC to 22 Roshan et al. (2017a) reported that only 18% of studied stations in Iran had thermal comfort bioclimatic conditions. Wachs and Singh (2020) projected that due to climate changes, residential cooling energy demand in northern cities of the state will increasingly grow by 2080.
The trend analysis and change point detection of climate parameters provide precious information for a physical understanding of climate change and its impact on our living environment. One of the most important challenges caused by global warming and climatic change is an energy consumption shift. Iran is a major energy supplier: the second oil and gas supplier and exporter in the world (before the USA sanctions against Iran in 2019). Yet, having more than 80 million population along with its old infrastructure made Iran also one of the major energy consumers. Energy consumption in the country is 4.4 times above the global average, making it one of the top ten greenhouse gas emitting countries (Gorjian et al. 2019). Achieving a sustainable development without considering energy sector is impossible. One of the most important potential factors of energy consumption increases in any region is temperature changes. Measuring cooling and heating degree-days can present a clear and exact image of building, city and region's thermal demands. It can also play an important role in reforming energy consumption patterns. In this study, monotonic analysis was performed on HDD, CDD and HDD + CDD time series for the rst time. Previous studies (Kadioğlu et al. 2001;OrtizBeviá et al. 2012;Shi et al. 2018;Spinoni et al. 2015) have investigated the long-term changes in HDD and CDD time series, but not abrupt changes. For an accurate uctuation analysis of these time series, pre-and post-breakpoint series trends were obtained and analysed. Return periods were also calculated using the best statistical distribution for the studied time series. Spatial analysis of the studied indices was carried out with Cokriging interpolation method. It should also be noted that in this study, the Pettitt test was used for the rst time in detecting multiple breakpoints in the time series. Our ndings can reveal climate change impacts on buildings' cooling and heating demand, thus contributing to decision-making and planning process for energy management.

Study area and data
The study area is the whole area of Iran with the total area of 1.648 million km 2 , located between 25º 00t o 39º 47´ N latitudes and 44º 02´ to 63º 20´ E longitudes. Iran has a complex topography with divers climates including hot deserts in the centre, high altitude mountainous regions in the north and west, at fertile plains in the southwest (and many other locations), and lowland humid areas in the north. Figure 1a demonstrates the location of the studied stations as well as Iran's digital elevation model (DEM) map. Figure 1b represents the spatial distribution of the country's average air temperature. The annual temperature difference in some areas reaches 25 ℃.
Iran's long-term records for 60 years (1960-2019) were collected from Iran Meteorological Organization. For this purpose, all the synoptic stations that had data since 1960 were selected to be used. The stations with more than 10% of missing data were excluded. Daily minimum and maximum air temperature longterm data for 60 years from 31 stations distributed over the country were used.

Data quality control
To assess the quality of the used data, some traditional tests were used. Observations indicated that in  (Grubbs and Beck 1972) was used to detect outliers in the time series. This test was hence selected for this purpose. In this procedure, when one or more data points at one station were detected as outliers, the data of the same year from one or more neighbouring stations were investigated. In the case, the data points of neighbouring stations did not match the outlier, the outlier was removed from the time series; otherwise, it was considered as a valid data point.

HDD + CDD index
Time series of HDD were prepared using the mean air temperature threshold equal to 18 ºC. Indeed, the HDD in a year is the cumulative absolute values of differences between mean daily air temperature from 18 ºC in the days when the mean air temperature is below 18 ºC. Here, the mean air temperature (T mean ) was calculated by averaging the minimum air temperature (T min ) and maximum air temperature (T max ) recorded for the same day.
In addition to the individual analysis of cooling and heating demands, their joint spatio-temporal distribution over Iran is also investigated through heating and cooling degree-day sums (HDD+CDD). The HDD+CDD index corresponds to the total amount of heating and cooling demand and overall outdoor thermal comfort. We use this index as an easy-to-understand indication of differences in heating and cooling demand in residential sector ((Petri and Caldeira 2015; Sivak 2008; Sivak 2009)). Following many studies carried out around the world to use the degree-day method to quantify the buildings' energy

Monotonic trend analysis
In order to carry out the monotonic trend analysis for annual time series of the three indices of HDD, CDD and HDD + CDD, the following approaches were used.
1. To investigate HDD, CDD and HDD + CDD time series trend for the last 60 years, the non-parametric Mann-Kendall trend test was used (Kendall 1975;Mann 1945). For the series with signi cant autocorrelation coe cients, the modi ed M-K method by Hamed and Rao (1998) was utilized. This is due to the fact that the existence of positive (negative) serial correlation will increase (decrease) the possibility of rejecting the null hypothesis of no trend while it is actually true (Tabari and Hosseinzadeh Talaee 2011). M-K test was calculated at 0.05, 0.01 and 0.001 signi cance levels. Trend line's slope was also obtained via the non-parametric Sen's slope estimator (Sen 1968), which is resistant to the effect of extreme values in time series and thus gives a robust estimate of trend magnitude (Tabari et al. 2015).
2. In order to detect the breakpoint in time series, Pettitt (1979) test was employed. For cases with signi cant abrupt changes at a 0.05 level, change percentage was determined for post-breakpoint against pre-breakpoint. Also for stations that had experienced abrupt changes at a 0.01 level, the trend line slope for pre-and post-breakpoint series was obtained using the simple linear regression.
3. Since the Pettitt test is only capable of detecting a single breakpoint during a time series, therefore, some modi cations were used here to develop this approach for multiple change points as follows: For those series that had one signi cant breakpoint (at a 0.01 level), the time series was divided into two subseries in the detected point. Then, the Pettitt test was applied again for both sub-series. This is repeated if there was a signi cant change point in any sub-series.

Return period
A frequency analysis was performed to estimate the magnitudes of the variable under consideration with respect to the following return periods: 2, 10, 25, 50, and 100 years. Sixty-ve statistical distributions were tested for the goodness of t using the Kolmogorov-Smirnov (K-S) test and the best distribution was used to estimate the return levels. The non-parametric K-S test does not depend on cumulative distribution function being tested. This test is also independent of the sample size. K-S test is a widely used method for evaluating the performance of marginal distribution functions in tting different meteorological variables

Spatial analysis
First of all, correlation among some available meteorological parameters (e.g. T min , T max , T mean , frost days and rainfall) was calculated using Pearson method to nd two variables with a strong correlation with the studied parameter in order to improve the accuracy of the interpolation results. For creating spatial distribution maps of the studied indices, Inverse Distance Weighting (IDW), Kriging and CoKriging geostatistical techniques were employed for the interpolation of the three studied indices over Iran. In CoKriging, we can bene t three auxiliary variables to improve the accuracy. In this study, the country's DEM as well as two other variables that had a strong correlation with the studied parameter were used as auxiliary variables. By applying cross-validation, performance evaluation criteria including Normalized Root Mean Square Error (NRMSE), Mean Absolute Percentage Error (MAPE) and Mean Bias Error (MBE) were obtained. The mathematical details of these equations can be found in related research references (Bayrakçı et al. 2018;Côté et al. 2018). Finally, the best model which can reasonably represent the region's topographical pattern was selected for preparing the nal map for spatial analysis of the studied indices.

Spatial analysis
Results of Table A. 1 shows the performance measurement of IDW, Kriging and CoKriging obtained by the cross-validation technique. In order to conduct the spatial analysis of the studied indices, three geostatistical methods of IDW, Kriging and CoKriging methods were utilized. The results revealed that Kriging and IDW methods were not su ciently accurate for the purpose of this study, as detailed below. In CoKriging, three auxiliary variables can be utilized. In this study, Iran's DEM as well as two other variables that were highly correlated with the studied index were utilized, giving highly accurate results.
Since only 31 synoptic stations were quali ed for the long-term analysis and this number of stations may not be su cient for a country with 1.648 million km 2 area and with a complex topography. That is why in regions with complex topography (i.e. a region including vast mountain ranges, valleys, plains and deserts.) IDW and Kriging methods do not yield results consistent with the region's topographic and climatic conditions. It should also be pointed out that air temperature and indices based on that (e.g. HDD and CDD studied in this research) depend largely on the topography of the studied region. They have indeed an inverse relation with elevation, being lower in highly elevated mountainous regions than plains, deserts and its neighbouring coastal atlands. The ndings of our study, however, demonstrate that IDW and Kriging methods are not capable of presenting accurate results and distinguishing between high and low elevation areas. IDW and Kriging have been used in many studies because of its simple usability and without considering its error and related compatibility with actual conditions of that area which may lead to unreliable results. Although using auxiliary variables in CoKriging method increase its accuracy compared to Kriging and IDW methods, the model still may not be able to distinguish the difference of degree-days between elevated areas and lowlands. Using the region's DEM as one of the auxiliary variables signi cantly enhances model's logical compatibility by identifying the topographical pattern of degree-days especially when the number and distribution of the stations are not su cient. The region's DEM was therefore used to improve the accuracy of the spatial analysis of the temperature-based indices. Our results show that in Kriging method, there were no signi cant differences between ordinary, simple, universal and disjunctive models, while in CoKriging method, some differences were observed.
The results indicated that for the spatial distribution of the studied indices, the disjunctive CoKriging model offered relatively better results compared to the other models. That's to say, the error for CoKriging method was 3 times smaller compared to IDW and Kriging methods. This model was then used for the spatial distribution analysis of the studied indices.
3.1.1. Spatial distributions of HDD, CDD and HDD + CDD Figure 2 shows the spatial distribution of annual HDD, CDD and HDD + CDD in Iran and Table 1 gives the average monthly HDD and CDD. The results show that there is no need for heating from June to September. As expected, the greatest demand for heating is during winter months (from December to February) and the greatest demand for cooling is in the summer. In general, the demand for heating in Iran is more than cooling. Figure 3 represents the box-plots of HDD, CDD and HDD + CDD. As can be seen in Fig. 3a, stations that are located in the southern tropical regions of the country had a higher energy demand for cooling than the other stations. In these locations, the highest CDD change range was also observed. On the other hand, the largest HDD change range were observed in the cold and mountainous regions, which are mainly located in the west of the country. According to Fig. 3b, the greatest range of HDD + CDD changes are seen in areas with cold or warm climates, such as the northwest and southwest of the country. In contrast, the change range is less in areas with relatively mild climates, such as coastal regions and some central areas of Iran. Table 2 represents the magnitudes of trends estimated by Sen's slope for monthly and annually HDD and CDDin selected stations. Table 3 gives the number and percentage of the stations with positive and negative trends as well as the number and percentage of stations with upward and downward abrupt changes. According to these tables, 19% of the stations showed signi cant annual trends in HDD time series at the 0.05 level. Three stations located in the west and northwest of the country experienced a signi cant positive trend in HDD time series, while three other stations located in tropical regions had a signi cant negative trend. About 80% of the stations had slightly positive or negative annual trends that were not signi cant, implying that the station's trend is approximately zero. It seems that positive HDD slopes mostly exist in mountainous regions that are speci cally located in the northwestern part of the country, while the negative HDD trend was commonly observed in tropical regions. The analysis of monthly trends showed that in the cold months of the year, especially November to January, the HDD trend was positive, while in the warm months of the year (April to June, September and October) the trend was negative in most stations. Cvitan and Sokol Jurković (2016) found a positive HDD trend for December in Croatia.

Trend analysis
According to the results of Table 3, 87% of the sites showed signi cant (α < 0.05) trends for CDD time series at the 0.05 level, which decreases to 74% for the higher signi cance level of 0.001. Among these stations, except Shahrekord station which had a signi cant negative trend, all the other stations experienced signi cant positive trends. The steepest CDD trend line belonged to the southwest of the country which is about 19.7 ºC-days per year. Table 2 also shows the output of the trend analysis for HDD + CDD time series in the selected stations in Iran. It can be suggested that all stations (except Zanjan) experienced a positive trend and 81% of them were signi cant at the 0.05 level. However, inspection of results revealed that such an increase is a result of increasing trends in CDD series. Increasing trends in HDD + CDD time series were intense in some of the stations. For example, trends in HDD + CDD time series in 42% of the selected stations were signi cant at a 0.001 level. The steepest trend line slope of 18.1 ºC-days was also calculated for the southwest of the country (Abadan and Ahvaz stations). Figure 4 shows the results of Pettitt test conducted for the three time series at the selected stations in Iran. The percentage change between the mean values of post-breakpoint and pre-breakpoint series along with the direction of abrupt changes are shown in the Fig. 4. Figure 5 shows the frequency of abrupt change years in 5-year periods. According to the results, all upward abrupt changes in HDD occurred before 1985, with 80% in 1981. All downward abrupt changes in HDD happened in the 90s. 32% of the stations had no abrupt changes in their annual HDD time series.

Breakpoints
Concerning CDD, the ndings denoted that (except Shahrekord) all the other studied stations experienced upward abrupt changes in CDD values. About 70% of the experienced abrupt changes happened in the 90s. The percentage of the CDD increase in post-breakpoint compared to pre-breakpoint is higher for the northern half of the country than for the southern half (Fig. 4b). From gures presented by Anjomshoaa and Salmanzadeh (2017) realized that from 1996 to 1998, Kerman city experienced a downward abrupt change during the heating season and an upward abrupt change during the cooling season, corroborating our results. They have also reported that the heating season begins later and nishes sooner in the recent years due to the impacts of global warming. Similarly, the cooling season begins later and nishes sooner than before.
According to Fig. 4c all the stations except two stations in the northwest of the country experienced an upward abrupt change in HDD + CDD time series, which in the 80s. Speci cally, 81% of the stations experienced their upward abrupt changes in 1981. It is noteworthy that at 25 stations, the abrupt change was so high that made their Pettitt test signi cant at the 0.001 level.

Time series trends before and after the abrupt change
Pre-and post-breakpoint series trends are also investigated via a linear regression. For the sake of brevity, we only considered abrupt changes that were signi cant at the 0.01 level and the gures are presented in Appendix. Four patterns are distinguishable in the trends of HDD, CDD and, HDD + CDD time series. The rst pattern is a positive trend in both pre-and post-breakpoint sub-series. The second one is a negative trend in both pre-and post-breakpoint sub-series. The third pattern is a positive trend in pre-breakpoint sub-series and a negative trend in post-breakpoint sub-series. In fact, the the breakpoint truned the direction of the trends from positive to negative. The fourth pattern is the opposite form of the third one. These four patterns for the three studied indices at some selected stations are presented in Fig. 6. A close investigation of the results revealed that among 16 stations that had a breakpoint in their HDD series at a 0.01 level, 13 stations had a negative trend after the breakpoint. Actually by air warming, the heating demand was decreased after breakpoints at these stations, but the reason behind why these stations experienced an upward abrupt change before 1985 is still unidenti ed for us. At one station (Shahrekord), pre-and post-breakpoint series trend is close to zero. Closely investigating the ndings of this research, it was found that at more than 90% of the selected stations, HDD trend was decreasing after 1980. In fact, Our results also showed that pre-and post-breakpoint CDD series trends were positive for the majority of the stations. Tavakol et al. (2020) also claimed that in Mississippi, USA, heat waves and hot days experienced an upward abrupt change in 1997; the trend before abrupt change was negative while the abrupt change caused a turning point in the trend and made it positive. It should also be mentioned that both Shiraz and Shahrekord stations experienced positive trends before breakpoints, but their trends after the breakpoints were approximately zero.
As explained before, HDD + CDD trends for all the studied stations (except Zanjan) over a 60-year period (1960-2019) are positive; however, the trends before and after the abrupt change were negative at 14 stations. In fact, upward abrupt change caused a positive trend for the period at these stations, while the trend was negative before and after the abrupt change. Anyhow, for some stations in different parts of the country, the trend line slope after the abrupt change was approximately zero (< 1.5 ºC-days). Only charts with an intense trend line slope after and before abrupt changes are shown. It seems that the HDD + CDD trend for most of the stations in the southern half of the country is positive before and especially after abrupt changes. In contrast, the majority of the stations with a negative HDD + CDD trend after abrupt changes are located in the northern half. At some stations, a turning point in trends occurred in HDD + CDD series.

Multiple abrupt changes in time series
The stations experiencing more than one abrupt change are investigated using a modi ed Pettitt test. Table 4 shows the year of abrupt changes and direction of changes for the stations experiencing more than one breakpoint. HDD time series had more than one breakpoints at 9 stations. The rst abrupt change for all 9 stations occurred in 1981 and had an upward direction. The second abrupt change was downward and generally occurred in 1997. Based on the results of Table 4, CDD time series in four stations experienced three breakpoints. At almost all stations experiencing more than one breakpoint in CDD, the direction of abrupt changes was upward. In other words, CDD increased step by step at these stations. The exceptional stations located in the west and northwest of the country experienced a downward rst abrupt change. Concerning HDD + CDD, both abrupt changes were upward for most of the stations, while stations in different parts of the country, the rst abrupt change was downward and the second one was upward. Conversely, what happened in Ramsar in the north and Mashhad, and Torbat Heydariyeh in the northeast of the country was an upward rst abrupt change and downward second one.
Although a few studies from all over the world have been found to discuss and analyse HDD and CDD trends, no researches are found to investigate abrupt changes in HDD and CDD time series for comparing with the ndings of this paper.

Return period
After a consecutive analysis procedure applied for the calculation of return periods, the spatial distribution of the magnitudes with respect to the four return periods are plotted in Fig. 8. The frequency and return period analyses of HDD indicated that the southern coastal strip of the country had the highest degree-days increase (hereafter de ned as a deviation from the average value of the baseline period) in return period. On average, HDD values of 10-, 25-, 50-and 100-year return periods for all over the country increased by 20%, 28%, 33% and 38%, respectively. Analyses of frequency and return periods for CDD showed that stations located in the northwestern part of the country experienced the largest increase with respect to their return periods among the other stations. The smallest increase in CDD return periods was calculated for southern tropical regions where about 20% increase is for CDD of a 100-year return period. The joint HDD + CDD index value sharply increases from 10-to 100-year return periods, while it is less the case for the two single indices.

Discussion
Although urban areas cover a small portion of Earth's surface, they generally have the highest energy consumption. Energy consumption is in turn one of the main reasons of climate change, however a large portion of energy consumed for heating and cooling is due to climate change (Wachs and Singh 2020). Besides, growing urbanization has made heat islands and led to signi cantly increased energy demands and greenhouse gas emissions. Population growth in Iran has been very rapid. In the rst o cial census in 1956, the population of Iran was 19 million, while in the last census in 2016, it reached 80 million. Also in 1995, only 32% of the population lived in cities, while in 2016, it was estimated that about 74% of the population lived in cities. An increase in the cooling demand in Iran by an average magnitude of 70°Cdays per decade was found, which can be attributed to anthropogenic in uences (e.g., climate change, urbanization). In contrast, HDD didn't change signi cantly. HDD + CDD shows a considerable increase due to the CDD increase. If global warming continues with the current trend, the demand for the energy supply of the residential sector will increase over the next decades. In developing countries such as Iran that use less renewable energies and rely on their fossil fuels, this level of increase in the demand for energy supply will cause various problems in the energy supply and impose great expenses to the country. Energy consumption in Iran is much higher than the global standards (Gorjian et al. 2019). Lack of renewal and enhancement of old, ine cient heating and cooling systems in buildings, a rapid population growth and continuation of global warming trend can cause an energy supply crisis in Iran. The country may encounter critical problems in energy supply, despite being one of the major energy (oil and gas) exporters. Tian et al. (2014) reported that between 1971 and 2010, gas consumption has tripled in buildings sector while oil consumption has slightly reduced. In the same period using coals for the heating in buildings was also halved. However, heating comes normally from using gas or other fossil fuels while cooling generally requires electrical power. It is also necessary to note that since Iran has enormous gas and oil sources, using fossil fuels for generating electricity is more economical, as 93% of the resources for electricity power plants in the country are supplied by fossil fuels (gas, gasoline and mazut), and the rest is supplied by hydropower. Therefore, any increase in the energy consumption, whether for heating or cooling, practically causes more emission of greenhouse gasses. Changing demand from fossil fuels to electricity can have considerable outcomes for energy policies and the environment (Li et al. 2012

Conclusions
This study investigated historical changes in heating degree-days, cooling degree-days, and their simultaneous combination. The results revealed that the highest heating demand and the lowest cooling demand belong to cold and mostly mountainous regions in the northwestern part of the country where the annual average HDD goes beyond 2380 ºC-days while annual average CDD is 276 ºC-days. On the other hand, for southern tropical regions, the coasts, the Lut desert and the central region of Iran (regions around Yazd), observations showed the lowest heating demand and the highest cooling demand. In Ahvaz and Abadan that located in the southwest of the country, the cities with the highest cooling demand, this demand has increased up to 3000 ºC-days due to global warming over recent years.
The temporal changes and their spatial distribution in HDD, CDD and HDD + CDD indices in Iran were also investigated in this study. HDD changes at the studied stations represent different patterns; some stations experienced upward abrupt changes, while others had downward abrupt changes. Although for the majority of the stations with abrupt changes, post-breakpoint series trend was negative, pre-breakpoint series trend was positive for some stations and negative for some others. 32% of the stations also experienced no signi cant abrupt changes in HDD series, while all the studied stations except Shahrekord experienced upward abrupt changes in CDD values. It is noteworthy that 94% of the stations experienced upward abrupt changes in HDD + CDD time series. Investigating CDD values, the trends for the majority of the stations were positive after abrupt changes. Regarding pre-breakpoint trends, with except of 6 out of 29 stations that had a negative pre-breakpoint trend, the trend line slope for the other stations were positive or approximately zero. The statistical signi cance level of 90% of the stations that experienced abrupt changes in CDD time series was the 0.001. All abrupt changes that happened in HDD + CDD series were upward and 93% of them occurred in the 80s. Most of the stations with a negative HDD + CDD trend after abrupt changes were located in the northern half of the country. By modifying Pettitt test, we found that 29% and 65% of the stations experienced more than one abrupt change in HDD and CDD series, respectively. Also, in almost 35% of the studied stations, more than one breakpoint was observed in HDD + CDD time series.
In Iran, energy is mainly supplied from fossil fuels and the country's average energy consumption is fourfold the global average. Our ndings showed that the energy demand for buildings' cooling signi cantly increased due to global warming. If global warming trend continues, electricity demand will therefore increase in the future. Taking the advantage of solar energy which has high potentials in Iran for generating electric power instead of fossil fuels can lead to reduced emissions of greenhouse gasses as well as controlling negative effects of climate changes. In addition, renewal and enhancement of buildings' energy sector in Iran that are generally old and low-e cient can play an important role in achieving an optimum energy consumption and reducing air pollution. The results of this study are bene cial for investors, designers and managers in the residential sector, as it provides a basis for analysis and prediction of buildings' energy performance. They can also be used for a better understanding of the impact of global warming on buildings' energy consumption.

Declarations
Declaration of Competing Interest   Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.  Time series trends before and after the abrupt change for a) HDD, b) CDD and, c) HDD+CDD at some stations. The rst case is for the stations that had positive trends in both pre-and post-breakpoint subseries. The second one is for the stations that had negative trends in both pre-and post-breakpoint subseries. The third case is for the series had a positive trend before the breakpoint and a negative trend afterwards. The fourth case is the opposite form of the third case. Mu denotes the average of sub-series.