A New Modeling Approach for Low-Carbon District Energy System Planning

Designing district-scale energy systems with renewable energy sources is still a challenge, as it involves modeling of multiple loads and many options to combine energy system components. In the current study, two different energy system scenarios for a district in Montreal/Canada are compared to choose the most cost-effective and energy-efficient energy system scenario for the studied area. In the first scenario, a decentral energy system comprised of ground-source heat pumps provides heating and cooling for each building, while, in the second scenario, a district heating and cooling system with a central heat pump is designed. Firstly, heating and cooling demand are calculated in a completely automated process using an Automatic Urban Building Energy Modeling System approach (AUBEM). Then, the Integrated Simulation Environment Language (INSEL) is used to prepare a model for the energy system. The proposed model provides heat pump capacity and the number of required heat pumps (HP), the number of photovoltaic (PV) panels, and AC electricity generation potential using PV. After designing the energy systems, the piping system, heat losses, and temperature distribution of the centralized scenario are calculated using a MATLAB code. Finally, two scenarios are assessed economically using the Levelized Cost of Energy (LCOE) method. The results show that the central scenario’s total HP electricity consumption is 17% lower than that of the decentral systems and requires less heat pump capacity than the decentral scenario. The LCOE of both scenarios varies from 0.04 to 0.07 CAD/kWh, which is cheaper than the electricity cost in Quebec (0.08 CAD/kWh). A comparison between both scenarios shows that the centralized energy system is cost-beneficial for all buildings and, after applying the discounts, the LCOE of this scenario decreases to 0.04 CAD/kWh.


Introduction
With the dramatic increase in the world's population in the last two centuries, the energy demand has also increased [1]. The U.S. energy information administration reported that the contribution of heating and cooling buildings to the total energy consumption is 40%, and mostly depends on fossil fuels [1]. Utilizing renewable energy resources leads to a decrement in the consumption of fossil fuels and, thus, the related emissions. Facing climate change, the general attitude toward energy consumption is to lower it, and also to reduce the emission of greenhouse gases. For example, the European Council decided to reduce the EU's energy consumption by 32.5% in 2030 compared to the 1990 levels by increasing energy efficiency [2].
Integrating renewable energy resources with different energy systems is one of the leading solutions to changing energy production systems [3,4]. Following many cities' sustainability goals, the share of renewable energy, especially in district heating, is growing [21]. Besides environmental aspects, the economic feasibility of energy systems is a crucial factor [22]. Tanguay [23] analyzed the economic performance of residential GSHP in North America's market. The study includes comparisons of equipment and energy costs in four provinces: Quebec, Manitoba, British Colombia, and Ontario. Results show that the average cost of GSHP for residential buildings in Quebec is 164 CAD per square meter of surface area. Aditya et al. [24] conducted a cost comparison of using GSHP compared to other technologies. An essential parameter in this study is the climate. Montreal, London, Singapore, and seven Austrian cities were chosen as case studies to show the effect of climate. A two-bedroom residential building was chosen as a case study and modeled in different cities. The results show that, due to the high heating demand in Montreal, and its low electricity costs, the city ranks as the cheapest city for GSHP installation among nine other case studies [24]. Furthermore, a comparison between Levelized Cost of Energy (LCOE) and Energy System Analysis (ESA) methods was made for decentralized heating and district heating systems [25]. Both methods show similar ranking results for district heating networks, whereas the results of decentralized heating systems vary.
Centralized district heating and distributed renewable energy systems have been widely studied in recent years. Dolla Rosa et al. [26] studied the effect of human behavior on load patterns in a low-temperature district heating. They also carried out a socio-economic comparison between district heating and distributed GSHP for a low-heat-density area. They showed that the low-energy district heating's levelized cost of energy is competitive with the GSHP-based scenario for low-heat-density areas. Moreover, according to this study, the optimal design of the network is of the utmost importance, since an optimal design decreases the average pipe size required, reducing the costs. Rämä and Mohammadi [27] studied centralized and distributed solar collectors in an existing district heating system. The showed that both scenarios are feasible; however, integrating renewable energy with a district heating system has higher economic feasibility than the distributed one. Although many works have been done in the field of renewable district heating and cooling networks in urban areas, there is still a gap in integrating building energy modeling with energy system sizing and energy distribution network design. In this paper, a novel methodology integrating these three modules and automizing the whole process has been presented. The current study is based on the Dominion Bridge area in Lachine-Est in Montreal, Canada. Figure 1 shows the location of the case study and the newly proposed designs. Developers intend to design an eco-quartier with an energy-efficient and cost-effective energy system. Two scenarios have been considered for the energy system of the mentioned area. Both scenarios share the same sets of buildings with the same characteristics, schedule, usage, and material. In the first scenario, a decentral energy system comprised of GSHP provides heating and cooling demands for each building, while, in the second scenario, a district heating and cooling system has been designed. Automatic urban building energy modeling (AUBEM) was used to calculate the energy demand of buildings. Then, a model was prepared in the INSEL simulation environment and implemented to design the energy systems for both scenarios. This model provides heat pump capacity and the number of required HPs, the number of PV panels, and AC electricity generation potential using PV. After designing the piping system and calculation of the heat losses for the centralized scenario, an economic assessment was carried out to choose the most cost-effective energy system scenario for the studied area.

