Energy, Environmental, and Economic Analyses of Geothermal Polygeneration System Using Dynamic Simulations

: This paper presents a thermodynamic, economic, and environmental analysis of a renewable polygeneration system connected to a district heating and cooling network. The system, fed by geothermal energy, provides thermal energy for heating and cooling, and domestic hot water for a residential district located in the metropolitan city of Naples (South of Italy). The produced electricity is partly used for auxiliaries of the thermal district and partly sold to the power grid. A calibration control strategy was implemented by considering manufacturer data matching the appropriate operating temperature levels in each component. The cooling and thermal demands of the connected users were calculated using suitable building dynamic simulation models. An energy network dedicated to heating and cooling loads was designed and simulated by considering the variable ground temperature throughout the year, as well as the accurate heat transfer coe ﬃ cients and pressure losses of the network pipes. The results were based on a 1-year dynamic simulation and were analyzed on a daily, monthly, and yearly basis. The performance was evaluated by means of the main economic and environmental aspects. Two parametric analyses were performed by varying geothermal well depth, to consider the uncertainty in the geoﬂuid temperature as a function of the depth, and by varying the time of operation of the district heating and cooling network. Additionally, the economic analysis was performed by considering two di ﬀ erent scenarios with and without feed-in tari ﬀ s. Based on the assumptions made, the system is economically feasible only if feed-in tari ﬀ s are considered: the minimum Simple Pay Back period is 7.00 years, corresponding to a Discounted Pay Back period of 8.84 years, and the maximum Net Present Value is 6.11 M € , corresponding to a Proﬁt Index of 77.9% and a maximum Internal Rate of Return of 13.0%. The system allows avoiding exploitation of 27.2 GWh of primary energy yearly, corresponding to 5.49 · 10 3 tons of CO 2 avoided emissions. The increase of the time of the operation increases the economic proﬁtability.


