Sustainability, climate resiliency, and mitigation capacity of geothermal heat pump systems in cold regions

Geothermal heat pump (GHP) systems have recently gained popularity in providing heat to buildings. In this study, the short-term and long-term performance of closed-loop horizontal GHP systems in cold regions and the eﬀectiveness of seasonal balancing were investigated. A period of 50 years was simulated to cover the entire life of a conventional system. A repeating block of heat exchangers at three depths was modelled in a series of 3D multi-physics ﬁnite element simulations. The models included the ﬂuid-pipe and pipe-soil heat transfer calculations, the porous nature of the soil, and the seasonal freezing and thawing at the ground surface. Two operation modes were deﬁned, representing the non-balanced mode and the seasonal balancing through the injection of lake water in summer. The climate conditions at the ground surface were explicitly modeled by deﬁning a temperature boundary condition. A novel AI framework was developed to generate the projections of ground surface temperatures from air temperature and other atmospheric variables. An artiﬁcial neural network was trained using the air and ground temperature measurements at the nearby weather station. Several training approaches were compared to minimize prediction errors. The model was then used to convert the air temperatures from downscaled outputs of the CanESM2 climate model into surface temperatures. Six ﬁnite element models, representing the two operation modes under three climate change scenarios, were simulated and veriﬁed by comparing the resulting ground temperatures with the local measurements. The outputs were processed into extraction power, energy, and carbon emissions. Results suggest that in order to attain a stable heat extraction throughout a year, heat exchangers should be placed at depths that are not aﬀected by seasonal variations of surface conditions. Thermal depletion was observed during the ﬁrst years of operation. The seasonal balancing method showed enhancement in total extraction capacity by 55%. The system was found to behave resiliently when subjected to the major climate pathways. Based on the average carbon footprint of electricity production in Canada, the studied system can save up to 800 tons and 1300 tons of equivalent CO 2 emissions over its life, under the non-balanced and balanced operation modes, respectively.


