Technical, economic and environmental optimization of district heating expansion in an urban agglomeration

In order to integrate large shares of variable renewable energy sources, district heating can play an important role. Furthermore, in order to increase the efficiency of district heating systems, interconnecting adjacent system could be socio-economically justified. In order to assess the economic and environmental consequences of the latter, a mixed linear integer optimization model was developed with the endogenous decision on the potential interconnectors. The case study was carried out for the city of Zagreb, Croatia. The results showed that all three studied interconnections are economically viable, while the socio-economic cost was 29.2% lower in the case of the implemented interconnectors, all other capacities being equal. Moreover, the optimal thermal energy storage capacity was found to be equal to 25 and 24 days of average heating demand in two alternative scenarios. Finally, compared to the reference case, the CO2 emissions could be lowered by 15.3%.CO2 savings derive mainly from better utilization of low carbon capacities after interconnecting the systems, as well as from installation of heat pumps and electric boilers. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Newly completed package Clean energy for all Europeans, presents an important step towards the implementation of the Energy Union strategy and it delivers to the Paris Agreement commitments for reducing greenhouse gas (GHG) emissions. In eight legislative acts, package addresses all 5 dimensions of the Energy Union: (i) energy security; (ii) international energy market; (iii) energy efficiency; (iv) decarbonisation of the economy; and (v) research, innovation and competitiveness. Increase in energy performance in buildings, 32% of renewable energy sources in EU's energy mix by 2030, building targets of at least 32.5% energy efficiency by 2030, robust governance system for the energy union and electricity market design are building elements of the Clean energy for all Europeans package [1].
Heating and cooling sector accounts for 50% of the EU's annual energy consumption. The largest primary energy source for heating and cooling is natural gas with 46% share, followed by coal, biomass and fuel oil with shares of 15, 11 and 10%, respectively. Renewable energy sources account for 18% of primary energy consumption, however, its potential is not even close to being fully utilized. Considering EU climate goals, heating and cooling demand is expected to fall by 42% by 2050 with a proportional reduction in CO 2 emissions. To achieve the latter mentioned, the EU identified actions in making the building renovations easier, increasing the share of renewables, higher use of waste heat from industry and involving consumers and industries. Mentioned actions were implemented in the first initiative in the EU that considered the energy used for heating and cooling in buildings and industry called Heating and Cooling Strategy 1 presented in February of 2016 by the European Commission. The strategy identifies actions that will transform heating and cooling sectors to smarter, more efficient and sustainable sector, while at the same time decreasing the energy imports, costs and CO 2 emissions. Mentioned actions are key actions incorporated in the Energy Union Framework Strategy 2 [2].
The basic idea of district heating (DH) is to use local waste heat sources in order to satisfy the local heating demand [3]. Traditional heat sources are combined heat and power, waste to heat and industrial plants. However, in the last decade, the higher share of renewable heat from biomass, geothermal wells and solar collectors was introduced in the DH system [4]. Smart energy system is the future energy system that challenges problems of integrating district heating with electricity sector and transport sector [5]. Smart energy systems can also integrate district cooling in the system [6]. Such a future system will be composed of smart district heating, electricity and gas networks which would be mutually coordinated to achieve an optimal solution for each sector and for the overall energy system [7]. List of studies mentioned in the review [5] concludes that district heating can, and will have an important role in the future energy systems, but they must undergo transition to low-temperature DH networks that will be part of the smart energy system. Hence, the needed transformation of the DH system can be defined under the term of 4th Generation District Heating Technologies and Systems (4GDH) [5]. As reviewed in Ref. [8], unlike the previous three generations, the 4GDH includes energy supply and conservation balancing, so the challenge of supplying more energy efficient buildings with space heating and domestic how water could be met. Moreover, the investments in reducing the district heating network losses should also be included, as well as the innovative and strategic integration of DH into the smart energy system [7]. Abilities needed for a future DH systems to be part of the smart energy system include: (i) supplying existing and newly built buildings with low-temperature DH; (ii) low grid losses; (iii) integration of renewable heat sources and utilization of recycle heat from low-temperature sources; (iv) participation in solving the problems of integrated fluctuating renewable energy sources; and (v) ensuring suitable planning and incentives for DH system operation, as well as strategic investment planning [8].
A study in Ref. [9] included design and operation analysis of eight efficient district heating markets in the EU. The study explored eight case studies from national policy frameworks to specific local business models to determine key factors that make them efficient, and to investigate if the best practices and business models could be transferred to the other Member States. The study has shown that district heating and cooling systems represent a powerful and cost-efficient measure to ensure the development of low-carbon and resilient local energy system. However, the study also concludes that there is no universal model to develop efficient DHC system, due to its distinctive architecture and reliance on the local dynamics. Hast et al. [10]. also adds to latter point concluding that it is not enough to perform country-level decision-making, since every municipality has its own local specific features. Therefore, case studies focusing on the municipality level are needed for a future development. Thus, case studies for several European countries have been tested.
This paper further investigates the impact of interconnecting geographically distributed district heating grids and expands methods used for the analysis of DH expansion in Sonderborg municipality [11]. Based on [3], proposed Sonderborg municipality interconnection represents tree-like structure in the development of district heating structures. Initially, district heating systems start in the form of a small tree-like structure, and by means of interconnections or expansion, it grows to ring-like structure, before it reaches the final fully meshed structure. Two main advantages of growth from tree to ring-like, or meshed structure are an increase in heat transfer capacity, and reduction of risk that large supply areas will not be served in case of major pipeline rupture [3]. Also, as addition to the previous research, impact of thermal energy Nomenclature C VO&M l;n;t Variable operating and maintenance (O&M) cost of plant n in hour t at location l EUR/MWh heat C fuel l;n;t The cost of fuel for plant n in hour t at location l EUR/ MWh fuel C CO2 t The cost of CO 2 emission in hour t EUR/ton CO2 K l;n The CO 2 intensity of the plant n at location l ton CO2 / MWh heat R ele l;n;t Revenues from electricity sales of plant n in hour t at location l EUR/MWh ele L n Electricity-to-heat generation ratio of plant n MWh ele /MWh heat q l;n;t εR þ Heat power of plant n in hour t at location l MW heat C cap n;l εR þ The lower bound of the connecting pipe capacity at location l MW q cap conPipe;l εR þ The connecting pipe capacity at location l MW E l;t District heating energy demand in hour t at location l MWh heat s ch l;t εR þ Heat charge of heat accumulator in hour t at location l MWh heat s dis l;t εR þ Heat discharge of heat accumulator in hour t at location l MWh heat h st;n;l Heat storage efficiency x n;l;t εf0; 1g Binary variable: 1, if unit n in hour t is on, otherwise 0 c n;l;t εf0; 1g Start-up cost activation binary variable; 1 if the unit n in hour t was turned on, otherwise 0 a l εf0; 1g Binary variable (do invest, don't invest) S l Annualized sunk investment cost in the case of piping investment decision at location l EUR H One hour (used for relating the energy generated (MWh) by using a certain power (MW) of the plant during the period of 1 h) T A set of hours in a year L A set of district heating locations N A set of different energy plant generators storage placement is investigated, with a goal of lowering costs in a form of smaller pipe size for interconnection of multiple district heating networks. Thermal energy storage acts as buffer between heat demand and supply, and allows for maximization of DH system flexibility and performance, as well as utilizing the potential of connecting more buildings to existing system, which leads to decrease in operating costs and emissions [12]. Furthermore, the impact of interconnecting multiple DH networks and the use of thermal storage as a buffer to an operation of Combined Heat and Power (CHP) units was investigated. Proper thermal energy storage placement and design increases operating hours of a CHP unit, as well as allowing for a better management to shift production of electricity to hours when electricity price is higher [12]. Moreover, steam demand for industrial purposes has been added to a model and tested on a case study for a City of Zagreb and surrounding regions of Zagreb County. Further three paragraphs describe the current state-of-the-art in the research on district heating. Two recent papers from Telebi B. et al. [13], and Sarbu I. and Mirza M [14]. provide extensive research on modelling and optimization of DH systems. In Ref. [13] authors provide a new approach in DH categorization including geographical variables, size of the DH system, DH heat density and consumer demand. The study also includes a research review of system and component modelling approaches, as well as a review of the state of the art in optimization-based. For each component in the DH system, the study provides equations and assumptions on how to model the whole DH system. Both, the deterministic and heuristic optimization solving techniques are reviewed, with deterministic methods being divided into complex, simplified and predictive models. Complex models include software such as Energy Plus 3 or TRNSYS 4 used for a building level modelling for a generation of consumer demand profile, but they are dependent on a quality of an input data and come with a high computational cost. Simplified models often use the assumption that simplifies building characteristics, while predictive models are used for demand profile generation using predictive methods such as Kalman filter or artificial intelligence methods [13]. Similarly to Ref. [13], authors in Ref. [14] divide optimization techniques into deterministic and heuristic, where heuristic techniques include genetic and evolutionary algorithms. Main advantages of heuristic approach are nondependency on restrictive assumptions of the objective function and that results are not dependent on a initial values of the decision variables. However, this advantages come to a cost of a longer calculation times. Authors conclude that researchers generally take the approach of DH system optimizations under steady-state conditions using deterministic approaches, which generally don't guarantee a global optimal solution. However, recent research shows that researchers tend to use heuristic optimization techniques such as simulated annealing (SA), particle swarm (PSO) and ant colony optimization (ACO) [14].
Related to the modelling for the optimal configuration of the system, a number of researches have used deterministic or heuristic approaches. Sanaei S. M. and Nakata T. in Ref. [15] used generalised reduced gradient method (GRG) to calculate optimal choice of the system components solving minimum total cost equation. Boro D. et al. in Ref. [16] created mathematical mixedinteger linear programming (MILP) model that resulted in the optimal configuration and operation of the DH system. The objective function included total annual cost (TAC) and total CO 2 emissions. EnergyPLAN modelling tool and Matlab algorithms were developed to optimize the expansion of a district cooling system in the tropical region [17]. The authors in the latter case used a simulation tool for modelling the expansion of the district energy system. Haikarainen C. in Ref. [18] developed a model as a MILP problem for the optimization of a DH network, including the objective function made out of all construction and operation costs of new network components. Mertz T. et al. in Ref. [19] created a model described as mixed-integer nonlinear programming (MINLP) problem that optimized configuration and design of the DH system using objective function summed up of investment costs and operational costs, and used GAMS 5 software to solve optimization problem. Li. L et al. in Ref. [20] created a multi-objective MILP model to optimize distributed energy sources (DER) with neighbouring DH network. Authors used Matlab 6 software for solving the objective function composed of the sum of TAC and CO 2 emissions. Delangle A. et al. in Ref. [21] developed a MILP model in GAMS and solved for optimal expansion of the existing DH system with the objective of maximising cost savings or minimising greenhouse gas emissions using CPLEX solver.
Regarding the modelling of the DH network, researchers used both the heuristic and deterministic approaches. Earlier research by Weber C. et al. [22] divided multiobjective optimization problem into two problems, optimizing based on the minimum costs and CO 2 emissions and for the optimal choice of the DH system such as new heat pumps and characteristics of the pipe dimensions and pipe water temperatures. Authors used multiobjective evolutionary algorithm for solving the objective functions. Craus M. et al. in Ref. [23] developed a hybrid genetic algorithm to solve for the optimal network extension taking into the account the most profitable consumers and using the constraints on optimal pipe path. Zeng J. et all in Ref. [24] used hourly substation demand to develop an optimization model for DH network. Authors used a genetic algorithm to solve the objective function composed of the total annual costs and using results for a real DH network determining the optimal diameter of pipes. Dobersek D. and Goricanec D. developed a new approach for the optimal DH network design using a simplex algorithm to solve objective function composed of a sum of capitalised costs [25]. Morvaj B. et al. created a MILP model that included a selection of optimal technologies, network design and optimal operating conditions to satisfy the demand for electric and thermal energy. Authors implemented the model in AIMMS 7 and used CPLEX solver to find the optimal solution for the objective function composed of total costs and CO 2 emissions [26]. Bordin C. et al. proposed a new methodology by using graph theory and MILP model, using objective function to maximise earnings of connecting new consumers to existing DH network, simultaneously minimising the operation costs of the DH network [27]. Doroti c H. et al. in Ref. [28] analysed benefits of integration of DH and cooling systems. Authors developed hourly based multiobjective optimization model capable of defining supply capacities and their operation on a one-year period. They conclude that integrated DH and cooling systems can operate with the same CO 2 emissions, but with a lower total discount cost. Vesterlund M. and Toffolo A. developed a multiobjective evolutionary algorithm for the minimisation of DH network expansion costs. Authors used the real case of the DH network expansion in the town of Kiruna. Results of the case study show that the optimal solution is compromise between investment costs for the new pipes and the new technologies, and operating costs in the DH network [29].
The mentioned paper rarely address the interconnection of already existing district heating systems. Furthermore, no mentioned paper in this area considered steam production optimization, on top of heat and electricity generation considerations. Finally, models from the papers mentioned in the state-of-the-art analysis do not take into account the Not-in-my-backyard movement due to which people resist having nearby installed generators that have flue gases that can pollute the air. In order to fill the mentioned gaps, and to improve the current models that were used for calculating the feasibility of interconnections between different district heating systems, this paper brings three main novelties.
To sum up, the main novelties in this paper are: (i) developed mixed-integer linear optimization model with endogenous decision-making on investments in interconnecting different district heating systems, (ii) considering thermal energy storage at each of the potentially interconnected district heating systems in order to potentially reduce the needed diameter of the connecting piping (iii) heat generation supply from the newly connected areas fulfilled only by heat with no emissions produced at the location, forestalling any potential air pollution sources To investigate aforementioned impacts, the case study for a City of Zagreb and surrounding regions of Zagreb County, the city of Velika Gorica and Zapre si c was carried out. The aim is to show if the interconnection of the multiple district heating systems can result in fuel-saving and/or reduced CO 2 emissions, while simultaneously ensuring the economic competitiveness of the connecting energy system.
Subsequently to the introduction, a model description is presented in chapter 2, while the case study is described in chapter 3. In chapter 4 one can find results and discussion on case study for the City of Zagreb, followed by conclusion in chapter 5.

Methods
A linear mixed-integer programming approach was used to develop a model used in this paper. This model presents an improvement in the model developed in Ref. [11]. The main improvement deals with the connecting piping dimensioning, which is now automatically calculated by the model. Furthermore, investments in the connecting piping between different district heating systems were modelled as a mixed-integer problem now, with binary variables used for the options of having or not having the investment in the connection piping. Moreover, steam demand for industrial purposes is introduced into the model.
The predefined technologies in the model were gas combined heat and power plants, gas boilers, electric boilers, electric heat pumps, heat accumulators and ancillary steam boilers. The model allowed for capacity extension and/or building interconnectors between different district heating systems. Sector coupling between power and heat sectors was modelled via CHPs, heat pumps and electric boilers. Revenues from electricity sales of CHPs on a day-ahead power market was modelled as an income in the district heating system, represented as the term.
R ele l;n;t in the objective function (1). Consequently, the total CO 2 emissions and fuel costs from CHPs was assigned to the district heating system, when calculating the costs of CO 2 emission allowances. Including electricity sales as revenue in the system avoided the need for complicated calculations of splitting the cost of fuel and CO2 emissions on the power and heat systems. Electricity needed to power the electric heat pumps and electric boilers were priced according to the day-ahead power prices, with additional transmission and distribution fees. Due to the combined production of electricity, heat and steam in CHP units, the coefficient that calculates the loss of electricity produced per each MWh of steam sold for industrial purposes, was used to account for the amount of electricity lost as an additional variable cost. The objective function was set to minimize the total socioeconomic costs of the system. It included the annualized investment cost, fixed and variable operating and maintenance costs (O&M), unit start-up costs, carbon and fuel costs, as well as the annualized sunk costs associated with the investment decision for building the piping connection (1).
Where N represents a set of technologies installed at a location l, while T represents a time set (8760 h of a year). The case study developed for this paper had four different district heating locations.
Several different constraints were used in this model. The heat energy balance is represented by (2). q n;t;l H ¼ E l;t þ s ch l;t À s dis l;t ; ðc t 2 T; c n 2 N; c l 2 LÞ; Where H ¼ 1 hour. Hence, all the generators n located at a certain location l needed to generate enough heat to cover the total heat demand plus the additional demand for charging the heat accumulator minus the heat discharged from the heat accumulator. The heat balance equation for heat accumulators is represented by (3). s level n;l;t ¼ s level n;l;tÀ1 þs ch n;l;t Às dis n;l;t . h st;n;l ; ðc n2N; c l2L;c t 2TÞ (3) The capacity of each plant needed to be sufficient to generate the heat in every hour of the year (4) 0 q n;t;l q cap n;l ; ðcn 2 N; ct 2 T; c l 2 LÞ Units needed to produce within the minimum and maximum power (lower and upper bounds) [30]. The minimum power for CHP units, if turned on, was 20%, while for boilers the minimum production was 0. q minimum n;l ,x n;l;t q n;l;t q cap;maximum n;l x n;l;t ; ðcn2N; c l2L;ct 2TÞ Binary variable c l;n;t was activated if the plant n was just turned on in hour t (6). c l;n;t ! x l;n;t À x l;n;ðtÀ1Þ ; ðcn 2 N; c l 2 L; ct 2 TÞ A maximum discharge (and charge) of the storage during 1 h minimize X l2L X n2N X t2T h C VO&M l;n;t þ C fuel l;n;t þ C CO2 t K l;n À R ele l;n;t L n þ c l;n;t , C startup l;n q l;n;t H þ C cap l;n þ C FO&M l;n q cap l;n þ a l S l i ; ðcl 2 L; cn 2 N; ct 2 TÞ was constrained to 20% of the total storage capacity (5).
s dis l;t H q cap storage;l ; ðc l 2 LÞ Investments in the connecting piping between different district heating systems were modelled with binary variables. More precisely, there was a sunk cost S l associated with the investment as soon as the decision was made. Technically, the binary variables a l were modelled as shown in (6).
The model described here was implemented in Python 3 with Gurobi solver. The optimization was run on a personal computer with Intel Core i7 4-core processor, 16 GB RAM and 500 GB SSD disc.

Case study
The City of Zagreb is the capital and the largest city of Croatia. The total city area is 641 km 2 with a population of 790,017, based on the last census from 2011 [31]. Heating for Zagreb district heating network is being produced in CHP units TE TO Zagreb and EL TO Zagreb. District heating network is divided into three hot water networks and two steam networks. Steam network is mostly being used for industrial purposes with its small part being utilized for space heating purposes. The total length of the Zagreb district heating network is 274.41 km, while the total installed heating power is equal to 1,349.23 MW h . Zagreb district heating network is supplying around 100,000 customers [32]. Production units are gas-fired units with the possibility of using extra-light fuel oil as a replacement fuel [33].
TE TO Zagreb power plant has seven electricity/heat-producing units, and those can be observed in Table 3. Units K and L are combined cycle gas turbine power plants (CCGT) commissioned in 2001 and 2009, respectively. Additionally, Unit C is also utilized for electricity and heat production. Units VK3, VK4, VK5 and VK6 are heat generation boilers for a hot water network. TE TO power unit is also equipped with thermal energy storage of 750 MWh capacity [34]. Additional to electricity and heat production units, five units produce steam for industrial purposes. TE TO power plant produces steam with constant parameters of 230 C and 9 bar, while the produced amount depends on the consumers' demand [35]. Mentioned units K, L and C, besides electricity and hot water production, produce steam with a capacity of 140, 70 and 70 t/h, respectively. Two additional steam boilers, units PK3 and M produce steam with a capacity of 64 and 60 t/h, respectively.
EL TO Zagreb power station has a nominal power output of 90 MW e and 491 MW h and numbers four units for joint production of electricity, how water and steam, two units for hot water generation and two additional units for steam generation. Unit specifications can be seen in Table 3. Units A, B, H and J are used for joint production of electricity, hot water and steam, while Units G and K are utilized for heat generation only. Units M and N are utilized for steam production for industrial purposes. Steam produced has constant parameters of 225 C and 16 bar, while the produced amount depends on the consumers demand. Since the first of January 2018, units A and B are not being utilized since they are not in compliance with the regulations on emissions [36]. Due to the non-profitable revitalization of the Units A and B, planning and construction of additional CCGT CHP unit with a power output of 150 MW e and 114 MW h started in 2018. Commissioning is expected to happen in 2021 [37]. The main goals of the revitalization of the existing system are guided by the efforts to increase energy efficiency, especially in the part of production and distribution. It is of high importance to lower generation and distribution losses, which are especially high considering its value of 16.32% [33]. In Ref. [38] authors made detailed comparison between district heating systems in Zagreb and Aalborg, using comparative analyses. They conclude that a number of good practices used in Aalborg case, such as the importance of lowering specific heat and water losses and charging customers according to the cubic metres of hot water used, could improve Zagreb DH system operation.
The City of Velika Gorica is located 16 km south of the Zagreb. The city area is 328.65 km 2 with a population of 31,553, while the municipality counts 63,517 inhabitants [39]. Heat consumption of a city is 197.34 GWh and the district heating system supplies 32% of the demand [40]. Close urban area is covered with 14 district heating networks, only partially interconnected. The total length of the DH network is 9.84 km, while its total heating power output is 69.6 MW, and heating units are mostly gas-fired. In 2016 DH system produced 59.18 GWh of heating energy. Idea of interconnecting Zagreb and Velika Gorica DH networks was studied more than ten years ago, and conclusions was that interconnection was not profitable. However, market prices of the equipment and changes in technologies, as well as results of this study, might be used to revaluate the proposed interconnection [33].
Zapre si c is a town in northwest Croatia and is part of the Zagreb County. It covers an area of 52.60 km 2 , while the inner-city area is 18.96 km 2 . Population equals 25,223 inhabitants and it is the thirdlargest city in terms of population in Zagreb County [41]. District heating network consists of 8 individual heating systems which are not connected by a common district heating network. Boilers are mostly gas-fired with a possibility to use extra-light heating oil as fuel. Heating systems counts 2,372 customers which represent a share of 9.4%. Total "district" heating system is 2.37 km long and has 20.36 MW of installed heating power. Yearly heat production for the year 2016 was 16.53 GWh [33].
Following the validated methods used in analysis on Sonderborg municipality [11], a significantly updated model used in this paper focuses on the analysis of the district heating networks of City of Zagreb and two cities in Zagreb County, City of Velika Gorica and City of Zapre si c (see Table 1). Three scenarios were investigated. First one represents a baseline scenario with current installed capacities, while two additional scenarios represent future district heating system that will boost renewable energy utilization for district heating systems. Hence, the second scenario represents future district heating networks with no interconnections between Zagreb (Zagreb South and Zagreb North), Velika Gorica and Zapre si c, while the third scenario considers the possibility of mutual interconnections. Future scenarios represent the year 2025. Scenarios analysed are listed in Table 2.
Starting capacities in future scenarios are constrained by some production units not following the regulations on emissions. Data on development plant of EL TO Zagreb from document [37] suggests that only units commissioned after 2000 will be producing. According to decommissioned units for EL TO power plant case, which was commissioned before 2000, and due to missing data on TE TO power station development plan, we assumed that TE TO units commissioned before 2000 are also decommissioned by 2025. Details on production units can be seen in Table 3. Dual values for heat and electric capacities suggests the difference in starting values for the scenario I, and scenarios II and III, where first values represent starting capacities (sunk costs or forced investments) for the scenario I, while second value represents starting capacities for scenarios II and III. The difference between the latter values arises due to different plants decommissioning.
Furthermore, in all the scenarios the assumed CO 2 price was 25 EUR/t, while both electricity transmission and distribution fees were set to 20 EUR/MWh, respectively. The discount rate for all the technologies was set to 6%, except for connecting piping which was set to 4%, as it was considered to be a long-term infrastructure investment that does not share the same risk as to the investments in different generation plants. Finally, the economic figures on different technologies can be seen in Table 4.
The potential investment in connecting piping consisted of the fixed part and the linear part dependent on the installed capacity, as seen in equation (1) and Table 5.

Results and discussion
Results on heat generation for three different scenarios can be observed in Table 6. In scenarios II and III, natural gas boilers were replaced by electric boilers and/or heat pumps in both Zapre si c and Velika Gorica. The latter was the result of the implemented policy Table 1 General data on district heating systems in Zagreb, Velika Gorica and Zapre si c for 2016 [33].  Table 2 Description of investigated scenarios.
Scenario Description I Baseline, business as a usual case scenario II Future district heating networks with the possibility of using new thermal energy storage, electric boilers and heat pumps, but not allowing for interconnections III Future district heating networks with the possibility of using new thermal energy storage, electric boilers and heat pumps and allowing for interconnections

Table 4
Economic parameters used in the study [42]. that no emissions should be allowed in the newly connected areas, in order to remove any potential source of the air pollution, a common complaint nowadays of the citizens. Units that were generating in scenario II are shown in Table 7. One can observe that only the unit K, M and N from EL TO, and units K, L and M from TE TO were modelled as sunk costs (forced investments). To compensate for the missing capacity for EL TO, based on [37] and assumptions for loss of capacities for TE TO, the addition of CCGT CHP unit (Block KKF) in EL TO was the optimal result of the model . Based on the model results for the EL TO power station we can observe: (i) the addition of CCGT CHP unit with heat capacity of 114 MW and electric capacity of 150 MW; (ii) the increase in heat capacity of the Block K from 121 MW to 249.4 MW; (iii) the addition of 66.4 MW heat capacity electric boiler (North ZG EB). It is worthwhile to mention that HEP group, a heat and power generation company, is planning the construction of the new CCGT CHP block in EL TO power station, and our model showed that the investment is socio-economically justified.
Results on the TE TO power station capacities showed: (i) the addition of the 51.9 MW heat capacity electric boiler (South ZG EB); (ii) the increase in thermal energy storage capacity from 750 MWh to 4381 MWh and (iii) the addition of 294.5 MW heat capacity heat pump (South ZG HP).
As previously mentioned, to boost electrification and lower fossil fuel usage, the model considered the use of electric boilers and/or heat pumps with the possibility of additional thermal energy storage for Zapre si c and Velika Gorica locations. Hence, the results on generation capacities for Velika Gorica showed the replacement of natural gas-fired heating boilers with an electric boiler and a heat pump, with installed capacities of 9.7 and 16.6 MW, respectively. The case for the city of Zapre si c showed the replacement of 9 MW heating boiler with an electric boiler and a heat pump with both having the heat capacities of 4.2 MW. Additionally, the thermal energy storage of 2.8 MWh capacity is added in Zapre si c.
Similar to scenario II, the model allowed additional installed capacities in scenario III, as well as the investment in additional interconnectors between district heating systems. Hence, the potential connections between Zagreb North and Zapre si c, Zagreb South and Zagreb North, and Zagreb South and Velika Gorica were assessed (Table 8).
Based on the results for the EL TO power station we can observe: (i) the addition of CCGT CHP unit with the heat capacity of 114 MW and electric capacity of 150 MW; (ii) the increase in heat capacity of Block K from 121 MW to 415.7 MW.
Results on TE TO power station capacities showed: (i) the increase in thermal energy storage capacity from 750 MWh to 4124 MWh and (ii) the addition of 51.36 MW heat capacity heat Table 5 Economic cost parameters for district heating interconnectors. Cost included a fixed part of the investment (a sunk cost) and incremental cost dependent on the optimal interconnector capacity.

Fixed annualized investment (EUR)
Piping annualized cost (EUR/MW)   pump (South ZG HP). Results on generation capacities for Velika Gorica showed the replacement of natural gas-fired heating boilers with an electric boiler and a heat pump, with installed capacities of 2.8 and 13.5 MW, respectively. The case for the city of Zapre si c showed the replacement of 9 MW heating boiler with an electric boiler and a heat pump with heat capacities of 2.4 and 3.7 MW, respectively.
Regarding the interconnections, one can notice that all proposed connection were economically viable, and their capacities were (i) 430.21 MW capacity connection between Zagreb North and Zagreb South; (ii) 10.54 MW capacity connection between Zagreb South and Velika Gorica and (iii) 2.98 MW capacity connection between Zagreb North and Zapre si c. Heat flows between Zagreb South and Zagreb North aggregated on a daily basis can be seen in Fig. 1, while heat flows between Zagreb South and Velika Gorica can be observed in Fig. 2. Additionally, a smaller capacity connection between Zagreb South and Zapre si c can be seen in Fig. 3. One can see that all three of the interconnections are used for heat transfer in both ways, depending on the different period of the year. Furthermore, it can be noticed that interconnections are used less during the summer season.
Regarding the economical results from the model, one can observe that total system cost decreases from 117 m EUR for the scenario I to 105.8 and 82.8 m EUR for scenarios II and III, respectively. These values can also be expressed as savings for future scenarios (II and III) regarding the business as usual scenario (I). In Table 9, the additional economic and environmental data is shown. One can note that future systems, both in the case of interconnecting district heating systems, and not interconnecting them, can be significantly cheaper than the current one. Results show that total system cost for scenario II is 11.2 m EUR lower, while the savings for the scenario III are 34.2 m EUR. Moreover, CO 2 emissions are lower in both future scenarios. In scenario II CO 2 emissions are 48.9% lower, while for the scenario III reduction in CO 2 emissions are 15.3%. The main reason why the CO 2 emissions are much lower in the second scenario, is the difference in installed capacities between the second and third scenario. In the second scenario there is total of 447.5 MW of installed capacities in form of heat pumps or electric boilers, while in the third scenario only 74 MW of newly installed heat pumps or electric boilers are available. That difference is also due to the model allowed connection in the third scenario, which results with the higher CCGT CHP production, which in turn gives higher CO 2 emissions. The CO 2 emissions for CHP units were calculated using the efficiency method [44].

Conclusion
In this paper, a mixed-integer linear program with endogenous decision-making on building district heating interconnectors was successfully developed. The model was applied for the case of the City of Zagreb and the two adjacent smaller cities. There are three main conclusions that arose from this study.
First, when allowed in the model, it was optimal to install both heat pumps and electric boilers in all three geographic areas when interconnectors were not allowed. Electric boilers were covering peak demand while heat pumps were covering a significant share of the heat demand. In scenario III, the one in which interconnectors between the regions were allowed, electric boilers generation was almost negligible. On the other hand, the interconnectors allowed for much larger utilization of cogeneration units, making the system more profitable.
Second, it was socio-economically optimal to invest in all three possible interconnectors, both within two district heating regions of Zagreb that are located next to each other, as well as interconnector between the cities of Zagreb and Zapre si c that are distanced almost 19 km. The latter points out that urban agglomerations should more thoroughly assess the possibility of interconnecting its systems as those could bring significant benefits.
Third, the thermal energy storage significantly increased its size in both scenarios II and III. The optimal capacities amounted to the 25 days of the average heat demand in Scenario II and 24 days of the average demand in Scenario III.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.