Automatic Urban Building Energy Modeling (AUBEM)
One way to store the detailed geometrical information of buildings so that it is readable by computer programming is to structure the data into the City Geography Markup Language (CtiyGML) format [28]. CityGML is a text-based format similar to the XML format, which uses the hierarchal structure to store objects and attributes related to buildings. Since CityGML is readable using Python programming language, a CityGML file of the new design is used in the present study. In the current study, the term automatic UBEM refers to our software workflow that carries out the whole process of 3D building modeling for the energy demand analysis and the entire process of energy system modeling in an automatic manner. A Python code was written to perform the following tasks:

•
Extracting the building coordinates and building characteristics (building use-type, year of construction, etc.) from CityGML format; • Merging the building surface coordinates with their related characteristics in the same list and organizing the data to be suitable for energy simulation software; • Assigning the building materials and constructions to each surface based on building use-types and surface types, whether they are walls, roof surfaces, or ground surfaces; • Reading and organizing the occupancy, electrical equipment, lighting, and ventilation schedule and assigning them to each building based on building usage; • Feeding the data into EnergyPlus [29] for energy demand calculation. Figure 2 shows the AUBEM system workflow. The CityGML is parsed through Python to extract each building's ID and their related information and organizing them with regard to a hierarchical structure which is callable in the next steps. The hierarchical structure starts using building IDs as a root, followed by the building use-type and then their surface information as a subcategory. This hierarchical structure paves the way for accessing detailed building geometry and attribute information in the next steps for each building. EnergyPlus is a physics-based model that captures the full dynamic of the building performance to provide a detailed analysis of buildings [29]. Due to the sufficient flexibility and accuracy of EnergyPlus, it is used for building energy demand calculations.
EnergyPlus first defines the building zones and then creates 3D building models. Building IDs are used to name the zones. It is necessary to assign each building's surface to its related building zones. Therefore, each surface has a zone name which is defined based on the building ID and should be defined automatically before adding the surface coordinates into EnergyPlus.
The building materials and constructions are extracted from the National Renewable Energy Laboratory (NREL) website [30] and stored in a JavaScript Object Notation (JSON) format. The JSON format is parsed, and the building constructions and materials are categorized based on the building use-type, stored in a Python dictionary, and named as NREL building construction archetypes. The building construction archetypes are assigned to each building's surface based on the building use-type. Building use-types are extracted from CityGML. The occupancy, electrical equipment, lighting, and ventilation schedules are extracted from the Department of Energy (DOE) building archetypes website [31] for a large office, secondary school, small office, and midrise apartment buildings. Since DOE has limited building archetypes, we had to consider the most similar building archetype for each building use-type in the Lachine district. The large office archetype is considered for the Civic Center building type, the secondary school for school, the large office for a commercial, the midrise apartment for a residential, and the small office for office. The extracted archetypes are stored in a text file and are read through Python to be assigned to each building. In the last step, each surface, with all its related information, along with the building occupancy, electrical equipment, lighting and ventilation schedules, are fed to the EnergyPlus to calculate the heating and cooling demands. EnergyPlus positions the surfaces in their specific location and connects them to form a 3D building urban model. The other parameters used in the EnergyPlus model are summarized in Table 1. Moreover, the building surfaces are categorized into three groups, walls, roof surfaces and ground surfaces, and their related archetypes are assigned to each group accordingly. The other required parameters for accurate building energy demand calculation, such as occupancy, electrical equipment, lighting, and ventilation schedules, change dynamically based on building use-types, and are all automatically assigned to each building.