Introduction
Residential buildings and the commercial sector in Canada consume about 27% of secondary energy, which accounts for 23% of the country's greenhouse gasses (GHG) emissions.About 60% of this energy is used for heating and air conditioning (Natural Resources Canada, 2016).The shares are similar in the USA, with 38% of energy being consumed by households and businesses (US Energy Information Administration, 2019).Although the energy efficiency of buildings is continuously being improved, the residential and commercial demand in Canada is expected to increase by 0.3% and 0.6% per year, respectively, due to population growth (National Energy Board, 2018).Of all the G20 members, Canada and the USA are among the top producers of per capita CO 2 emissions.Substituting fossil fuels with cheap, clean and sustainable sources of energy is now a priority in developed countries.As a signing member of the Paris Climate agreement, Canada plans to cut down GHG emissions to 30% below 2005 levels by 2030.The development and promotion of renewable energy-sourced HVAC systems can be a sustainable approach to reduce the country's carbon footprint.
Geothermal heat pumps (GHPs) have recently grown in popularity for heating and cooling of buildings (Omer, 2008).The ground temperature below a depth of a few meters remains almost constant throughout the year.The soil is usually warmer than the ambient weather during winter.This energy can be harvested in several ways.In a conventional closed-loop GHP system, fluid is circulated through ground heat exchangers (GHEs), which transfer the heat between the fluid and adjacent soil.The temperature of the heat carrier fluid increases as it absorbs heat from the ground due to conduction.The fluid is then directed to a heat pump to extract its energy and provide heat to the building at the desired temperature.A similar but reverse approach can be taken for cooling of buildings during summer since the soil is colder than the ambient temperature.
In cold regions, where the warm season is rather short, there is often an imbalance between heating and cooling demands.The excessive heat extracted from the ground through a GHP system may not be replenished naturally and may result in excessive cooling of the ground over the years, known as thermal depletion.While a GHP system may perform well during the first year of operation, depletion can significantly reduce the output of the system in the long run (Saaly and Maghoul, 2019).This has been the subject of a few studies on the sustainability of vertical heat exchangers, such as borehole heat exchangers and thermal piles (Hein et al., 2016;Cai et al., 2019).However, no sustainability study is conducted on horizontal closed-loop heat exchangers, which are a favourable alternative in North America due to the availability of land and lower investment cost.The horizontal and vertical arrangements behave differently.Unlike the vertical exchangers, which are installed under or integrated into building foundations, the horizontal arrangement can be affected by climate conditions at the ground surface (Desideri et al., 2011).In cold regions in Canada's north, a part of the soil domain is seasonally or permanently frozen.The latent heat plays a significant role in the surface energy balance, and hence should be taken into account.Also, the ongoing global warming is affecting atmospheric temperature and precipitation patterns, which all affect surface energy balance (IPCC, 2013).Therefore, a full-scale sustainability study requires the incorporation of atmospheric conditions at the surface, climatic variations throughout the lifespan of the system, and the multi-physics aspect of the problem, such as soil's freeze-thaw cycles near the ground surface.
Efforts have been made to mitigate thermal imbalance due to disproportionate seasonal demand, and also to improve the long-term sustainability of GHP systems (Saaly and Maghoul, 2019).Assisting GHPs with captured solar energy is widely discussed in the literature for vertical heat exchangers (Kjellsson et al., 2010;Wang et al., 2012).This method is also used for the storage of heat (Chapuis and Bernier, 2009).However, using solar collectors may not be cost-effective in high latitudes due to lower solar flux.On the other hand, shallow lakes and ponds absorb significant solar energy.It has been generally shown that water bodies can be incorporated into GHP systems as ground-lake hybrid loops (Chiasson et al., 2000).Freshwater bodies cover about 9% of Canada, including up to 15% of Quebec, Ontario, and Manitoba (Statistics Canada, 2012).The proximity to lakes and ponds makes freshwater bodies a suitable alternative for seasonal balancing.Despite being frozen for several months a year, temperatures in shallow lakes of arctic regions can reach up to 20 • C during the warm season.The warm lake water can be used, directly or indirectly, to replenish the ground energy, or in other words, to store heat for the next cold season (Fathalikhani et al., 2019;Fatolahzadeh et al., 2019).However, the long-term effectiveness of this method is not yet known.
This study aims to investigate both short-term and long-term performance of GHP systems, and the resilience and mitigation potential of the technology considering climate change impacts through a series of multi-physics 3D finite element (FE) simulations of a sample GHP system.The simulations include thermal analyses of heat exchangers and the soil body, taking into account the porous nature of soils, groundwater, and seasonal freeze-thaw cycles.The daily performance of the system is analyzed over a study period of 50 years.In order to incorporate the climate trends and the interaction between atmospheric elements and the ground surface in the model, a novel deep learning approach is proposed to convert projections of air temperatures into ground temperatures.The short-term performance, i.e., the extraction power within a year, is conducted through the analysis of heat exchangers located at different depths.In order to study the long-term sustainability, the annual extracted energy over the study period is evaluated, and the effectiveness of the seasonal balancing using warm lake water is examined as a means against thermal depletion.The resilience of the system is investigated by studying its response to different climate change pathways, and based on the current CO 2 intensity of electricity generation, some remarks on the carbon footprint of the system under different climate pathways are given.

GHP System Specifications
The sample GHP system is based on a monitored system in the city of Winnipeg, Canada.However, the study is not specific to the site since the majority of subarctic Canadian midlands and the northern states in the US prairies have a relatively similar climate.These regions have long and cold winters, with average daily lows that can reach to -20 • C in January.Summers are short and warm, with average daily highs that exceed 25 • C in July (Environment and Climate Change Canada).The system has 30 heat exchangers, each 100 meters long, implemented by horizontal (directional) drilling.The heat exchangers are arranged in parallel and interconnect two main pipes laid down in trenches.The outer and inner diameters of the exchangers are 37.5 mm and 30 mm, respectively.Water circulates through the system, and it is assumed that it does not freeze since the heat exchangers and the main pipes are implemented below the freezing depth.The ground profile consists of a silty clay layer, reaching down to at least 15 meters, and the bedrock underneath.The heat exchangers are located at the depths of 4.5 m, 9.5 m, and 13.5 m below the ground surface, all within the clayey soil layer.Due to the heating demand in the region, the system operates from the beginning of October until the end of May every year.During this period, water is injected into the heat exchangers at 1 • C.

