Design of district heating networks in built environments using GIS: A case study in Vitoria-Gasteiz, Spain

Journal of Cleaner Production 349 (2022) 131491 Available online 22 March 2022 0959-6526/© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/bync/4.0/). Design of district heating networks in built environments using GIS: A case study in Vitoria-Gasteiz, Spain Mikel Lumbreras a , * , Gonzalo Diarce a , Koldobika Martin-Escudero a , Alvaro Campos-Celador b , Pello Larrinaga a a ENEDI Research Group, Energy Engineering Department, Faculty of Engineering of Bilbao, University of the Basque Country (UPV/EHU), Pza. Ingeniero Torres Quevedo 1, Bilbao, 48013, Spain b ENEDI Research Group, Energy Engineering Department, Faculty of Engineering of Bilbao, University of the Basque Country (UPV/EHU), Avda. Otaola 26, Eibar, 20600, Spain ARTICLE INFO Handling Editor: Mingzhou Jin Keywords: GIS Industrial waste heat Data-driven model LiDAR District heating ABSTRACT The efficient integration of high levels of industrial waste heat in low temperature district-heating networks is a promising technique that requires specific methodologies for its satisfactory implementation. This paper presents a novel methodology for assessing the energy and economic feasibility of new district-heating networks in existing urban areas for the integration of industrial waste heat sources. The methodology consists in an innovative multistep procedure using geographic information systems and data analysis tools, combining georeferenced data about buildings, industries and roads. The spatial distribution of the analysis area is divided into smaller buffers and grids, as a result, the routing design of the pipelines that makes up the district-heating topology is obtained under several assumptions. The methodology provides the most suitable area choice for the deployment of a district-heating, also implemented with a multi-step algorithm for routing the pipelines of the network. This methodology is applied to a particular case study located in Vitoria-Gasteiz (northern Spain). Different configurations for the district heating network are obtained with lengths of the network varying from 8 to 27 km. Payback values near to six years are achieved in most of the district-heating network configurations. The maximum payback period obtained within the configurations is 8.5 years. An economic sensitivity analysis is presented for the proposed optimal district-heating network configuration. The proposed methodology could be replicated for different case studies as long as the input data is available to the user. 1. Introduction Energy consumption in buildings currently accounts for around 40% of the total energy consumption in the European Union (EU) (P´ erez-Lombard et al., 2008). In the particular case of residential buildings, 57% of the total final energy consumption is used for space heating and 25% for domestic hot water (DHW) (Balaras et al., 2005). More than 50% of this energy consumption is nowadays fulfilled with natural gas and electricity (European Commission, 2019). Therefore, the implementation of alternative energy sources in buildings is vital to maintain a sustainable environment in cities and achieve the objectives of carbon neutral environment for 2050 (European Commission, 2017; Zhang et al., 2020) in the EU. District Heating (DH) networks can provide an efficient alternative to individual installations in densely populated urban areas (Christian Holmstedt Hansen, 2018). Today, DH networks provide more than 13% of the heating energy to buildings in the EU (Lund, 2007). Besides, new DH networks enable the connection of low grade and decentralized renewable energy sources (Lund et al., 2018; Wen et al., 2021; Yılmaz Balaman and Selim, 2016), which offers the possibility to employ waste heat streams as a sustainable heat source. There are different sources of waste heat available in cities or near them. Industrial residual thermal energy —commonly referred to as industrial waste heat (IWH) — is one of the most common. Recent studies show that the amount of IWH currently rejected to the environment accounts for 20–50% of the industrial energy consumption across EU (Brueckner et al., 2014). From that, 18–30% of it could be re-used in a technically feasible way (Brueckner et al., 2014; Vance et al., 2019). Since the emission temperature of around 65% of this IWH remains below 200 ◦C (Ankur Kapil, Igor Bulatov, Robin Smith, 2017), its reutilization for electricity production is hindered; however, these * Corresponding author. E-mail address: [email protected] (M. Lumbreras). Contents lists available at ScienceDirect Journal of Cleaner Production journal homepage: www.elsevier.com/locate/jclepro https://doi.org/10.1016/j.jclepro.2022.131491 Received 18 August 2021; Received in revised form 8 February 2022; Accepted 20 March 2022 Journal of Cleaner Production 349 (2022) 131491 2 temperatures are feasible for DH systems (Fit´ o et al., 2020; Moser and Lassacher, 2020), as demonstrated in the literature. For instance, in a case study (Ziemele et al., 2018), the use of IWH was able to cover from 4% to 12% of the system load in the heating season of a DH network located in Latvia. Similarly, the potential for integrating IWH (and solar thermal systems) into the DH of Germany was studied (Pelda et al., 2020). A great potential for exploitation was identified. Nonetheless, due to the spatial dimension of the problem, the use of georeferenced data is essential for the assessment and preliminary design of DH solutions. The spatial configuration of the buildings connected to the DH and the road network will define the length and, consequently, the investment cost for the DH deployment. These specific design requirements can be faced by Geographical Information Systems (GIS), which has been proven to be efficient for different applications at urban scale (Ali et al., 2018; Marques-Perez et al., 2020; Torabi Moghadam et al., 2018), including DH network deployment. As representative examples, the potential of extension of DH networks in United States based on density demand analysis and high-resolution GIS data was analyzed (Gils et al., 2013). This study concluded that heat distribution costs and, consequently, economic feasibility of these system, are strongly dependent on the heat demand density applied. Another reference (Chicherin et al., 2018) proposed a GIS-based model to combine annual cost analysis of DH systems with georeferenced data in order to help in the decision-making process of DH network planning. The method was applied in the DH network in Omsk (Russia) and the paper highlighted the advantages of using GIS for urban planning. Another study (Nielsen and M¨ oller, 2013) developed a GIS-based analysis to examine the potential for expanding DH in Denmark by developing a detailed cost analysis. The paper concluded that some high-density areas in Denmark are suitable for installing DH networks. Finally, in the last reference analyzed (Untern¨ ahrer et al., 2017), a GIS based methodology was proposed to spatially assess the integration of DH networks in urban energy systems, based on georeferenced data of the geothermal energy resource and the road network, among others. All these studies showed the advantages of combining GIS tools with traditional cost analysis in DH networks. Nonetheless, no studies are found that present the complete methodology for an economic feasibility study of a new DH system. When the related literature is assessed, it can be observed that most of the articles are devoted to non-constructed areas. To date, DH systems have been usually designed for new urban developments (Nguyen et al., 2020; Tol, 2015), where the buildings have low heating energy consumption, and the dwellings are structured in an orderly manner. In these cases, the routing of the pipelines of a new DH network is usually very intuitive and manageable. However, the current rate of increase of new buildings’ construction is around 2% per year in the EU (Eurostat, 2020). This means that more than half of the building stock in a mid-term range is already constructed. Accordingly, energy supply scenarios for the following years need to consider energy efficiency actions in already occupied areas and, in turn, simple design procedures that are capable to take those aspects into account are also required. As a result, this study presents a GIS-based multistep methodology suitable to design new DH networks within a built environment. An industrial waste heat source is introduced as a fixed spatial constraint. The method considers the available road network for routing the DH pipes, a working scale that considers small sub-areas and a novel multistep methodology for graph creation and algorithm routing. The procedure is then applied in a particular case study in Vitoria-Gasteiz (Basque Country, northern Spain). It is intended to provide a simple method to preliminary assess the economic feasibility of new DH systems fed with IWH that contribute to more sustainable cities. 2. Methodology This paper presents a novel methodology for the optimization of the design of routing process of the pipelines in a new DH network. The proposed methodology starts with the determination of the two following boundary conditions: (i) Identification of the IWH available near to the new DH network and (ii) definition of the heating energy demand in the buildings to be covered by the DH network. There are various options for the calculation of each of the boundary conditions and the methodologies chosen for the case-study presented in this paper will be shown in Sections 3.1 and 3.2. If the boundary conditions are predefined, the methodology presented in the following lines would be valid for any other case-study. Thus, this section will consider that the IWH available is already characterized and that the demand in all the buildings in the region are also available. Then, Section 3 will explain the methodologies chosen for this case-study. The methodology starts with the definition of the areas in which the algorithms for the routing of the pipelines will be applied. This first step is detailed in Section 2.1. Then, Section 2.2 explains the multistep combination of the algorithms employed to define the optimal route of the network’s pipelines and finally, the economic metrics used to evaluate the different DH configurations are outlined in Section 2.3. Two main software were used in the study. An open-source GIS software, QGIS (QGIS Development-team, 2020), was employed for the activities related with georeferenced variables, while R Studio was applied for data processing and cost calculations (Core-Team, 2013). 2.1. Surface lay-out: buffers and grids Once the heat supply source (the industrial facility) was selected, a circular area of 1.5 km radius around the factory was set, to ensure that the heat demand within that area was larger than the waste heat from the factory, enabling all waste heat available to be exploited. This area was divided into smaller sub-areas in order to calculate the DH network. The reason for this division is that, although DH viability studies are usually made at regional or municipal scale, the design of the pipeline Nomenclature Acronyms EU European Union DHW Domestic Hot Water DH District-Heating RES Renewable Energy Sources WH Waste Heat IWH Industrial Waste Heat GIS Geographical Information System DSM Digital Surface Model DTM Digital Terrain Model PB Payback MST Minimum Spanning Tree Parameters Q n Building demand [MWh/year] dNPV Discounted net present value [EUR] C i Cash-flow r Discount-rate [%] DEM Nominal demand density [kWh/m 2 ] F Number of floors in a building [-] S Horizontal surface of a building[m 2 ] P Ratio between heated surface and the total surface LN Linear Heat Density [MWh/m] Ø Pipeline Diameter [m] M. Lumbreras et al. Journal of Cleaner Production 349 (2022) 131491 3 routing for the network is developed at a lower scale. However, reducing the working scale to individual buildings would require an excessive computational cost, so a middle ground solution was taken defining the mentioned sub-areas within the area of study. In order to produce them, two different types of partition were tested: buffers and grids (Fig. 1). Regarding buffer distribution, the operation consisted of estimating the heating demand of circular buffers, with the IWH source at the center of the buffer (green triangle in Fig. 1a). Buffers of various sizes were calculated, in order to optimize the result. The size (diameter) of these buffers differed in 100 m between them. The use of these buffers implied that the resulting DH network was evenly distributed around the IWH source, covering all the residential buildings inside the circular buffer. This approach intends to find correlations between the economic results and the geographical variables in the system but is not focused on finding the optimal configuration of the DH network. For grid distribution typology (Fig. 1b), the main area of 1.5 km was divided into a grid formed by squared cells. Each cell in the grid comprises a building or group of buildings. The heat demand of every cell was estimated by Eq. (2) (Section 3.1), in order to subsequently define the route of the DH network. Since the result depends on the employed cell size, different values were tested: from 50 ×50 m to 500 ×500 m, with steps of 50 m. For each grid size, a DH network configuration was proposed. Every cell was linked with the adjacent cells and the path of the network only connected those buildings showing the best economic results. A grid distribution of 150 ×150 m is shown in Fig. 1b. For the grid distribution case, different adjacency levels were considered, in order to find the optimal path for the DH network. Starting from a particular cell, the path that the DH network might follow (i.e., the next cell) depends not only on the characteristics on the immediate surrounding cells, but also on cells that are far from it. However, this implies that, in order to perform the calculations for each specific cell, all the remaining cells that comprise the total surface should be assessed, which is not feasible due to computational reasons. Therefore, calculations with different adjacency levels were performed. This idea is illustrated in Fig. 2 and explained next. The starting cell represents the industry location (red cell in Fig. 2). The first adjacency level includes the eight cells that directly surround it. The second adjacency level is formed by the cells that surround the first level. The subsequent levels are defined in the same way. Now, in order to determine the next cell that the path of the DH will follow, the first adjacency level assesses eight different possible options. The incorporation of a second adjacency level analyzes 32 different options and so on. Note that those cells that have been covered once by the network are considered to be not useable for the next steps. According to performed preliminary attempts, four levels of adjacency were included in the present study; the incorporation of higher levels entailed an excessive computational cost. As a result of this section, the buildings that may be connected to the new DH network are established. In the case of buffer distribution, this means all the buildings within these buffer areas, and, for the grid distribution, it means all the buildings included in the optimal path that is defined by cells’ algorithm. Different configurations will result in different buildings. The next section explains the multistep methodology followed for the routing of the pipelines of the network. 2.2. DH network definition using routing algorithms The goal of this section is to select the optimal DH network configuration. For this purpose, the length of the network is estimated using a multi-step methodology combining different routing algorithms. Both primary and secondary sides of the DH network are defined for each configuration proposed in the previous section. Thus, the lengths are computed based on graph theory methods (Berge, 2001). The objective of this routing step is to determine the minimum distance connecting all the buildings selected from previous section and using the existing road network. Under this framework, the buildings under study are transformed into vertices corresponding to the centroid of the surface of each building. The pipelines of the primary side of the DH are usually constrained by the road network, as pipe ditching is thus facilitated. The used algorithm comprises the following three steps that are sequentially applied: Fig. 1. Selected area of study divided in (a) buffers; (b) in a grid of 150 ×150m. Fig. 2. Adjacency levels concept: grid showing 4 levels of adjacency. M. Lumbreras et al. Journal of Cleaner Production 349 (2022) 131491 4 1. Graph Creation (Delaunay triangulation). This algorithm enables an optimal connection of all the buildings/points. 2. Routing (Johnson algorithm) of the network connecting the buildings. This algorithm enables the identification of road segments to join all the connections from the triangulation process. 3. Calculation of the Minimum Spanning Tree (MST), enabling the calculation of the minimum DH network lengths. 4. Calculation of Economically Optimal DH Configuration The combination of these four steps is illustrated in Fig. 3. Each step is detailed in the following sections. 2.2.1. Graph creation The buildings marked as optimal for DH deployment in Section 2.1 will serve as input for this point. These buildings are converted into points or vertices, using the centroid of the shape of each building. In this first step, the objective is to create a graph in which all the buildings are inter-correlated. For this purpose, the Delaunay triangulation (Delaunay B., 1934) is used to connect the buildings in each of the cases, forming a set of vertices and edges that make up a set of triangles. Delaunay triangulations maximize the minimum angle of all the angles of the triangles in the triangulation, avoiding triangles with one or two extremely acute angles that are not desirable during some interpolations and mathematical processes. The edges of these triangles join each pair of points that correspond to the centroids of the buildings, which must be analyzed in the following steps. The set of vertices and edges in the system forms a planar graph, which ensures that buildings that are far away from each other are excluded. Thus, each triangle connects three buildings, forming three pairs of buildings. The output from this triangulation process is multiple pairs of buildings connected by the Delaunay triangulation that will be used as input for the following process. Fig. 3. Routing Algorithm process scheme: Graph Creation (a), routing between buildings (b) and calculation of the minimum spanning tree (c). M. Lumbreras et al. Journal of Cleaner Production 349 (2022) 131491 5 2.2.2. Calculation of routing between buildings The objective of this section is to define the real distance by the road network for all the pairs of buildings connected by the Delaunay triangulation of Section 2.2.1. The weight of the edges from the graph created using the Delaunay triangulation corresponds to the Euclidean distance between the nodes or buildings. However, this distance does not correspond to the actual road network. Accordingly, Johnson’s algorithm (Johnson, 1977) is used to find the shortest path between all pairs of buildings in the affected area, taking into account the road network. In real circumstances and taking into account the fact that the distribution pipelines (primary side) of the DH network will couple with the current road network, the distance to get from one point to another could be much longer than the Euclidean distance. The real routing distance between two buildings will be the sum of the road segments between the two points and the distance between the building centroid to the closest road. Thus, a new weight is obtained for each pair of buildings in the planar graph, with the Euclidean distance between them. This will be the input for the following step. 2.2.3. Calculation of minimum spanning tree The objective of this final step is to define the configuration of the DH network by the calculation of the Minimum Spanning Tree (MST) from the set of vertices and edges resulting from Section 2.2.2. For this purpose, Kruskal’s algorithm (Kruskal, 1956) is applied, which enables the connection of all the buildings (vertices) with the shortest path weighted by the real distance between the different buildings connected by real roads. The calculation of the MST ensures that the DH configuration obtained from this multistep algorithm is the shortest network and, consequently, the optimal one in economic terms. The primary network connects the heat source with the substations of the building, whereas the secondary network in a DH connects the substations with the facilities of the building. The length of the primary side of the network (L PRIM ) is made up of the sum of the lengths of the different arcs in the MST. Furthermore, the sum of the distance of the translation of the centroids corresponds to the length of the secondary side (L SEC ) of the network. This connection can be made using an algorithm in (QGIS Development-team, 2020) introduced with the plugin from (Conrad et al., 2015). 2.2.4. Calculation of Economically Optimal DH configuration This final section aims to identify the optimal DH configuration among the different scenarios proposed in previous section. For this purpose, a simplified economic assessment is proposed in which only the largest expenses and cash flows are included. According to (ETI, 2017), the largest costs are reached for civil works, with 36% of the total costs (including the planning and development of the project). Moreover, the pipelines that need to be introduced in the trench after the civil works also take up an important part of the total costs; both adding up to over 50%. Since this analysis does not include operation costs, the results should be considered as preliminary. However, this approach is accurate enough for an initial feasibility study of any new DH network. For the comparison of all the cases, a simplified payback (PB) period is proposed. According to (Harris, 2018), the PB period of this kind of system is calculated as shown in Eq. (1): Payback =Initial Investment Energy savings Eq. (1) where the energy savings represent the costs saved by not using the heat supply system that is currently in use in those buildings. As there is no exact information of what type of heat supply technology exists in each building, it is assumed that all the buildings used a natural gas condensing boiler, with a nominal thermal efficiency of 90%. This approximation value for the efficiency is a standard value for boilers working at their nominal value. 2.3. Economic sensitivity analysis of selected DH network In this final step and with the objective of analyzing the economic sensitivity of the selected network, an extended economic assessment is carried out for the optimal network. This economic assessment is based on the general economic metric, discounted net present value (dNPV) which is calculated by the following equation (Eq. (2)). dNPV =∑ T t=1 Ci (1+r)iEq. (2) where F i is the net cashflow during period t, and r refers to the discountrate. The sensitivity analysis is based on the variation of some of the most fluctuating operational variables and studying the effects of these changes in the dNPV. Positive dNPV indicates that the discounted economic returns are greater that the investment required. So, for this economic sensitivity analysis, discount-rate and natural-gas price are parametrized for the first 20 years. The values used for the case study are presented in Section 3.3 Thus, the selected configuration for the deployment of the network is economically analyzed in order to avoid interpretation errors and to cover a wide range of scenarios. 3. Description of the case study The methodology proposed along the article was applied to estimate the basic design of a DH network in the city of Vitoria-Gasteiz (249,176 inhabitants in 2018), located in the administrative region of the Basque Country (northern Spain). According to K¨ oppen-Geiger classification, this location corresponds with a C fb climate, referring to oceanic climate. The case study was applied to an area urbanized over various decades, with a predominance of buildings constructed in the 1970s. The study area comprises a 1.5 km radius circular buffer around the industry under study (specified in Section 2.1). This area is chosen to ensure that the heat demand is larger than the waste heat from the factory. As detailed in Section 2.1, this area was divided into two types of smaller sub-areas, buffers and grids. 3.1. Industrial waste-heat sources The first step in the presented case-study was the characterization of the waste heat sources available in the factories and the identification of a suitable industrial emplacement. The characterization of the IWH available in the framework of this case study was performed by using the results from a previous work done by the authors in the Basque Country (Spain) (Larrinaga et al., 2021). This work presents the characterization of the most important factories in the region by the application of a bottom-up methodology using commonly available measures, such as heat sourced burned in the factory. A factory is appropriate for the exploitation of its waste heat by a DH network if some conditions are met. First, high waste heat energy has to be available and high temperatures provide a higher-grade energy. The factory may be located near an urban core with high density demand. The low distance between the factory and the buildings connected to the DH network reduce the heat losses in the distribution pipelines and optimizes the use of the IWH. Finally, no large physical barriers have to be found between the factory and the buildings. These barriers would drastically increase the initial investment required for the network. Examples of physical barriers could be mountains or rivers. The chosen industry was a foundry (Fundici´ on Olazabal y Huarte, 2020) that has been in operation for more than fifty years. It is located very near to the city center, which shows a large energy density demand, and currently employs 54 workers. The estimated technical waste heat available from the different manufacturing processes is around 158 GWh/year (Larrinaga et al., 2021). The temperature distribution of the heat available in the foundry is shown in Table 1. From all the factories M. Lumbreras et al. Journal of Cleaner Production 349 (2022) 131491 6 analyzed, the factory chosen for this case study presented meets all the conditions mentioned above. Around 60% of the waste heat available in the foundry is at high or very high temperature range (>500 ◦C); whereas the rest of the energy is available at low-medium temperature. Even though high temperature range IWH is typically used for electricity production due to its higher quality, all the heat available in the factory is adequate for injection into DH networks. Since this is a preliminary study, it was opted to include also the high temperature IWH into the available residual heat for the DH network. 3.2. Buildings energy demand estimation The heat demand (Q n ) in each building was calculated trough Eq. (3): Qn[kWh] = DEMn⋅Fn⋅Sn⋅PEq. (3) There, DEM n is the nominal heat demand in each building [kWh/ m 2 ]; F n , the number of floors [-]; S n , the total surface [m 2 ] and P [-] is the ratio between the heated surface and the total surface of the building. Detailed information to calculate the parameters in Eq. (3) is provided next. In order to estimate the number of floors (Fn), the height of each building had to be first calculated. Typically, GIS information of buildings is available in shape files provided by the local authorities. They include location, total surface and polygonal shape of buildings; however, they do not usually provide information on the useful area or net height of buildings, as required in Eq. (3). To obtain them, LiDAR (Laser Imaging Detection and Ranging) data were herein employed (Sharma et al., 2021). LiDAR is an optical remote sensing technique that uses laser light to obtain a dense sample of discrete values of the height of earth’s surface. Amongst the different LiDAR data typically available, the Digital Surface Model (DSM) and the Digital Terrain Model (DTM) were herein applied. DSM presents real elevations on the surface, including manmade features, such as buildings or infrastructure, as well as natural features, such as trees. DTM is a model derived from DSM, where all the above-mentioned surface features have been removed, leaving only the bare terrain. Starting from these two datasets, the height of the involved buildings can be estimated with a raster calculator by means of Eq. (4): Net Heights =DSM −DTM Eq. (4) Since LiDAR data are usually available as *.las or *.laz files, they need to be converted into VRT files (raster files), which enable the raster operations between the DSM and DMT raster files. To do so, it is required that the distance between points (or cells in raster) is the same in both DSM and DTM data. This renders a new raster layer that includes the net height of all the features above the terrain. In the present case, the obtained grid size in both raster layers was 5 m. Then, a plug-in from SAGA for QGIS (Conrad et al., 2015) was used to assign location attributes from the shape file of the buildings to the heights file of the features. An average height was obtained for each building. Fig. 4 shows the main difference between these two datasets for the involved location (the specific details of the buildings included in the present case-study are provided in Table 1). GIS information of the buildings of the case-study was obtained from (Gobierno Vasco, 2020), while LiDAR data were obtained from CNIG (CNIG, 2020). The number of floors of each building (F n in Eq. (3)) were calculated dividing the net height obtained from LiDAR data by the distance between floors. The distance between floors was herein fixed to 2.7 m, slightly higher than the minimum distance between slabs fixed for residential buildings (2.5 m) (CTE, 2019). The result from the division was then rounded to its lowest value. Finally, the horizontal projection of each building (S n in Eq. (1)) was obtained from the information gathered in the GIS files. Therein, buildings are available as polygonal elements; thus, S n was obtained by calculating the area of those polygons. This horizontal projection was multiplied by a reduction factor (P)of 0.83, in order to obtain the useful heating area in each building (D’Alonzo et al., 2020). This way, external or internal walls and non-heated spaces such as elevators are considered. The building layer from (Gobierno Vasco, 2020) was clipped by the 1.5 km buffer, identifying 4555 buildings in the region. The buildings in the location included various types of constructions, which are detailed in Table 2. Table 1 IWH distribution by temperature range in the selected foundry, retrieved from (Larrinaga et al., 2021). Temperature[◦C] 200–300 ◦C 300–400 ◦C 500–1000 ◦C >1000 ◦C Total IWH [GWh] 48.77 14.10 73.14 21.57 157.55 IWH [%] 30.95 8.94 46.42 13.69 100 Fig. 4. (a) DSM and (b) DTM for Vitoria-Gasteiz (Spain). Table 2 Distribution of the buildings in the area of study divided by the type of buildings. Type of Building Number of Buildings Percentage [%] Generic Construction 3733 81.95 Lightweight Construction 150 3.29 Greenhouse 38 0.83 Warehouse 484 10.63 Buildings in ruin 2 0.04 Sacred Building 71 1.56 Canopy 68 1.49 Singular Building 9 0.20 M. Lumbreras et al. Journal of Cleaner Production 349 (2022) 131491 7 The total heat demand of the buildings within the area was calculated using LiDAR data and based on building heights, following the methodology explained in Section 3.1. The proposed methodology is only valid for residential buildings. The rest of the types of building shown in Table 2 might have different heating demand profiles. Regarding the nominal heat demand (DEM n ) for the buildings under study, two alternatives were employed to estimate it: 1. DEM1: Nominal density values provided by the Spanish Institute for the Diversification and Saving of Energy (IDAE) (IDAE, 2011). The data are based on various simulations of average building types in different climates; thus, different demand densities for space heating and DHW are achieved. According to (IDAE, 2011), the nominal energy demand for space-heating is 163.6 kWh/m 2 and 13.5 kWh/m 2 for DHW. 2. DEM2: In order to include heat consumption of newer buildings, a lower nominal heating demand density was also employed, according to (Ter´ es-Zubiaga et al., 2015). The used value was constant: 60 kWh/m 2 (40 kWh/m 2 for space heating and the rest for DHW). Thus, the obtained total heat demand of the residential buildings in the area reaches 2528.82 GWh/year with DEM1 heat density and 1273.97 GWh/year with DEM2 density. These results represent at least 10 times the total IWH available in the considered factory. This excess demand reaffirms the idea that the buildings that might be connected to the DH network have to be adequately chosen in order to obtain optimal economic assessment. One of the variables that best describes the energy needs of a specific area is the energy density in kWh/m 2 . The heat density in each of the buffers is calculated as the sum of all the demand in the buffer divided by the surface covered by the buffer. Within buildings with similar energy consumption, the higher buildings will show higher energy density. In fact, Fig. 5 shows the buildings layer classified by the energy density (Fig. 5a) and a 3D view (Fig. 5b) of the buildings in the area using a plugin by (WebGL technology and three.js JavaScript library, 2020), also classified by the heat consumption density. Thus, of all the buildings in the location, only those of generic and lightweight construction (Table 2) are included in the analysis. Other building typologies would require an additional demand characterization due to their singular heat consumption profiles. Regarding these two types of constructions, Fig. 6 shows the distribution of the residential buildings according to the energy demand density. 3.3. Economic assessment & sensitivity analysis On the one hand, the initial investment comprises the cost of trenching, pipelines and the connection to the installations of the existing buildings (substations). The equations Eq. (5) and Eq. (6) for this economic approach are taken from (Persson and Werner, 2011): Investment Pipelines =130 +2858⋅ ØEq. (5) Ø=0.0486⋅LN +0.063 Eq. (6) where Ø is the diameter of the pipeline in [m] and LN represents the linear heat density in [MWh/m]. In these equations, the investment costs are given by pipe length unit in [EUR/m]. Moreover, the cost for the DH substations, including the equipment, installation and connections is fixed at 2500 € /substation (Gudmundsson et al., 2013). On the other hand, the operational variables for the case study are set for the sensitivity analysis of the selected network. First, the natural gas price is expected to increase in the following years (Gao et al., 2021) and consequently values in Table 3 are proposed as for household consumers’ price. Even though the current natural gas price is affected by several factors, a mean value of 0.05 € /kWh is used as current reference. The other main parametric variable is the discount rate. This parameter is changed from 0% up to 10% using the values shown in Table 3. Thus, the selection of the optimal DH network configuration is carried out calculating the simplified PB period and the sensitivity analysis is calculated only for the selected network configuration. 4. Results 4.1. Buffer approach This section will present the results from the application of the methodology shown in Section 2 to the case study presented in Section 0. The application of routing algorithms explained in Section 2.1 allows the real path for the DH network configurations in each of the variants to be defined. The study area of the case-study has been distributed and divided in two different ways, representing two different concepts of network configuration: circular buffers and grids. The distribution of the study area into circular buffers aims to develop multidirectional DH networks from the IWH source as the initial point. Thus, all the residential buildings within the circular buffer are covered by the heating network, regardless of the overall heat demand or heat demand density of the buildings. As an example of the construction of the DH network using the routing algorithm presented in Section 2.1, Fig. 7 presents two network configurations for two different buffer sizes: 300 and 1000 m. The economic results obtained for the buffer’s method are presented in Table 4. The payback period for different buffer sizes and the two considered heat demand densities are outlined. Since the radial Fig. 5. Area under study. (a) Building planar layer; (b) 3D view. M. Lumbreras et al. Journal of Cleaner Production 349 (2022) 131491 8 character of the network makes its length increase exponentially with the size/radius of the buffer, the costs for the construction of a new DH network are completely influenced by the size of the buffer. A simplified payback period below 10 years is achieved within a buffer radius equal to or less than 600 m for the two heating demand density values, DEM1 and DEM2, as it can be observed in Table 4. With the DEM1 heating density, all the IWH available is consumed in the 400m radius buffer, whereas, with the DEM2 density, the waste heat from the industry manages to cover all the buildings within the 700 m buffer. No more demand can be covered in larger regions and, consequently, the economic metrics exponentially worsen from those buffer size values. From the buffer sizes indicated above, the payback variable leaves its linearity and begins to increase exponentially, mainly due to the large pipeline system necessary to cover all the buildings. In both cases, the payback shows a minimum in the 200 m buffer. From the analysis of the results of the buffer distribution, it can be concluded that the buildings that might be covered by the network have to be carefully and adequately chosen. Although the analysis made in the study based on buffer areas is interesting from a theoretical point of view, the real network would be constructed by the results obtained from the grid analysis. This algorithm allows directional DH network configurations to be built, with the economic optimization included. 4.2. Grid algorithm & definition of DH network configuration The grid distribution of the location aims to develop a DH network configuration with a unique main directionality, finding regions and buildings with optimal economic return. Unlike the buffer distribution, this algorithm finds the optimal direction of the network in function of the external conditions (heat demand, road length, etc.), so that only the Fig. 6. Number of residential buildings classified by the energy demand density for DEM1 and DEM2. Table 3 Economic variable Sensitivity analysis. Natural Gas Price [ € /kWh] Discount Rate 0.03 0% 0.05 3% 0.10 5% 0.15 10% Fig. 7. Configuration of the network. (a) Buffer of 300 m radius and (b) Buffer of 1000 m radius. Table 4 Payback period analysis for different buffer-radius [m] lengths. Buffer [m] 100 200 300 400 500 600 700 800 DEM1 4.5 3.5 4.3 4.5 5.4 8.5 13.1 17.9 DEM2 5.2 3.9 5.6 6.0 6.8 8.5 10.1 14.1 Buffer [m] 900 1000 1100 1200 1300 1400 1500 DEM1 24.5 30.5 40.6 50.9 59.9 70.7 80.1 DEM2 19.9 25.1 34.3 43.4 51.6 61.6 70.2 M. Lumbreras et al. Journal of Cleaner Production 349 (2022) 131491 9 most suitable buildings are connected to the heating network. The algorithm is developed in order to maximize the economic viability of the network, following Eq. (3), Eq. (4) and Eq. (5). As an example of these networks’ concept, Fig. 8 shows the layout of two DH configurations for two cell sizes. Overall, forty DH network configurations are modelled, including 10 grid sizes and 4 adjacency levels for each of the grid sizes. The distribution of the buildings varies among the different distribution sizes, and, in consequence, the network configuration will be different in each case. Even though each of the DH configurations is developed to be economically optimal, the change of boundary conditions by changing cell size will affect the lay-out of the network. Thus, simulating 40 different cases with different boundary conditions allow the most appropriate areas for the deployment of the network to be identified in order to match them up with the most repeated areas within all the simulated cases. Fig. 9 shows the heat map of all the 40 cases overlapped in the same image, identifying the most repeated buildings/regions with red marks. Fig. 9 shows that most of the DH configurations are conducted to the south and southwest of the industrial site. As it can be observed, some of the buildings turn out to be connected to the DH network in all the 40 configurations, whereas northern regions are not considered in any of the 40 cases. The initial distribution of the area of study into different sized grids defines the initial conditions for the simulations. Variation of the cell size affects the boundary conditions with different building distributions along the cells of the grid. Thus, the resulting path for the network may slightly vary from case to case (See Fig. 9). The roads used will vary from case to case and, as a result, the lengths and characteristics of the pipelines will also vary. As concluded from buffer analysis, the pipelines following the road network and their lengths and diameters determine the economic feasibility of the system. Thus, Fig. 10 presents the correlation between the lengths of the primary and secondary sides and the linear mean energy density of the network in each configuration. The primary side of the network connects the heat source, the factory in this case, and the substations, whereas secondary side of the network connects the substations with the final users/dwellings. From the study above was concluded that mean linear energy density is a critical parameter for economic feasibility of the system and this parameter results to be exponentially increase with long networks. Of the 40 configurations, the largest network would need around 27 km of pipelines, including the primary and secondary sides, whereas the shortest network would only result in 8 km. Two main curves are observed in Fig. 10 and some points present smaller networks for the same energy density. This is caused by the better road configuration of the cases with smaller lengths. Thus, heat network with large linear density and short networks are the most interesting from an economic view. The initial investment from the installation of the distribution pipelines completely limits the viability of the proposed system. The other variable affecting the payback period is the economic savings from not requiring the energy covered by the IWH. Fig. 11 shows the evolution of the two cash flows of Eq. (4), clustered by the grid size of the initial distribution. A clear divergence between the two linear tendencies of the cash flows considered in the study is observed. As long as the yearly demand covered by the network increases, the initial investment required for the installation of the network also increases, with a greater slope than economic saving function. The initial investment in the resulting configurations ranges from around 20 M € to around 72 M € ; whereas the economic savings from the primary energy savings range from 2.5 M € to around 10 M € . Besides, Fig. 12 shows a heat-map of the payback period, showing the 40 configurations under analysis for the proposed DH, distributed by Fig. 8. Configuration of the network for two grid sizes: (a) 250m and (b) 400m. Level 2 of adjacency in both cases. Fig. 9. Heat Map of the most covered areas by the different DH network configurations. M. Lumbreras et al.