Introduction
Industrialization has promoted the use of oil, natural gas, coal, and other conventional energy sources causing the risks of stock depletion and environmental pollution [1]. Thus, the environmental emergency is a priority on the policy agenda of different countries. Indeed, in the framework of the Conference of the Parties in Paris [2], a group of countries has signed and ratified different documental for the ORC ranges from 7 to 15% with a geothermal source temperature of 160 • C. Sometimes to improve the system the integrations with solar collectors are considered to keep the working fluid under the desired conditions at the turbine inlet. In addition to trigeneration use, the geothermal brine sometimes could be used for the desalinization of water [52][53][54].
In the present study, an ORC module for geothermal applications was designed and simulated to supply the energy demand of a district in Monterusciello, a district of Pozzuoli (Naples), in the geothermal area of Phlegraean Fields, South Italy. The ORC module uses a geothermal source in cascade application to supply thermal demands for space heating and cooling and domestic hot water demands of the district in a polygeneration application; otherwise, the electric output energy is partly used to supply the energy request of thermal network and the remaining part is sold to electric national grid. The selection of organic fluids is based on a literature analysis and depends on the thermophysical, economic, and environmental properties [55,56]. An appropriate selection of the working fluid of the cycle is crucial for optimizing the efficiency of the binary plant, to maximize the conversion efficiency or to determine the best configuration for a given plant capacity. The selection of the working fluid also significantly affects the costs of the heat exchanger.
Considering the previous literature [57], the choice of the working fluid depended on the source temperature and critical temperature of the fluid. The working fluid must satisfy the general technological and environmental criteria, which have been widely discussed in the literature, for example, suitable thermodynamic fluid properties, no toxicity, no or low flammability, low cost, a low Global Warming Potential, and no Ozone Depletion Potential impact.
By using the previous analysis [58][59][60], the organic fluid R245fa was selected because of its environmental and thermodynamic efficiency in the case of heat sources with temperatures close to 150 • C for geothermal applications and/or other renewable sources for ORC applications. The use of R245fa is also effective for temperatures below 170 • C [61]. Other studies that optimize the ORC efficiency, consider the zeotropic mixtures such as the combination of 70% R245fa and 30% R125a [62].
The analyzed polygenerative plant application in Monterusciello is used to satisfy the thermal loads of a small local energy community through a mini-grid system. The geothermal district was inspired by energy community based on geothermal activity in Japanese real cases [63]. Its simulative model is created to evaluate the economic and social benefits of a district energy system as suggested in the literature [64], getting closer to meet economic and environmental needs. The study wants to encourage the energy district diffusion. It defines the limits and advantages of small independent energy districts fed by renewables reservoirs.
The keys for comprehensively developing such district contexts and making them economically profitable are the integration of different sources, the simultaneous production of different energy vectors (polygeneration), and the "load plant sharing" approach, which maximizes the source exploitation. In this context, a pre-feasibility study cannot be neglected. The authors propose a district system based on the exploitation of geothermal sources integrated with an auxiliary biomass backup system for the simultaneous production of electricity, heating, and cooling for domestic air conditioning (using a district heating and cooling network (DHCN) and domestic hot water network (DHWN). The plant system was simulated in TRNSYS environment [65], and all the heat exchangers and the ORC were first developed by previous analysis in the AspenEDR [66] and AspenONE [66] environments, respectively. AspenEDR allows a detailed and in-depth design; according to its simulation results, the heat-exchanger geometry is developed in the TRNSYS environment, creating user-defined components.
Because the ORC built-in model is absent in the TRNSYS library, it was first developed and simulated in AspenONE to extrapolate working maps as functions of the inlet heat source temperature and mass flow rate and after inserted in TRNSYS environment to create a user-defined component following the approach adopted in [52][53][54]. A control strategy is developed in this study to manage the layout and the simultaneous production of multiple energy vectors.
The dynamic model of the entire system includes a typical building configuration in the real considered residential district realized by TRNBUILD (a TRNSYS tool). The novelty of the study is related to defining a possible geothermal energy district in a not yet used geothermal field. In particular, a combination of two flexible and stable RESs is considered (biomass and geothermal energy). A sustainable energy district in load sharing configuration is analyzed by using low-temperature geothermal source, to requalify a district of Naples by improving the local energy sources and the local energy network. Moreover, a sensitivity analysis is performed by varying the depth of the geothermal well, to take into account the uncertainty related to the depth where a suitable geothermal source temperature is available. An economic analysis, which is the most central aspect of the pre-feasibility study and based on accurate market surveys, is performed by taking into account different scenarios of the Italian market for electricity production (with and without feed-in tariffs).

Area of Interest
This study should be referred to the improvement of geothermal reservoirs exploitation in Italy. Considering the general Italian geothermal context, it has been highlighted from the geothermal maps (available at depths of 1000, 2000, and 3000 m) that the areas of interest concerning the temperature and geothermal flow are located in Tuscany, Aeolian Islands, and Neapolitan area [67,68]. In [69,70], it is demonstrated that the geothermal anomalies also in different Italian areas such as the Alps, Sicily, and the central Tyrrhenian and the Mediterranean Sea with a heat flux of 80, 40, and higher than 150 mW/m 2 , respectively.
The area of interest in this study is the Neapolitan area (Phlegraean Fields), which (similar to numerous wells already surveyed in the 1980s) has a highly aggressive geothermal fluid and temperatures of approximately 100 ÷ 150 • C at the land surface [71][72][73]. Interest in the geothermal area of Naples started in 1930 and grew rapidly in the mid-1980s. A total of 117 geothermal wells were investigated, with a maximum depth of 3046 m. The results were particularly encouraging for the Phlegraean Fields [74] and Ischia, where high geothermal gradients have been recorded, owing to the presence of localized high enthalpy fluids (T > 150 • C) at low depths (hundreds of meters), and both steam and water dominated [75][76][77].
In the Neapolitan contest, the Phlegraean Fields represent the area of greatest interest. According to the previous investigation for some reservoirs of these sites, the average geothermal gradient is 0.170 • C/m (average reference value of 0.03 • C/m) [68] and the average geothermal flow is 149 mW/m 2 (average reference value of 63.0 mW/m 2 ) [68]. In this area, buildings with different intended uses (such as residential, commercial, and also industrial buildings) can be found and they all can be connected by a thermal grid fed by geothermal sources.
In the past, the local geothermal applications in the Neapolitan area provided to direct use of the source in thermal employment using available reservoir sites in medium-low enthalpy. Nowadays, the technological maturity of the ORC component and the large diffusion of polygeneration systems to define energy district can be made usable these low/medium enthalpy sources. This study analyses the possibility to employ the large low/medium temperature geothermal sources available in Phlegraean Fields in polygeneration approach to satisfy the energy loads of a district. In addition, this work focuses on the replacement of pre-existent wells realized from 1930 until 1985, during which a large amount of data were collected firstly by the SAFEN Company and successively in ENEL-AGIP Joint Venture for geothermal exploration [77].
In a previous study about Monterusciello [78], it was analyzed only a system, without ORC module, meeting thermal loads using geothermal fluid at a temperature of 50 • C available at a depth of 100 m. In the current analysis, the geothermal fluid could be obtained from two wells of extraction (really investigated from the previous AGIP campaigns [77]) at the desired temperature and flow rate. In this upgrade study, an analysis of the existing geothermal wells of Phlegraean Fields was conducted as reported in Table 1. According to geological maps, the geothermal analyzed wells are collected and geolocated. For each well, the geothermal gradient is estimated by considering a linearity approximation as reported in the last column of the table. For Monterusciello, the gradient considered Energies 2020, 13, 4603 6 of 34 is 0.1 • C/m in the simulation at depth of 1500 m that guarantees a geothermal brine temperature of 150 • C. This value represents a possible real gradient in the Monterusciello by considering the near wells.

Buildings and Heating and Cooling Network Characterization
A previous analysis for the building modeling was performed based on a field investigation to define both for the envelope characteristics and users/consumers typologies. The definition and modeling of the building are set up using a tool of TRNSYS software (TRNBUILD). Each building of the district, which is an ensemble of social housing, consists of four floors and eight apartments and is in a residential zone. In the SKETCHUP environment [79], the rendering of the building is represented in Figure 1. Each building consists of four floors with eight thermal zones of 198 m 2 and two apartments for each thermal zone; such a configuration was implemented by considering, on one hand, a plausible reproducibility of the considered context and, on the other hand, the computational burden of such a complex layout. The building model reproduces the typical council housing of 1960s. The stratigraphy respects the Italian building regulations for its specific age of construction [80]. All 90 buildings of the district led to a total of 1440 apartments and 1980 habitants, covering~6% of the Monterusciello population. The opaque and transparent building envelope characteristics and the geometrical parameters for each building are reported in Table 2. The transmittance of windows is referred to a single glass component with aluminum frame. Space heating and cooling loads were evaluated considering an occupancy schedule, represented in Figure 2, based on four different cases (a, b, c, and d): • case a: two working adults and two children, • case b: two adults, a child, and an elderly person, • case c: two elderly habitants, and • case d: two working adults.   While during the winter period the heating load depends on the desired comfort condition (indoor temperature at 20 °C) required in the whole apartment, during the summer season the energy delivered by cooling network depends by activation of terminal cooling units that commonly serve only the rooms of the apartment effectively occupied. The yearly trends of heating and cooling load for the whole building are presented in Figure 3. The 1-year simulation was based on the weather file available in Trnsys library, which refers to [83].     While during the winter period the heating load depends on the desired comfort condition (indoor temperature at 20 °C) required in the whole apartment, during the summer season the energy delivered by cooling network depends by activation of terminal cooling units that commonly serve only the rooms of the apartment effectively occupied. The yearly trends of heating and cooling load for the whole building are presented in Figure 3. The 1-year simulation was based on the weather file available in Trnsys library, which refers to [83].  Each one of these schedules refers to two apartments of each considered building. For the domestic hot water, load is defined as a normalized hourly profile for different months of the year [81] and a daily volume consumption of 50 L per district resident [79].
According to Italian normative [82], the heating period for Naples in climate zone C goes from 15 November to 31 March. While the cooling period is from 1 June to 30 September. The maximum number of heating and cooling systems operating hours is established to 10 h per day.
All the electric energy produced by the plant is sold to the power grid excluding the part needed for auxiliaries of the thermal district network.
While during the winter period the heating load depends on the desired comfort condition (indoor temperature at 20 • C) required in the whole apartment, during the summer season the energy delivered by cooling network depends by activation of terminal cooling units that commonly serve only the rooms of the apartment effectively occupied. The yearly trends of heating and cooling load for the whole building are presented in Figure 3. The 1-year simulation was based on the weather file available in Trnsys library, which refers to [83].  The buildings of the Monterusciello District area, modeled through TRNBUILD tools and based on the previous information, are connected by means of two energy networks: DHCN and DHWN.
The district area is divided by considering the constraints of the DHCN and DHWN network installation. Figure 4 reports the pipelines segments (in light blue color) running along the main roads, the spatial distribution of the buildings and a viable energy conversion system plant location (in red color). The pipelines segment are labeled as "A" or "B", referring to different blocks of the district or "MAIN". The lengths of the pipeline are consequently calculated by considering the real building location and the plausible location of the plant. The diameters sizing takes into account the thermal demand distribution and a fluid velocity of 2.0 m/s is considered. The DHWN sizing was performed similarly to the DHCN one. The geometric features of piping are listed in Table A1 in Appendix. The DHCN considers a time operation of ten hours (from 10 a.m. to 8 p.m.). A parametrical analysis was also conducted to consider the cases in which DHCN is turned off 2 and 4 h later than the base case (10 p.m. and 12 p.m.). The DHW is guaranteed for all 24 h. The set point of DHW temperature is ensured by the biomass boiler. The overall head and thermal losses of the distribution networks takes into account the ones related to the DHCN and the one related to the DHWN.  The buildings of the Monterusciello District area, modeled through TRNBUILD tools and based on the previous information, are connected by means of two energy networks: DHCN and DHWN.
The district area is divided by considering the constraints of the DHCN and DHWN network installation. Figure 4 reports the pipelines segments (in light blue color) running along the main roads, the spatial distribution of the buildings and a viable energy conversion system plant location (in red color). The pipelines segment are labeled as "A" or "B", referring to different blocks of the district or "MAIN". The lengths of the pipeline are consequently calculated by considering the real building location and the plausible location of the plant. The diameters sizing takes into account the thermal demand distribution and a fluid velocity of 2.0 m/s is considered. The DHWN sizing was performed similarly to the DHCN one. The geometric features of piping are listed in Table A1 in Appendix. The DHCN considers a time operation of ten hours (from 10 a.m. to 8 p.m.). A parametrical analysis was also conducted to consider the cases in which DHCN is turned off 2 and 4 h later than the base case (10 p.m. and 12 p.m.). The DHW is guaranteed for all 24 h. The set point of DHW temperature is ensured by the biomass boiler. The overall head and thermal losses of the distribution networks takes into account the ones related to the DHCN and the one related to the DHWN.  The buildings of the Monterusciello District area, modeled through TRNBUILD tools and based on the previous information, are connected by means of two energy networks: DHCN and DHWN.
The district area is divided by considering the constraints of the DHCN and DHWN network installation. Figure 4 reports the pipelines segments (in light blue color) running along the main roads, the spatial distribution of the buildings and a viable energy conversion system plant location (in red color). The pipelines segment are labeled as "A" or "B", referring to different blocks of the district or "MAIN". The lengths of the pipeline are consequently calculated by considering the real building location and the plausible location of the plant. The diameters sizing takes into account the thermal demand distribution and a fluid velocity of 2.0 m/s is considered. The DHWN sizing was performed similarly to the DHCN one. The geometric features of piping are listed in Table A1 in Appendix. The DHCN considers a time operation of ten hours (from 10 a.m. to 8 p.m.). A parametrical analysis was also conducted to consider the cases in which DHCN is turned off 2 and 4 h later than the base case (10 p.m. and 12 p.m.). The DHW is guaranteed for all 24 h. The set point of DHW temperature is ensured by the biomass boiler. The overall head and thermal losses of the distribution networks takes into account the ones related to the DHCN and the one related to the DHWN.

System Configuration and Layout
The simplified system layout is shown in Figure 5, where the main components are presented but the DHCN and the DHWN are omitted. For simplicity of the scheme, only one production well is represented, even if such an analysis two production wells are supposed to be exploited.

temperature.
GHE2 is designed to provide a suitable HW mass flow rate basing on the temperature approach. If the temperature of the HW is higher than that of the geothermal source entering GHE2, the thermostatic valve (through D1/M1) diverts the HW bypassing GHE2, to prevent heat dissipation. Diverter D2 diverts the HW flow to keep TK1 and TK3 thermally loaded, giving priority to TK1. Different opportune set temperatures on the top of the tanks are provided for both the heating and cooling operation modes. A thermostatic valve (through D3/M3) ensures a constant temperature (45 °C) of the DHW for the networks.
A thermostatic valve (through D4/M4) ensures a constant temperature of the HWDHCN for the DHCN, which is sent directly to the network during the heating operation while it is sent to the First, geofluid powers an ORC module through a heat exchanger (GHE1) heating the hot water (HW ORC ); the HW ORC heats the ORC working fluid in the evaporator.
The ORC is calibrated to ensure a constant power production; it is condensed with cooling water (CW), which is cooled through the cooling tower (CT ORC ) and whose mass flow rate is adjusted using a variable speed pump to keep the temperature difference at the condenser constant. The ORC module is intended to operate under steady-state conditions with stable evaporation and condensation pressures. This condition at the evaporator is guaranteed by fixed temperature and flow rate of geothermal brine. While to assure the stationary condition at the condenser a variable cooling water flow rate from the water supply net is considered.
A little part of the electricity available from the ORC plant is used to cover the network auxiliaries (pumps, etc.), while the most is sold to the national grid. After delivering thermal energy at ORC evaporator the geothermal fluid is used to provide space heating, cooling, and DHW by means of DHCN and DHWN. A brief description of the system operation and the control strategy is presented herein following a previous study [78].
The Hot Water (HW, red line in the layout) feeding tanks (TK1 and TK3) are produced through the GHE2 heat exchanger. One part of HW is stored in a tank (TK1) to recover thermal energy for space heating and cooling. Through a control strategy, HW is alternatively sent to storage tank TK3, and it is used to store and heat the DHW. Then, a dedicated subsystem is used to ensure a constant temperature.
GHE2 is designed to provide a suitable HW mass flow rate basing on the temperature approach. If the temperature of the HW is higher than that of the geothermal source entering GHE2, the thermostatic valve (through D1/M1) diverts the HW bypassing GHE2, to prevent heat dissipation.
Diverter D2 diverts the HW flow to keep TK1 and TK3 thermally loaded, giving priority to TK1. Different opportune set temperatures on the top of the tanks are provided for both the heating and cooling operation modes. A thermostatic valve (through D3/M3) ensures a constant temperature (45 • C) of the DHW for the networks.
A thermostatic valve (through D4/M4) ensures a constant temperature of the HW DHCN for the DHCN, which is sent directly to the network during the heating operation while it is sent to the adsorption chiller (ADS) during the cooling operation. The HW temperature is set to 50.0 • C during the heating operation and 100 • C during the cooling operation for feeding the ADS. If the set point temperatures cannot be ensured, the biomass boiler (BB) fed with wood chip, is activated until the indicated temperature is reached and until the thermal storage TK1 is thermally loaded. Similarly, if the DHWN set point temperature cannot be ensured, the DHW flow is diverted to the biomass boiler trough thermostatic valve D7/M7. During the cooling operation, the ADS is activated (by cooling tower, CT ADS ), and the chilled water (ChW ADS ) produced is sent to storage tank TK2, where the water ChW DHCN is stored at 10.0 • C and is sent to the DCHN.
To fed DHCN and DHWN the GHE2 is designed to provide a suitable mass flow rate of HW in a temperature range of 40 ÷ 75 • C.
The modeling of the plant has considered an ad hoc defined flow rate for geothermal fluid and a recommended minimum rejection temperature of 70.0 • C as the literature suggestions [84,85]. This value represents a good compromise for avoiding excessive depletion of the geothermal source and moderately mitigating problems involving scaling at the heat exchangers and the mechanical apparatus of the rejection process. Once the geothermal source is set, according to the mass flow rate and temperature range, the available thermal energy is calculated, and the entire system is calibrated accordingly. The ORC is calibrated to ensure constant power production. It is intended to operate under steady-state conditions of working fluid (R245fa); the evaporation and condensation pressures are stable.

Models
The development of the whole simulation model has been carried out in three steps, namely, the development of the ORC module in AspenPlus environment, the heat exchanger in AspenEDR environment, and the dynamic simulation model in TRNSYS environment. The TRNSYS library does not include the ORC module; the ORC was implemented in AspenPlus software and simulated by varying the inlet temperature of the source; there were obtained working maps reporting the main output parameters as a function of the inlet temperature of the thermal source, given a constant mass flow rate.
These functions have been implemented in TRNSYS environment adding a user-defined component. The entire dynamic simulation model also includes the ORC heat exchangers modules (evaporator and condenser) and geothermal heat exchangers (GHE1 and GHE2) designed, as previously defined in the introduction section, in the AspenPlus and AspenEDR environments, respectively.
All the heat exchangers of the plant have been developed and designed in AspenEDR environment. After this, the geometry and suitable heat transfer coefficient correlations have been implemented in TRNSYS software to create a user-defined component: this approach is more rigorous with respect to using the built-in ones since it takes into account the real instant operation of the heat exchanger in terms of overall heat transfer coefficient and efficiency. In the Table A2 in Appendix A all the input parameters of the dynamic model are listed for ORC, working fluid, and heat exchangers. In particular, the heat exchanger GHE1 and GHE2 are simulated considering shell and tube heat exchanger model in titanium material. The resulted characteristics of GHE1 and GHE2 show an external diameter of tubes of 19.05 mm, thickness of tube of 1.2 mm, pitch of 23.8 mm. For the ORC module, the isentropic efficiencies are fixed: 70% and 85% for the pump and turbine, respectively. The ORC pinch point temperature differences are 7 and 5 • C respectively at the evaporator and the condenser. The working fluid used for the ORC plant is the R245fa. The evaporator and condenser for ORC module are two AISI306 shell and tube heat exchangers. The output parameters of the heat-exchanger design and the ORC module obtained as results are presented in Tables A3 and A4, respectively. In the ORC evaporator (not showed in the layout figure), the inlet and outlet temperatures of R245fa are, respectively, 51.08 and 120.1 • C at 19.35 bar, while at the condenser, they are 66.05 • C for the inlet and 50 • C in the outlet section at 3.43 bar. The thermal power exchanged in the evaporator and condenser is 4698 and 3835 kW, respectively.
The nomenclature used to define each parameter in Tables A2-A4 is referred to system layout scheme represented in Figure 5.
The ORC is calibrated considering a nominal power production of 500 kW el . All the heat exchangers are simulated by calculating the outlet streams' conditions using the heat and mass balances and the surface area. The inlet and the outlet temperatures of the condensation process are fixed. The mass flow rate is determined accordingly.
Once the cycle is completely defined and simulated in the AspenONE environment, the evaporator and the condenser are designed in AspenEDR.
Both the geothermal heat exchangers (GHE1 and GHE2) are designed by considering plausible values from an ORC market survey of the cold-side mass flow rate for feeding the ORC module, a plausible value of the cold-side temperature difference, and a hot inlet-cold outlet temperature approach. The correlations adopted for the overall heat-transfer coefficient to define heat exchangers parameters are based on the fully developed laminar flow inside the duct with an isothermal wall [86] and the fully developed turbulent flow inside the duct with an isothermal wall [87].
The calibration of the overall system is based on a preventive dynamic simulation of the building equipped with a plausible number of fan coils per apartment, which are also implemented in the TRNSYS environment. Then, according to the simulation results, the heat-exchanger geometry is developed in the TRNSYS environment (with the creation of user-defined components using macros and Calculator blocks) to dynamically simulate the real heat-transfer performance with regard to the instant value of the overall heat-transfer coefficient, efficiency, and number thermal Unit (NTU).
The simulation model of the building is linked to the simulation model of the fan coil in the same TRNSYS environment. The calibration of the fan-coil simulation model is based on a market survey.
The seasonal data about fan coil are listed in the Table 3, such as the nominal power of fan coil (P nom,fc ), flow rate water and air of fan coil (m fc,w , m fc,a ), the comfort temperature set (T set,amb ), and the ratio between sensible and total thermal power (P th,sens /P th,TOT ). This approach allows us to determine the effective thermal and cooling loads under dynamic operating conditions and then take into account variable weather conditions, the thermal inertia of the building envelope, and the variability of the space occupation.
Once the time-dependent thermal and cooling loads are extrapolated and the maximum thermal energy available from the geothermal source is considered, it is possible to determine the maximum number of apartments served by the system.
Finally, the variable parameters are calibrated in the TRNSYS environment such as mass flow rates, set point temperatures, and characteristics of components under nominal conditions.
The system plant simulation model is linked to the DHCN one, which takes into account the heat losses occurring in the network piping. In turn, the DHCN simulation model is linked to the thermal and cooling loads, allowing the return water conditions to be determined.

Methodologies
The system described in the previous sections has been analyzed from an economic an environmental point of view according to the following methodology.
As regards to economic analysis, the total plant investment cost, Z TOT , is given by the sum of the costs of all the modules composing the plant, the cost of the distribution networks, and the BOP (Balance of the Plant) cost Z BOP , to take into account all the auxiliary systems and supporting components of the plant. Table 4 reports the parameters used in the economic analysis, as well as the main parameters of the reference scenarios for each energy vector and emission factor. The performances of the renewable ORC coupled with DHCN system are evaluated by assessing the energetic and economic savings, and by comparing the proposed system (PS) with the reference one (RS). In particular, in case of the reference scenario • a reversible electric-driven air-to-water heat pump for space heating and cooling is supposed to be installed in each building; the coefficient of performance (COP) in heating mode is supposed to be 3.0, while in cooling mode it is supposed to be 2.5 [88]; • a natural gas boiler for DHW production is supposed to be installed in each dwelling; its efficiency is supposed to be 0.85; • power production is given by the national grid characterized by an overall average efficiency, η grid , equal to 46% [89] The same space heating and cooling terminal units (fan coils) are adopted both in case of PS and RS.
All the electricity prices (for the cases with and without feed-in tariffs) are based on the Italian real market trends from January 2018 to August 2019 [31]. Different prices are provided for the three time-dependent partitions of the Italian market (F1, F2, and F3). The prices of heating and cooling energy are based on market surveys.
The Z BOP is calculated as 3.00% of the overall plant cost. An economic analysis is performed by taking into account two different scenarios for electricity sales: with and without feed-in tariffs provided by the Italian market [31].
In the case of feed-in tariffs, the revenue related to electricity sales R el,feed-in is defined as follows, where E el represents the amount of yearly electricity produced, and c el,feed-in represents the price of electricity under feed-in tariff market conditions. In the case where no feed-in tariffs are provided, a guarantee minimum price tariff system is assumed: regardless of the effective electricity market price, the plant always sells energy above a certain threshold price established by the GSE (Gestore dei Servizi Energetici, i.e., Italian energy services management institution). Moreover, three selling prices are considered based on the hourly partition of the time-dependent tariff system. Given the uncertainty of the results due to the variability of prices, a plausible range of the revenue (i.e., R el,min and R el,max ) is calculated according to the minimum and maximum values of the sale price. The minimum and maximum sale price values are established with consideration of the real price trend from January 2018 to August 2019 (reported in Figure A1 in Appendix A).
Then, the minimum and the maximum revenues are calculated as follows, R el,min = E el f F1 c el,F1,min + E el f F2 c el,F2,min + E el f F3 c el,F3,min where the subscript F i indicates the specific partition of the time-dependent tariff system and the subscript f Fi indicates the yearly fraction of the hours belonging to the specific F i partition.
The revenue related to the selling of thermal energy for space heating R th is calculated as follows, where E th represents the amount of thermal energy provided to the network yearly, and c th represents the price of thermal energy based on an Italian market survey. Similarly, the revenue related to the selling of cooling energy R cool is defined as where E cool represents the amount of cooling energy provided to the network yearly, and c cool represents the price of cooling energy based on an Italian market survey. The simple payback (SPB) is defined as follows, SPB = Z TOT CFS (6) where CFS represents the yearly cash-flow statement, which is defined as R TOT represents the sum of all the revenues, and C TOT represents the sum of the yearly maintenance cost C O&M and the yearly operational cost C Op , which are defined as follows, where P AUX,TOT represents the total thermal energy provided by the auxiliary boiler, LHV biom represents the lower heating value of the biomass (wood chips), and c biom represents the unit cost of the biomass. Finally, the main economic indicators are calculated to assess the economic profitability of the system, i.e., the discounted payback period (DPB), net present value (NPV), profit index (PI), and internal rate of return (IRR): where a represents the discount rate, AF represents the annuity factor, and N represents the service life.
Regarding the environmental analysis, the saved primary energy source (PE) is calculated by considering a specific reference scenario for each produced energy vector (electric energy and thermal energy for air conditioning).
Then, the PE is calculated as follows, where  (15) where EF grid represents the emission factor related to the national grid and EF nat,gas represents the emission factor of the natural gas. The implemented economic model calculates the investment and the operating costs of both PS and RS. The cost functions adopted for PS and RS components are taken from scientific and technical literature. The cost of GHE1 and GHE2 are obtained by previous literature studies as a function of the heat exchanger area [51]. The tank cost is a function of occupied volume [90]. The cost of line networks for DHCN and DHWN is the functions of diameters and length of pipes [91]. The cost of pumping depends on flow rates which cross the pumps [92]. The study realizes a parametric analysis with different depths of the geothermal well. The depth of the geothermal well affects both the operation of the system plant and economic performance. The aim of the analysis is the prefeasibility study by considering the uncertainty related to the depth where the geothermal source is available and assessing a range for each main performance indicator where the system can be considered feasible and profitable.
A second parametric analysis is performed with different operation times of the DHCN, without changing the operation time of the DHWN; this ensured a constant temperature of 45 • C. When the DHCN operation time is reduced, a larger amount of electric energy is available, because all the auxiliaries of the DHCN subsystem are off. Consequently, higher revenue related to electricity production is expected.
Additionally, the DHWN subsystem is forced to use the auxiliary biomass boiler to ensure the appropriate temperature of the DHWN; thus, a higher operational cost related to boiler fuel (wood chip) is expected.
The objective is to evaluate the effects of the production strategy on profitability and to determine the optimum operation point.

Thermodynamic Analysis
First, the daily results are presented and discussed through two days representative of the system functioning for each operation mode: one for heating (winter day) and one for cooling (summer day). Figures 6 and 7 show the main system outputs, respectively, for the winter day where the nomenclature and numeration refer to Figure 5.  When the DHCN is on, ORC power is employed to feed all the system auxiliaries (pumps and the overall control and monitoring system). Then, the net power sent to the grid decreases to 380 kWel. Figures 8 and 9 show the main system temperatures for the characteristic representative summer day, and Figure 10 presents the main system powers.
In particular, Figure 9 indicates the temperature trend of TK3 (T15, T16, T19, T20) and the ADS module (T13, T14, T15, T16, T17, T18). As shown, the feed temperature T13 of the ADS module is constant at 100 °C, except at the initial time, where no cooling load is observed ( Figure 10) and the feed water is available at 120 °C. Moreover, the results indicate that the system is perfectly capable of ensuring a constant temperature T19 (10.0 °C) of the cooled water to be sent to the network, which corresponds (taking into account the losses in the network piping) to the set temperature of the fan coil in the cooling mode.  When the DHCN is off, T 3 coincides with T 2 , whereas it continuously varies depending on the inlet condition of the cool side at GHE2, i.e., T 10 . T 3 never decreases to 70.0 • C, due to the calibration of the overall system and the control strategy. T 10 and T 8 exhibit the trends of TK1 to and from GHE2, respectively; T 11 and T 21 represent the trends of the water temperatures to and from the DHCN, respectively. The system can ensure a constant temperature of 50.0 • C to the DHCN and the set point temperature of the fan coils installed in the apartments taking into account the small heat losses occurring in the piping. Moreover, the system is capable to give the DHW at a constant temperature value of 45.0 • C, as shown in Figure 6, independently by the load, shown in Figure 7.
From 6:00 to 10:00, the geothermal source is mainly employed to thermally load TK3, which belongs to the DHW subsystem. From 10:00 to 20:00, the source is mainly employed to thermally load TK1, which belongs to the DHCN.
Simultaneously, the geothermal source is constantly employed to feed the ORC module.
The foregoing is clearly shown in Figure 7, which presents the main system powers. When the DHCN is off, the amount of achievable thermal power by the geothermal fluid is equivalent approximately to 4.50 MW, whereas it increases to 13.0 MW when the load required by the network is maximized.
The ORC net power production depends on the operation time of the DHCN. When the DHCN is off, the ORC power (approximately 430 kW el ) is employed for geothermal fluid pumping and the ORC cooling tower auxiliaries, and the remaining part is sent to the grid.
When the DHCN is on, ORC power is employed to feed all the system auxiliaries (pumps and the overall control and monitoring system). Then, the net power sent to the grid decreases to 380 kW el . Figures 8 and 9 show the main system temperatures for the characteristic representative summer day, and Figure 10 presents the main system powers. When the DHCN is on, ORC power is employed to feed all the system auxiliaries (pumps and the overall control and monitoring system). Then, the net power sent to the grid decreases to 380 kWel. Figures 8 and 9 show the main system temperatures for the characteristic representative summer day, and Figure 10 presents the main system powers.
In particular, Figure 9 indicates the temperature trend of TK3 (T15, T16, T19, T20) and the ADS module (T13, T14, T15, T16, T17, T18). As shown, the feed temperature T13 of the ADS module is constant at 100 °C, except at the initial time, where no cooling load is observed ( Figure 10) and the feed water is available at 120 °C. Moreover, the results indicate that the system is perfectly capable of ensuring a constant temperature T19 (10.0 °C) of the cooled water to be sent to the network, which corresponds (taking into account the losses in the network piping) to the set temperature of the fan coil in the cooling mode.  Table 5 presents the main results of the thermodynamic analysis on a monthly and yearly basis. As shown, the system is perfectly capable of ensuring the temperature levels in each subsystem and it is perfectly capable of covering the total network load required (sum of thermal total power required by DHCN and DHWN) and the losses occurring in the network. During the cooling season, the overall cooling energy provided by the system for district cooling, ESUPPLIED,DHCN is higher than the cooling energy required Ereq.SUMM,DHCN; even in this case the difference is caused by the losses occurring at the network piping.
Similarly to the heating season, during the cooling season the cooling energy exiting the TK2 ETK2,out is higher than the supplied one ESUPPLIED,DHCN because of the losses occurring at the system plant.     In particular, Figure 9 indicates the temperature trend of TK3 (T 15 , T 16 , T 19 , T 20 ) and the ADS  module (T 13 , T 14 , T 15 , T 16 , T 17 , T 18 ). As shown, the feed temperature T 13 of the ADS module is constant at 100 • C, except at the initial time, where no cooling load is observed ( Figure 10) and the feed water is available at 120 • C. Moreover, the results indicate that the system is perfectly capable of ensuring a constant temperature T 19 (10.0 • C) of the cooled water to be sent to the network, which corresponds (taking into account the losses in the network piping) to the set temperature of the fan coil in the cooling mode. Table 5 presents the main results of the thermodynamic analysis on a monthly and yearly basis. As shown, the system is perfectly capable of ensuring the temperature levels in each subsystem and it is perfectly capable of covering the total network load required (sum of thermal total power required by DHCN and DHWN) and the losses occurring in the network.
During the heating season, the overall thermal energy provided by the system for district heating, E SUPPLIED,DHCN , is higher than the thermal energy required, E req,WINT,DHCN ; the difference is caused by the losses occurring at the network piping.
Moreover, the sum given by the thermal energy exiting the TK1, E TK1,out , and the one coming from the auxiliary boiler for the DHCN, E AUX,DHCN , is higher than E SUPPLIED,DHCN : this is due to the thermal losses occurring at the plant piping.
Thermal losses occurring at the thermal storages are given by the difference between the thermal energy entering and exiting the tanks.
What has just been discussed can be referred to the domestic hot water production system. During the cooling season, the overall cooling energy provided by the system for district cooling, E SUPPLIED,DHCN is higher than the cooling energy required E req.SUMM,DHCN ; even in this case the difference is caused by the losses occurring at the network piping.
Similarly to the heating season, during the cooling season the cooling energy exiting the TK2 E TK2,out is higher than the supplied one E SUPPLIED,DHCN because of the losses occurring at the system plant. The gross amount of electricity produced, P el,ORC , remains constant throughout the year, and P el,net , which represents the amount of electricity produced and sold to the grid, is lower during the cooling operation than during the heating operation. This is because the entire ADS subsystem, which is off during the heating operation, must be powered.
The constant average first-law efficiency of the ORC is 10.6%, owing to the steady working condition of the module.
The ADS, whose COP is 0.700 under nominal conditions (temperature of the feed hot water equal to 100 • C, the temperature of chilled water equal to 7 • C, and temperature of the cooling water equal to 22 • C), operates during the year with an average COP of 0.660: this is due to fluctuations in operating conditions (in terms of the feed water temperature, cooling water temperature, load), which move the module from the optimum point of operation.

Environmental and Economic Data Analysis
In this section, the main results of the economic and environmental analyses are presented and briefly discussed. The costs of the components are presented in Table 6, along with the cost function adopted. The total investment cost takes into account three different values of the depth of the geothermal well: the base case at 1500 m and the case of higher depth at 2000 and 2500 m by considering two production and one reinjection wells. Globally, the total investment cost (Z TOT ) ranges between 7.84 and 8.21 M€.  Table 7 presents the detailed yearly economic analysis results for each well depth and the different scenarios of the electricity market (with and without feed-in tariffs).
In the case of no feed-in tariffs, the total revenue R TOT ranges from 1.09 to 1.18 M€. With a total operational and maintenance cost C TOT of approximately 0.40 M€, the CFS ranges between 0.69 and 0.78 M€ per year.
All the economic indicators are negatively affected by the low CFS, and the system appears to be unprofitable. Table 7. Main results of the economic analysis (on a yearly basis).

R el (M€) R th (M€) R cool (M€) R TOT (M€) C O&M (M€) C Op (M€) C TOT (M€) CFS (M€)
No Feed-in Tariffs   deep  min  max  min  max  min  The SPB period exhibits a minimum value of 10.0 years and a maximum value of 11.9 years. These are far higher than the acceptable range of 5.00 to 7.00 years for private investments.
With a discount rate of 5.00%, the DPB period increases (minimum value of 14.3 years and maximum value of 18.5 years).
The NPV exhibits low values ranging between 0.388 and 1.88 M€, corresponding to a PI ranging between 4.73% and 24.0%. Acceptable values of the PI are between 60.0% and 70.0%.
Finally, the IRR exhibits low values between 5.51% and 7.64% The economic profitability is improved in the case of the feed-in tariffs for electricity selling. In fact, given the constant revenues related to the thermal and cooling energy for air conditioning, R th , and R cool , as well as the constant operational and maintenance costs C TOT , in correspondence of CFS value of 1.12 M€ the higher revenue related to power production R el are equal to 0.631 M€.
Consequently, acceptable values of the SPB and DPB are obtained. The SPB period ranges between 7.00 and 7.33 years, and the DPB period ranges between 8.84 and 9.36 years.
The NPV is increased to values between 5.74 and 6.11 M€, corresponding to PIs of 69.9% and 77.9%, respectively.
Finally, the IRR ranges between 12.2% and 13.0%. Globally, as shown in Table 8, such a renewable system allows avoiding the employment of 27.2 GWh of primary energy and (depending on the reference scenario adopted for the energy and material outputs) avoiding 5.49 × 10 3 tons of CO 2 emissions. The main results of the thermodynamic and economic analyses performed by varying the operation time of the DHCN are presented from Tables A5-A9 in Appendix A. In detail, the data about electricity production (R el ), the CFS, the SPB, the NPV, and the PIs in the case of feed-in tariffs (green) and the cases of the minimum and maximum selling prices of electricity are reported.
Longer operation time of the DHCN subsystem corresponds to higher economic profitability: all the economic indices are affected by lower DHCN operation time.
In fact, in the case of the longest time of operation, a slightly lower amount of electricity is sold to the grid because the auxiliaries are on, and the revenue related to electricity is lower as well.
On the other hand, a higher amount of thermal energy is available from the geothermal source, then avoiding the use of biomass and reducing the operational cost of the biomass boiler to supply the DHW.
In case of the shortest time of operation, a higher amount of electricity is available (auxiliaries are off) for the selling. At the same time, any load from DHW network is cover by the biomass boiler, and a higher operational cost related to boiler fuel is obtained.
The higher operational cost of the biomass boiler prevails on the higher revenue related to electricity when the time of operation is reduced and the global profitability decreases.
The revenues related to space heating and cooling are constant. It is worth noticing that, despite the analysis was carried out with the greatest possible precision, the results strictly depend on the specific case study analyzed and on the considered market context; moreover, they are affected by unavoidable uncertainty given by weather conditions, which affect the electric and thermal loads, and the ORC operation (including the cooling tower); -real availability of geothermal source, which can be assessed only after a specific exploration campaign; -fluctuations of the electricity prices and variation of the feed-in tariffs system; and variations of prices related to district heating and cooling (even though they are characterized by more stability if compared to electricity one).

Conclusions
A thermodynamic, economic, and environmental analysis of a renewable polygeneration system connected to a DHCN was performed.
The system is designed for a suburban area of the metropolitan city of Naples (South of Italy) and it is powered by geothermal sources and biomass, producing electricity, thermal energy for space heating and cooling, and DHW production.
The entire dynamic simulation model (which includes an organic Rankine cycle module, an ADS, an auxiliary biomass boiler, geothermal heat exchangers, thermal storage tanks, the DHCN, residential systems of space conditioning including terminals, and a suitable model of the building envelope) was developed and implemented in AspenONE and TRNSYS environments.
The layout and the control strategy were implemented to match the appropriate operating temperature levels for each component and to prevent the temperature of the geothermal fluid reinjected into the wells from decreasing below 70.0 • C.
The overall component calibration and the economic analysis were based on manufacturer data and market surveys.
The economical returns are estimated by considering the electric power sold to the national grid. Differently, the thermal energy for heating and cooling is employed to satisfy the thermal yearly load of the whole district.
An analysis was performed with different depths of geothermal wells where the source is available considering a geothermal temperature gradient of 0.1 • C/m. An economic analysis was performed for two different electricity purchase scenarios: with and without feed-in tariffs.
The system, whose investment cost ranges between 7.84 and 8.02 M€, is economically feasible only under feed-in tariff conditions.
In fact, without feed-in tariffs, economic indicators suggest that the system is unprofitable: the minimum SPB period is 10.1 years (corresponding to a DPB period of 14.3 years if a discount rate of 5.00% is applied), and the maximum NPV is 1.88 M€ (corresponding to a maximum PI of 24.0%).
Finally, the maximum IRR is 7.64%. Conversely, if feed-in tariffs are considered, the economic indicators suggest that for the specific case study analyzed, the system is attractive.
The minimum SPB period is 7.00 years, corresponding to a DPB period of 8.84 years. The maximum NPV is increased to 6.11 M€, corresponding to a PI of 77.9%. The maximum IRR is 13.0%.
The system allows avoiding exploitation of 27.2 GWh of primary energy yearly, corresponding to 5.49 × 10 3 tons of CO 2 emissions avoided yearly.
The keys to making the small-medium scale systems powered with geothermal sources feasible and attractive are the polygeneration and the load-sharing, which aim at maximizing the source exploitation by diversifying the energy and material vectors produced and by obtaining loads and requests that are distributed over time to the greatest extent possible.
In future studies, different typologies of final users and different energy vectors could be considered using the sources availability 24 h/24 h to feed the plant.

Acknowledgments:
The authors gratefully acknowledge the financial support of provided through the GeoGrid project POR Campania FESR 2014/2020 CUP B43D18000230007. Nicola Massarotti gratefully acknowledges the local program of the University of Napoli "Parthenope" for the support of individual research.

Appendix A.2. Cost Functions
In the following are the cost functions used in the economic analysis; it is also indicated the reference the function was taken from.
The cost of the well zwell, which is divided into three different segments, is a function of the diameter and depth of each segment. In this work, two production well and one reinjection well  (F1, F2, and F3)