Conceptual Model
The conceptual model consists of a two-layered domain of clay and bedrock, according to the stratification.It extends down to the depth of 30 meters to eliminate the thermal impacts of the lower boundary on the heat exchangers.A temperature boundary condition, T s (t) is applied to the top of the domain to model the equivalent interaction between weather elements and the ground surface.The bottom of the domain is set to a constant upward heat flux to represent the geothermal flux at the site.As the heat exchangers are implemented in a repetitive pattern, the model contains one block, in order to reduce redundancy and save computation effort.The subsurface domain, dimensions of the block, and the arrangement are shown in Figure 1.The heat exchangers, defined as pipes, can thermally interact with the surrounding soil through heat conduction.It is assumed that the flow of heat carrier fluid (water) is distributed equally between heat exchangers.However, to maintain the symmetric condition, i.e., eliminating any thermal interaction between two adjacent blocks, the sides of the domain are defined as adiabatic, and the heat exchangers on the sides carry half the flow circulating in the heat exchangers within the block.Assuming no heat loss at any part of the system other than the heat exchangers, the model includes the GHEs and the adjacent domain only.The thermal performance of a single heat exchanger can be calculated as: where Q [J] is the amount of heat extracted or released , V [m 3 /s] is the heat carrier fluid flow rate, C pf [J/(kg • K)] is the specific heat capacity of the fluid, and T i and T o [ • C] are the inflow and outflow temperatures, respectively.In this regard, the fluid at a specific temperature is injected into one end of pipes, and the outflow temperature at the outlets are to be logged.Since the hydraulic conductivity of the clay and bedrock is significantly low, convective heat transfer is negligible compared to conduction.Therefore, it is assumed that the heat transfer is merely through conduction.

Governing Equations
The governing partial differential equation of conductive heat transfer in a solid medium is: is the heat source or sink , and q [W/m 2 ] is the conductive heat flux from Fourier's law as: where λ [W/(m • K)] is the thermal conductivity and ∇T [K/m] is the temperature gradient.In a saturated porous medium, the governing equation for the transient conductive heat transfer, i.e., assuming no heat transfer due to water movement between pores, is written as (Liu et al., 2019): where ρ, C and λ are the bulk density, specific heat capacity, and thermal conductivity of the saturated freezing medium, respectively.The last term of the equation accounts for phase change, i.e., seasonal freezing and thawing in the ground, in which L f [J/kg] stands for the latent heat of freezing for pore water, ρ i [kg/m 3 ] is the density of pore ice, and θ i is the pore ice content.Again, Q is the heat source/sink, such as heat released from or absorbed by a ground heat exchanger.Eq. 4 can be written as: where the first term is known as apparent heat capacity: The density of the saturated frozen soil (ρ) can be obtained as: in which ρ s , ρ w , ρ i and θ s , θ w , θ i are the density and volumetric fraction of solid grains, pore water, and pore ice, respectively.Similarly, the specific heat capacity of soil is defined as: where C ps , C pw , and C pi represent the specific heat capacity of the solid grains, pore water, and pore ice, respectively.In a similar manner, thermal conductivity of a saturated freezing medium is assumed as: in which λ s , λ w , and λ i are the thermal conductivity of the solid grains, pore water, and pore ice, respectively.In a GHP system, the coupled conductive heat transfer between the heat exchanger and the surrounding medium can be written as (Dacquay et al., 2020): where ρ f , λ f , C pf are the density, thermal conductivity, and specific heat capacity of the heat carrier fluid, respectively.d h [m] is the hydraulic diameter, A [m 2 ] is the inner area of the heat exchanger, and the velocity of the fluid is denoted by u.