Energy System Model
A model for designing an energy system is prepared using INSEL 8.2. which is a simulation environment using a flexible graphical programming language comprised of ready-to-use blocks. The former INSEL distributor doppelintegral GmbH is liquidated and further development of INSEL is now pursued at Concordia University, Montreal, Canada. Moreover, new extensions or completely user-defined blocks can also be added and used. INSEL's main functionalities are meteorology, building modeling, and renewable energy systems' modeling [35,36]. The proposed energy system model obtains an hourly energy demand for both cooling and heating, the geometry of building (exclusively the roof dimension) and temperature as an input and provides heat pump capacity and the number of HPs required to best fit the demand curve, the number of PV panels and AC electricity generation potential using PV. The complete energy system model (ESM) includes different parts and features. However, in this study, only the heat pump section and the PV system will be discussed briefly. The ESM workflow starts by obtaining and analyzing the heating and cooling demand and ends with various results, including energy consumption, the system's coefficient of performance (COP), PV self-consumption and equipment selection. The summary of parts related to the scope of this paper is shown in Figure 3.

Heat Pump Model
The most common HP performance indicator is COP, a unitless parameter that can be interpreted as the thermal energy transferred by the HP to the HP electricity consumption. Since a heat-pump COP is highly dependent on heat-source and heat-sink temperatures, COP fluctuations during different conditions should be taken into consideration. In the case of GSHP, the heat transfer fluid temperature leaving boreholes and going to the building is used as a heat source input to the HP.
Based on the HP manufacturers' performance data, a correlation between HP energy generation, electricity consumption, sink, and source temperature can be found. Usually, reversible HPs' COP is higher than 1.0 due to the fact that they transfer heat from a heat source to a heat sink. That being said, by referring to the Carnot theorem (theoretical highest efficiency of an HP), it can be seen that by lowering the supply temperature or increasing the boreholes' return temperature, higher COPs can be achieved. In our case study, 40 and 12 • C are used in calculations for heating and cooling supply temperatures, respectively, which are consistent with the low-temperature heating/cooling concept.
In addition to improving COP, a lower heating supply temperature will lower heat losses in piping and storage and enables the heating system to achieve higher exergies. These values have been studied and investigated for low-temperature and ultra-lowtemperature district heating (DH) systems [7,[37][38][39]. Although 40 • C is sufficient for providing domestic hot water (directly or with supplementary heating) [40], the focus of the study is merely on heating and cooling demand and the electricity consumption of HPs' covering demand.
Previously, Weiler et al. [35] proposed a model for an Air Source Heat Pump (ASHP), which uses a polynomial fit to the manufacturer's data. In this study, the same approach is used and further developed for GSHP. In this method, two out of three HP parameters, including COP, electricity consumption and heat output, are defined as a function of supply temperature and heat source temperature (boreholes to HP), based on the HP performance provided by data manufacturers. Table 2 is extracted from Maritime Geothermal Ltd. W series GSHP catalogue [41], which shows the W900 70-ton HP's electrical consumption and COP for a constant supply temperature (heating) but various source temperatures. As a result, HP output heat can be determined, which will then be used to calculate the number of HPs required to cover the demand in each time step (1 h).

Hourly Analysis and PV System
As mentioned before, the ESM receives hourly demand data from urban building modeling. After preprocessing the load curves, some statistical analyses, including maximum demand, average demand (non-zero values), heating and cooling degree days, and different percentiles, will be calculated. Based on these data, the model will select a heat pump capacity among the available HPs in the energy library, which best matches the building load curve.
In each time-step (1 h), ESM takes heating and cooling demand (kW) as well as outdoor temperature ( • C) and solar irradiation, from an input file. The INSEL energy systems' library offers a comprehensive list of PV panels and inverters with various specifications and manufacturers, from which a 300 W panel with 17.24% efficiency and an appropriate inverter were selected. INSEL's available blocks were used to calculate direct and diffused solar radiation on a 31 degrees tilted panel, facing south, to feed the PV block. a maximumpower point tracker was added to ensure the system's optimum output. Finally, DC/AC electricity generation of PV system was used to cover HPs' electricity demand and the excessive portion will be exported to the grid. Figure 4 shows the setup for the PV system in INSEL.

Hydraulic and Thermal Modeling of the Distribution Network
As mentioned before, it is crucial to study the hydraulic and thermal behavior of the distribution network to understand the effects of the heat losses and pressure drops on the energy system. To do so, at first, a distribution network was designed for supply and return pipes. Figure 5 shows the configuration of the buildings and the designed network system for heating and cooling, supply, and return lines. A MATLAB code was then developed based on the graph theory to introduce the geometry to the code. Finally, the proper form of the related equation was derived and solved to obtain the network's mass flow rate and temperature levels.

Network Modeling Using Graph Theory
A graph represents a set of connected objects by nodes and edges. In the context of the district heating network (DHN), edges are interpreted as pipes, and nodes demonstrate the junctions and the consumer stations. To implement the concept of graphs in a mathematical model, an incident matrix A is defined, whose elements are 0, +1, and −1. The incident matrix shows the interconnection of the nodes and the edges. The rows and columns of the incident matrix represent the number of nodes and pipes, respectively. In the incidence matrix, +1 is assigned to the inlet node of a pipe, while −1 is assigned to the outlet node of the pipe. 0 is assigned to the nodes that do not have any interaction with a specific pipe [42]. The incident matrix is the mathematical representation of the network's configuration, which can be combined with the conservation equations to calculate the mass flow rates, and temperature distribution in the network (Equations (3) and (4)).

Hydraulic Modeling
For a pipe in the network, connecting two junctions, the one-dimensional, steady-state momentum equation was derived as Equation (1) [42] where P out is the outlet pressure, P in is the inlet pressure, g is the acceleration of gravity, . m is the mass flow rate in the specified pipe, f is the friction factor, L and D are the length and diameters of the specified pipe, respectively, and β is the total loss coefficient.
The associated equation for the calculation of the mass conservation could also be derived. Equation (2) shows the simplified mass conservation equation [42] ∑ where . m ext is the flow rate leaving each node. In consumer stations, . m ext is obtained by using the energy demand calculated in the previous sections ( . m ext = Demadn o f the building c p (T supply −T return ) ).
By employing the incident matrix, the above-mentioned equations will turn into the following forms [42] m D 4 ρ 2 π 2 ∑ k β k . A MATLAB code was written to solve the above equations using the Newton algorithm to solve the system of non-linear equations.

Thermal Modeling
Once the mass flow rate in each pipe is determined, the mentioned code uses the one-dimensional energy conservation law (Equation (5)) to calculate the temperature level in each node. The energy equation is discretized in time by the Euler backward method. An upwind scheme translates the temperatures at the boundary to the node upstream [42] where c p is specific heat capacity, S is the area of the pipe, L is the length of the pipe, P is the perimeter of the pipe, and U is the global heat transfer coefficient. In order to calculate the pipe diameters, Equation (6) can be used [42] where v is the design velocity that should not exceed specific limits in order to avoid noise and damage to pipes. Thus, the design velocity for the water system is considered to be 1.5 m/s for sizing calculations [42,43]. Once the mass flow rate in each pipe is determined, the pipe sizing can be completed.

Economic Assessment
Among various cost estimation methods, LCOE is a good indicator of the costeffectiveness of renewable energy systems. The LCOE of both scenarios is calculated to investigate the economic feasibility of central and decentral district heating and cooling systems in the current study. The LCOE is defined as the total lifetime cost of an investment divided by the cumulated energy generated by this investment. An alternative (but mathematically identical) approach is the definition by means of the net present value (NPV). The LCOE is the (average) internal price at which the energy is sold to achieve a zero NPV.
The formula to calculate LCOE is shown in Equation (7) where I t is investment expenditures in year t, M t is operations and maintenance costs in year t, F t is fuel expenditures in year t, E t is energy generation in year t, r is interest rate, and n is lifetime of the technology. As mentioned earlier, the location of the current study is in Montreal, Quebec. All the prices, including equipment and labor costs, are obtained from a local company named Marmott Energy (7106 St-Denis, Montréal QC, H2S 1S4, info@marmottenergies.com), based on their previous large-scale projects. Table 3 includes the detailed costs of different items involved in the economic assessments. Table 3. Costs of the items for economic assessment.

Demand Results
For the configuration of the buildings shown in Figure 5, the hourly heating demand is calculated. As mentioned above, different factors affect heating and cooling demands, such as building surface area, its application and location, etc. Figure 6 shows the heating and cooling curves of the buildings in the area. Based on this figure, building A has the highest heating and cooling demand, which is due to its heated surface area, which is the largest of all the buildings.  To better analyze the results, Table 4 shows the maximum heating and cooling for each building throughout the year.  Figure 6 shows that, during specific periods (around 3000 or 6500 h) there is a demand for both heating and cooling. Predefined heating and cooling and transitional seasons method can be implemented to avoid having heating and cooling demands at the same time, since this would not be feasible according to the high thermal inertia of the network. In this method, two specified periods are separately allocated to cooling and heating, and two other periods (one before autumn and one before spring) are defined as transitional seasons, which do not need heating or cooling. In this case, optic-variable walls (OVW) have high solar absorption during winter and low solar absorption during summer, and can be useful for providing thermal comfort during the transitional seasons [44]. However, in the current study, two separate distribution networks are considered for heating and cooling, so that the heating and cooling thermal demand is answered without causing problems. These calculated demands are the starting point for network analysis and energy system design. The following section discusses the energy system modeling results.

Energy System Results
ESM was used for each of the six buildings in the district plus the central district heating and cooling scenario. ESM results are summarized in Table 5. It is worth mentioning that supply temperatures for both heating and cooling are assumed to be constant, while fluid temperature leaving boreholes are dependent on the outdoor temperature. Moreover, ESM determines the number of HPs in each iteration that cover the maximum demand, and the required number of required HPs is added to the results. HPs are considered to have single-stage compressors and either work with 100 percent capacity or zero. Besides calculating metered parameters like electricity to/from the grid, only HP electricity demand has been considered. Since a single COP value cannot reflect the fluctuations in COP during different weather conditions, the Seasonal COP (SCOP) was used to understand HP performance better. SCOP for heating and cooling seasons are defined as follows where Q H , Q C are HP heat output (kWh) and and E H , E C electricity consumption (kWh) in the heating and cooling cycle.
Since the selected HPs are single-staged and the model requires them to cover the maximum demand, in almost all time-steps, a surplus of energy will be available. It is clear that sizing the energy system relying on maximum demand results in an over-sized system for most operating hours, which will lower the system's total efficiency and lifespan by causing numerous on/off cycles.
Comparing the central scenario with the decentralized system output shows notable differences regarding electricity balance and the number of HPs required. The central scenario's total HP electricity consumption is 17% lower than that of the decentral systems, while the central scenario has a higher electricity export and lower import from the grid with 10 and 18 percent, respectively. The central system requires two HPs for heating and three HPs for cooling, less than the decentral scenario, because the demand profiles of the buildings are not similar, and, in the central scenario, the peak demand in one building might be covered by a simultaneous valley in another building's demand. In some hours, especially peak demand hours, the total amount of surplus energy of six buildings in each hour is higher than the capacity of a single 70-ton HP (responsible for the extra number of HPs in the decentralized scenario).

Hydraulic and Thermal modeling of the Distribution Network Results
In this section, the calculated mass flow rate, heat losses, and the piping system are presented. At first, the mass flow rate for each pipe is determined. Figure 7 shows the mass flow rate of the supply line of the heating system in the network.   Figure 7 shows the mass flow rates in the main pipes, with pipe 1 and pipe 5 having the highest value. Moreover, the flow rate in pipe 6, which transfers energy from the plant to building A, has a higher value compared to the other pipe, because building A has the highest energy demand in the configuration. This is also true for the cooling system operation. After calculating the mass flow rate in the system, the size of the pipes can be determined. Table 6 shows the piping results for the heating and cooling systems. The results will be used for the economic assessment. After obtaining the flow rates and sizes of the pipes, it is possible to calculate the temperature levels in the network. This calculation helps to evaluate the heat losses in the grid and, more importantly, to make sure that the demand will be properly responded to by the system. Figure 8 shows the temperature variation of the heating system in building A's consumer substation. The temperature setpoint for heating in Montreal is 22 • C, and Figure 8 shows that this energy system scenario can properly respond to the energy demands of this building. The minimum temperature of the supply line for the heating system was 28 • C, and the maximum temperature of the supply line for the cooling system was calculated as 19 • C. Thus, the centralized scenario can properly provide energy for consumption.
The calculation showed that the heat loss for this configuration of the network is 2.15 × 10 6 kWh, which is 10.1% of the total demand.

Economic Assessment Results
In this section, the results of the economic assessment are presented. LCOE is used as an indicator for evaluating the economic performance of both scenarios. Table 7 shows the results of calculating LCOE for the first scenario, where each building has an individual geothermal loop. The electricity cost in the province of Quebec is 8 CAD/kWh. Based on the table's data, the LCOE for all the buildings is less than this value. Therefore, using the geothermal loop for buildings in the Lachine area is cheaper than using the direct electrical heating systems commonly used in Montreal. Table 7. Economic assessment of the first scenario (decentralized energy system).  In Figure 9, the left axis and the bars represent LCOE. The right axis and the circles represent building surface area. In addition, the size of circles represents building surface area to show a comparison between them.

Building Name
It is shown in the figure that the bigger the building is, the cheaper the system becomes. The LCOE for building E, the smallest building, is 0.07 CAD/kWh, while this number for building A, the biggest building in the area, is 0.05 CAD/kWh.  Table 8 shows the results of the economic assessment of the centralized system under two different sets of assumptions. In the first case, the system's detailed costs are equal to the values used for the distributed system. The cost of heat pump, labor, and borehole drilling and filling are the same as the ones used for the decentral system. In the second case, the economies of scale are considered. Hence, either a discount on different items has been received and/or the project is subsidized. The results show that, in both cases, the LCOE is lower than the electricity price in Quebec, which proves that using a geothermal district heating and cooling network is cost-beneficial when compared to direct electrical heating. Moreover, in both cases, the LCOE of the central system is lower than or equal to the LCOE in decentral system. In the unsubsidized system, only for building A, the LCOEs of both central and decentral systems are the same, while, for all other buildings, the LCOE decreased. This reduction in LCOE is more significant in smaller buildings. On the other hand, in subsidized cases, the LCOE for all buildings decreased which proves that the heating cost could become even cheaper in case of governmental support.

Conclusions
The present paper proposed a new method for district-scale automated building energy modeling and energy system simulation. This method was implemented to design an economic and efficient energy system for a Montreal case study district called the Dominion Bridge. An automated urban building energy modeling (AUBEM) method was used to calculate the energy demand of the buildings. This method carries out the whole process of 3D building modeling for the energy demand analysis and the entire process of energy system modeling in an automatic manner. Using the energy demand results from the AUBEM, an energy system model for decentral and central reversible heat pumps was prepared in the INSEL 8.2 simulation environment. Moreover, a heating network model was designed and implemented for a centralized heat pump system in MATLAB. This code provides the pipe sizing, mass flow rates, temperature distribution and heat losses in the distribution network. Finally, LCOE was used as an indicator when evaluating the economic performance of both scenarios.
Although the centralized scenario experienced heat losses through the grid, according to the ESM results, this scenario required less electricity consumption and fewer HPs to meet the demands. Moreover, the economic assessment results revealed that the LCOE of both scenarios varies from 0.04 to 0.07 CAD/kWh, which is cheaper than the electricity cost in Quebec (0.08 CAD/kWh). The LCOE for the bigger buildings was lower compared to the smaller ones. A comparison between centralized and decentralized scenarios revealed that the centralized system is cost beneficial for all buildings. It was also shown that if discounts are received due to the economy of scale, or the government subsidizes the project, the cost of heating could decrease further. If the central project is subsidized and/or discounts are received on different items, the LCOE decreases to 0.04 CAD/kWh.
The ESM results also showed that system sizing for 100% of hours (8760 in a year) results in an over-sized system for most operating hours. Adding storage, short-or long-term, and system-sizing for lower percentiles of the maximum load (i.e., P = 0.98) will lead to a better outcome. The optimal percentile determination requires a detailed investigation into load profiles, peaks and valleys, and the number of hours with consecutive high demand.