Numerical Model Setup
In this study, the numerical models were simulated in COMSOL Multiphysics v5.4.The simulations included heat transfer in solids and pipes, with the implementation of latent heat and apparent properties of the saturated porous medium subjected to freeze-thaw cycles.The latter were defined as dependent variables, which were updated at each computation increment based on the current ice content.The thermal and hydraulic parameters of the clay and bedrock were taken from previous studies in the area (Ferguson and Woodbury, 2004;Saaly et al., 2020), or were assumed from similar materials.These parameters are summarized in Table 1, for the clay and bedrock layers.Previous borehole measurement in the area indicated a geothermal flux of 0.38 W/m 2 (Jessop and Judge, 1971), which was applied to the bottom of the domain as a heat flux boundary condition.Ground surface temperatures, as explained in the next section, are estimated from the projections of air temperature using an artificial neural network (ANN) curve fitting approach, and applied to the top surface of the model.A mesh of tetrahedral elements is defined, as shown in Fig. 2. Based on a preliminary sensitivity analysis, the maximum size of elements in the clay domain was limited to 3 meters, and the soil elements in contact with the heat exchangers were further refined.Larger elements were used at the base of the bedrock domain to improve computational efficiency.

Density
Heat Capacity Thermal Conductivity The study period consisted of a warm-up and the main operation phases.Since the initial ground temperature profile did not cover the entire depth of the model domain, a 5-year warm-up phase was defined before the main study period to allow the ground to reach a state of equilibrium.Verification was done by comparing the ground temperature profile at various instances with those measured at the site.The analysis proceeded through the study period of 50 years-from 2020 to 2070-to simulate the long-term performance of the system.The models were run in two operating scenarios and under three climate change pathways.In the non-balanced scenario, the injection temperature was set to 1 • C. The system operated only during winter, and the flow was set to zero during summer.The second operation mode incorporated seasonal balancing, with the same injection temperature (1 • C) during winter, but the system continued to operate during summer with water at lake water being injected into the heat exchangers.The measurements at Lake Devonian, a shallow lake located in the city of Winnipeg, were chosen for this purpose.The lake has a maximum depth of 8 m and reaches a maximum temperature of 22 • C in summer.The two operation scenarios are referred to as non-balanced and balanced, for reference.With the investigation of climate change impacts as one of the objectives, the incorporation of the surface energy system in the FE models was a major challenge of this study.Although the surface energy balance (SEB) analysis is the most realistic modelling approach to obtain ground temperatures, it has some practical limitations.A comprehensive SEB analysis requires complex boundary conditions reflecting climate elements at the ground surface, e.g., wind speed, solar radiation, vegetation, and snow depth.Also, the simulation becomes computationally intensive due to the added physics, up to the point of becoming impractical over a long study period.In such cases, ground surface temperatures (T s ) are often applied to the surface boundary to implicitly model the effects of surface energy components.Since one of the objectives of this study is to assess the performance of GHP systems over the next decades, the boundary condition includes future temperatures, accordingly.A few global climate models are able to predict the ground surface temperature through a basic analysis of energy balance.However, such predictions are not suitable for engineering applications.Ground surface temperatures depend on surface elements that vary on a local scale, while the spatial resolution of global climate models is often in the order of hundreds of kilometres.On the other hand, climate models have more focus on atmospheric predictions.High-resolution projections of air temperature (T a ) and other atmospheric variables can be obtained from downscaled regional climate models (Scinocca et al., 2016).In many engineering studies, ground surface temperatures are estimated from air temperature.A widelyused approach assumes a sine function for ground surface temperatures, based on the freezing and thawing n-factors that are defined as the ratio of the number of days with negative (or above-zero) temperatures in the soil and the air (Carlson, 1952;Lunardini, 1978).However, as global warming alters the freezing and thawing indices, the n-factors are subjected to change, and the method is not valid for the prediction of future climate trends.
Alternatively, if the air and soil temperatures at the site are available for a period of time, the data can be used to train an ANN to detect patterns, trends, and the correlation between the two variables, which can be incorporated into a transfer function.For the studied case, a dataset of daily records was available, which included air and soil temperatures, as well as other climate variables.The dataset contained 3,163 entries from 2011 to mid-2019 and was used for the training of an ANN model.The model had ten hidden neurons and was trained by Levenberg-Marquardt algorithm using the ANN fitting toolbox in MATLAB.Of about nine years of recorded data, the first six years were used for training/validation/test in a 70:15:15 ratio, and the rest was used to benchmark the extrapolation.Several training scenarios were defined and compared in terms of root mean square prediction errors (Table2).Scenario 1 represented the daily pairs of the ground and air temperature, as dependent and independent variables, respectively.Since the ground surface temperature at an instance can be a delayed response to the past changes in air temperature, Scenarios 2-4 involved a time series of changes in air temperatures, i.e., the past two, seven, and fourteen days, as independent variables.In Scenario 5, an additional climate variable-solar radiation-was included.Scenarios 6-7 involved the daily extremes in air temperature as independent variables.It should be noted that the resulting transfer function is only valid for the site since the model is trained by site-specific data.The resulting prediction errors indicate the effects of different training variables and including time series as the independent variable.While lower errors could be attained by defining new scenarios or through optimization, it was not in the scope of this study.Therefore, the transfer function trained by Scenario 4 was chosen to convert the projections of air temperature from downscaled results of Canadian Earth System Model 2 (CanESM2) into ground surface temperatures.The daily air temperatures from CanESM2 may differ with the temperatures measured at the station due to the temporal accuracy of the climate model (Fig. 3).Nevertheless, both show similar variations on a monthly scale.Surface temperatures, on the other hand, exhibit smaller daily fluctuation.The plots also demonstrate conformity between the ground surface temperatures predicted from the two sources, i.e., the measurements and and CanESM2.To take different cli-mate pathways into account, ground surface temperatures were generated under RCP2.6,RCP4.5 and RCP8.5 scenarios, and were applied to the finite element models, accordingly.Plotting the mean annual ground surface temperatures (MAGST) reveals the warming of the ground surface over the 21st century, though the trends diverge after 2070 (Fig. 4).

Validation
In total, six models were simulated, which account for the response of the system under the non-balanced and balanced operation modes and three climate pathways.No records of outlet temperatures were available at the site to validate the analysis of heat transfer in pipes.Nevertheless, heat transfer in the soil domain was validated through a comparison of measured and calculated ground temperature profiles in the last year of the warm-up cycle.The measurements of ground temperatures at the site were only available up to the depth of 1 m.Therefore, temperatures at the maximum depth were used for validation (Fig. 5).

Data Processing Framework
For each simulation, fluid temperatures at the inlet and outlet nodes of the top, middle, and bottom heat exchangers were recorded.These outputs were used to calculate the extracted energy, cost, or extraction power.According to Eq.1, the amount of heat extracted by an exchanger is proportionate to the density and heat capacity of the heat carrier fluid, pumping rate, and also the difference between the inflow and outflow temperatures.Assuming water as the heat carrier fluid and an equal distribution of constant fluid flow between pipes, the extraction power of exchangers were calculated.The power was averaged and multiplied by the duration of the operation season, to calculate the annual extracted energy.

Short-term Performance
Analyzing the extraction power from heat exchangers located at the three depth levels reveals a difference in thermal behavior within a year (Fig. 6).Focusing on the early years of operation, while the highest extraction at the beginning of the cold season each year was observed from the top level, it exhibited a drop from 300 W to less than 200 W per heat exchanger.On the other hand, the deeper heat exchangers showed less than 20% decline and produced steadier power throughout the operation season.This is due to the soil around the heat exchangers at the top level being close to the ground surface, and therefore being heavily affected by the seasonal fluctuations at the surface.The annual variation of ground temperature profile (Fig. 7) at the level of top exchangers (4.5 m) is about 3 • C, while it is less than 0.5 • C at 9 m.According to Fourier's law of heat conduction, higher thermal gradients result in higher heat fluxes within the soil medium surrounding the heat exchangers.Another explanation is that at the beginning of the cold season, when the shallow ground is still warm from the previous summer due to the ground thermal mass, the extraction power by the top heat exchangers is at its peak.It then rapidly declines because of the heat being extracted by the exchangers, and also the cooling of shallow soil due to atmospheric conditions in winter.However, the effect of surface conditions is more significant since the three depth levels have the same injection temperature.In a similar manner, the top level replenishes more quickly in summer due to its proximity to the surface, and in the case of the balanced operation, the contribution by the injection of warm lake water.
The peak heating demand in the Canadian Prairies occurs in January and February every year, while the calculated extraction power from the top layer is maximum at the beginning of the operation season, i.e., October.Although the observation of annual extracted energy from the heat exchangers (Fig. 8a) indicates higher extraction at the top level, the aforementioned short-term instability should be considered in the design of the system.The heat exchangers at middle and bottom levels, on the other hand, have more stable heat extraction.However, it should be noted that the deep implementation requires longer drilling and a higher chance of reaching bedrock, which may result in a higher implementation cost.

Thermal Depletion and Long-term Sustainability
The extraction of heat through a GHP system alters the ground's natural energy budget.Thermal depletion happens when the rate of heat extraction is higher than the rate of replenishment by geothermal flux or surface absorption so that the extraction declines each year.Reviewing the performance of single exchangers (Fig. 8a) reveals a decline in annual extraction during the first five years of operation despite an increase in MAGST due to climate warming (Fig. 8b).This initial drop, which indicates thermal depletion, was not noticeable for the heat exchangers at the top level, but it was about 8% at the middle and bottom levels over the course of five years.The results, however, showed a strong correlation between the annual extraction from the top exchangers and MAGST.This was also visible at the middle and bottom levels and indicates that the extracted energy from the heat exchangers is influenced by the annual weather patterns and the long-term impacts of climate change (Fig. 4).Thermal depletion cannot be simply evaluated by analyzing the annual thermal yield at the beginning and the end of the study period.Eliminating the climatic variations through simulating over a constant ground surface temperature for the duration will solve the issue.However, it would lead to erroneous conclusions as some years since it ignores the impacts of climate change.Therefore in this section, the results were analyzed under the three climate pathways to reveal the extent of thermal depletion while considering the possible replenishment due to global warming.Probing the annual thermal extraction of the block operating under non-balanced and balanced modes (Fig. 9, the depletion pattern was visible during the first years of operation under all climate pathways when the system operated in the non-balanced mode though it was slightly diluted as the block consists of heat exchangers at the three levels.However, the depletion will be counterbalanced after two decades by intensive warming of the ground surface due to climate change.The balanced operation mode, on the other hand, shows a significant improvement in the energy extracted per block.It starts with a rapid boost in the first decade and then gradually increases at a slower rate over the rest of the study period.Under the RCP 2.6 pathway, the balanced approach results in 55.7% higher total output than the non-balanced mode over the design life.About 70% of this improvement occurs within ten years from starting the balancing procedure (10).The method was also found effective in other climate scenarios than RCP 2.6, showing improvement by 55.3% and 54.0% under RCP 4.5 and RCP 8.5, respectively.

Climate Resiliency
With the looming impacts of climate change on built infrastructure and communities, it is important to highlight the influence of climate change on the performance of GHP systems, i.e., resiliency.Increasing atmospheric temperature and solar radiation alter the ground temperature, which can affect the performance of GHP systems over time.Analyzing the annual extraction under RCP 2.6, RCP 4.5, and RCP 8.5 scenarios (Fig. 9), though showing an increase, did not indicate a significant variation in extracted heat with respect to the climate pathways.In the case of non-balanced operation, the slight increase due to climate change barely surpassed the depletion.However, it should be noted that the climate pathways diverge in the second half of the century.The maximum difference between MAGSTs predicted by the optimistic RCP 2.6 and the worst-case RCP 8.5 during the study period is less than 1.5 • C, but it can reach up to 4 • C at the end of the century (Fig. 5).Moreover, the extracted energy was presented for a block of heat exchangers, in which only the shallow ones are significantly affected by surface temperatures.Although the increased thermal outputs due to global warming may be viewed as an opportunity, the maximum contribution due to global warming is 11% during the study period, which is about five times smaller, compared to seasonal balancing.

Carbon footprint and climate change mitigation
GHP systems, though a source of thermal energy, require electricity to operate.The Coefficient of Performance (COP), i.e., the ratio of extraction to consumption, is between 3 to 5 for conventional systems.It means that a GHP system can save electricity consumed in the heating of buildings, ranging between 66% to 80% of its nominal extraction capacity.Yet the equivalent reduction in GHG emissions depends on the carbon footprint of electricity generation in the region.According to Canada's Official Greenhouse Gas Inventory, the CO 2 emission intensity of electricity generation is 130 g CO 2 /kWh.However, it drastically varies between provinces due to different means of electricity generation.For instance, in Quebec, Manitoba and BC, where hydroelectricity is abundant, the intensity is much lower than the national average, less than 5 g CO 2 /kWh.On the other hand, in Alberta, Saskatchewan, and Nunavut, it reaches up to 790 g CO 2 /kWh, due to the use of fossil fuels (Canada Energy Regulator, 2017).Assuming the national average carbon intensity of electricity generation, the equivalent CO 2 savings were calculated for the studied system Fig. 11.Over a period of 50 years, the system can save more than 800 t CO 2 eq, when operating in the non-balanced mode.In case the system operates in the balanced mode, the reduction can reach up to 1300 t CO 2 eq.Since a significant portion of generated electricity in North America is consumed for residential and commercial HVAC systems, geothermal technology can contribute to the mitigation of climate change, along with other sources of renewable energy.

Conclusion
Based on the results, in cold regions, where the heating demand is significant, GHP systems can provide a stable yield throughout the cold season if the exchangers are implemented in a depth that is not affected by seasonal variations of surface weather.Thermal depletion is observed in all depth levels, mostly in the deepest exchangers during the early years of operation, and should be considered in performance studies.Seasonal balancing can relieve depletion in the short term, and also significantly enhance the extraction capacity in the long term.It was also shown that the GHP technology could remarkably reduce carbon emissions, and therefore, can be promoted by public policymakers to reduce the future carbon footprint, especially in remote communities that currently use fossil fuels for electricity generation.However, while the simulations can predict the system behaviour up to a predefined tolerance, the principal model assumptions and climate projections are associated with uncertainties.The most optimistic and worst-case climate scenarios were simulated to show the extent of the system responses.The results show slightly improved thermal outputs due to global warming, which may be considered a benefit for a system merely used for heating.The results also indicate the system to be resilient to climate variations.However, the rise of air temperatures shifts the energy consumption in buildings toward higher cooling demands, which can be critical to the performance of dual-purpose GHP systems.Furthermore, while the GHP technology can significantly reduce the carbon footprint of residential and commercial buildings, the calculations are based on the current emission intensity of electricity generation.Recent trends from the last few decades show emissions per kWh to be on the decline due to environmental pressures and technological advances.Therefore, the estimation of future carbon footprint may significantly vary between cases.

Figure 1 :
Figure 1: Model domain and the arrangement of heat exchangers within the block

Figure 2 :
Figure 2: Finite element mesh of tetrahedral elements

Figure 11 :
Figure 11: Carbon footprint reduction over 50 years, based on Canadian average carbon intensity of electricity generation

Table 1 :
Material properties used in finite element simulations

Table 2 :
Summary and the prediction errors of ANN training scenarios T aavg : Average daily air temperature T amax , T a min : Maximum and minimum daily air temperature RS avg t i : Daily average solar radiation