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 2emissions over its life, under the non-balanced and balanced operation modes, respectively. Keywords: geothermal energy, sustainability, climate change, adaptation, mitigation 1. 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 eciency 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 Preprint submitted to Journal As of: August 13, 2020 Board, 2018). Of all the G20 members, Canada and the USA are among the top producers of per capita CO 2emissions. 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, uid is circulated through ground heat exchangers (GHEs), which transfer the heat between the uid and adjacent soil. The temperature of the heat carrier uid increases as it absorbs heat from the ground due to conduction. The uid 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 rst year of operation, depletion can signi cantly 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 di erently. Unlike the vertical exchangers, which are installed under or integrated into building foundations, the horizontal arrangement can be a ected 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 signi cant role in the surface energy balance, and hence should be taken into account. Also, the ongoing global warming is a ecting atmospheric temperature and precipitation patterns, which all a ect 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. E orts 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-e ective in high latitudes due to lower solar ux. On the other hand, shallow lakes and ponds absorb signi cant 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 2 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 e ectiveness 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 nite 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 di erent depths. In order to study the long-term sustainability, the annual extracted energy over the study period is evaluated, and the e ectiveness 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 di erent climate change pathways, and based on the current CO 2intensity of electricity generation, some remarks on the carbon footprint of the system under di erent climate pathways are given. 2. Methods 2.1. GHP System Speci cations The sample GHP system is based on a monitored system in the city of Winnipeg, Canada. However, the study is not speci c 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 pro le 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. 2.2. Conceptual Model The conceptual model consists of a two-layered domain of clay and bedrock, according to the strati cation. 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, Ts(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 ux to represent 3 the geothermal ux 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 e ort. The subsurface domain, dimensions of the block, and the arrangement are shown in Figure 1. The heat exchangers, de ned as pipes, can thermally interact with the surrounding soil through heat conduction. It is assumed that the ow of heat carrier uid (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 de ned as adiabatic, and the heat exchangers on the sides carry half the ow 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: @Q @t=_VC pf(TiTo) (1) whereQ[J] is the amount of heat extracted or released , _V[m3=s] is the heat carrier uid ow rate,Cpf[J=(kgK)] is the speci c heat capacity of the uid, and TiandTo[C] are the in ow and out ow temperatures, respectively. In this regard, the uid at a speci c temperature is injected into one end of pipes, and the out ow temperature at the outlets are to be logged. Since the hydraulic conductivity of the clay and bedrock is signi cantly low, convective heat transfer is negligible compared to conduction. Therefore, it is assumed that the heat transfer is merely through conduction. Figure 1: Model domain and the arrangement of heat exchangers within the block 2.2.1. Governing Equations The governing partial di erential equation of conductive heat transfer in a solid medium is: Cp@T @t+rq =Q (2) 4 in which[kg=m3] is the density of the medium , Cp[J=(kgK)] is the speci c heat capacity , T [K] is the soil temperature, Q[W=m3] is the heat source or sink , and q [W =m2] is the conductive heat ux from Fourier's law as: q=rT (3) where[W=(mK)] is the thermal conductivity and rT[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): Cp@T @t=r(rT) +Q+Lfi@i @t(4) where,Candare the bulk density, speci c 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 Lf[J=kg] stands for the latent heat of freezing for pore water, i[kg=m3] is the density of pore ice, and iis the pore ice content. Again, Qis the heat source/sink, such as heat released from or absorbed by a ground heat exchanger. Eq. 4 can be written as: CpLfi@i @T@T @tr(rTs) =Q (5) where the rst term is known as apparent heat capacity: Capp=1  CpLfi@i @T (6) The density of the saturated frozen soil ( ) can be obtained as: =ss+ww+ii (7) in whichs,w,iands,w,iare the density and volumetric fraction of solid grains, pore water, and pore ice, respectively. Similarly, the speci c heat capacity of soil is de ned as: Cp=1 (ssCps+wwCpw+iiCpi) (8) whereCps,Cpw, andCpirepresent the speci c 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: =1 (sss+www+iii) (9) in whichs,w, andiare 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): AC app@T @t+fAC pfurT =rAfrT +A 2dhjuj3QA (10) wheref,f,Cpfare the density, thermal conductivity, and speci c heat capacity of the heat carrier uid, respectively. dh[m] is the hydraulic diameter, A[m2] is the inner area of the heat exchanger, and the velocity of the uid is denoted by u. 5 2.2.2. Numerical Model Setup In this study, the numerical models were simulated in COMSOL Multiphysics v5.4. The sim- ulations 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 de ned 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 ux of 0.38 W=m2(Jessop and Judge, 1971), which was applied to the bottom of the domain as a heat ux boundary condition. Ground surface temperatures, as explained in the next section, are estimated from the projections of air temperature using an arti cial neural network (ANN) curve tting ap- proach, and applied to the top surface of the model. A mesh of tetrahedral elements is de ned, 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 re ned. Larger elements were used at the base of the bedrock domain to improve computational eciency. Table 1: Material properties used in nite element simulations Density Heat Capacity Thermal Conductivity [kg=m3]Cp[J=(kgK)][W=(mK)] Soil Skeleton 2560a942a1.58a Water 1000 4180 0.56 Ice 917 2000 2.24 Bedrock 2600 1200 2.30b a(Saaly et al., 2020) b(Ferguson and Woodbury, 2004) The study period consisted of a warm-up and the main operation phases. Since the initial ground temperature pro le did not cover the entire depth of the model domain, a 5-year warm-up phase was de ned before the main study period to allow the ground to reach a state of equilibrium. Veri cation was done by comparing the ground temperature pro le 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 ow 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. 6 Figure 2: Finite element mesh of tetrahedral elements 2.2.3. Prediction of the Surface Temperature Boundary using Arti cial Neural Networks (ANN) 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 re ecting 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 ( Ts) are often applied to the surface boundary to implicitly model the e ects 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 ( Ta) 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 widely- used approach assumes a sine function for ground surface temperatures, based on the freezing and thawing n-factors that are de ned 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 7 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 tting toolbox in MATLAB. Of about nine years of recorded data, the rst 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 de ned 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 sur- face 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 radi- ation| 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-speci c data. Table 2: Summary and the prediction errors of ANN training scenarios Scenario Training Set Entry RMSE [C] Description 1Taavg ti2.8932 Current day, average daily air temp. 2Taavg ti;Taavg ti;ti1;Taavg ti;ti22.4449 Last 2 days, average daily air temp. 3Taavg ti;Taavg ti;ti1;:::;Taavg t1;ti72.1615 Last 7 days, average daily air temp. 4Taavg ti;Taavg ti;ti1;:::;Taavg ti;ti142.1044 Last 14 days, average daily air temp. 5Taavg ti;RSavg ti2.7227 Current day, average daily air temp. and radiation 6Tamax ti;Tamin ti2.9083 Current day, daily max/min air temp. 7Tamax ti;Tamin ti;Tamax ti;ti1;Tamin ti;ti12.5837 Last 2 days, average daily air temp. Taavg: Average daily air temperature Tamax;Tamin: Maximum and minimum daily air temperature RSavg ti: Daily average solar radiation The resulting prediction errors indicate the e ects of di erent training variables and including time series as the independent variable. While lower errors could be attained by de ning new scenarios or through optimization, it was not in the scope of this study. Therefore, the trans- fer 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 tempera- tures. The daily air temperatures from CanESM2 may di er 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 uctuation. The plots also demonstrate conformity between the ground surface temperatures predicted from the two sources, i.e., the measurements and and CanESM2. To take di erent cli- 8 mate pathways into account, ground surface temperatures were generated under RCP2.6, RCP4.5 and RCP8.5 scenarios, and were applied to the nite 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). Figure 3: Daily TaandTsconversions from CanESM2 RCP4.5 vs measurements at St. Adolphe Station, Manitoba, Canada 2.3. 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 tem- peratures 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 pro les 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). 2.4. Data Processing Framework For each simulation, uid 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 uid, pumping rate, and also the di erence between the in ow and out ow temperatures. Assuming water as the heat carrier uid and an equal distribution of constant uid ow 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. 9 Figure 4: MAGST pathways calculated from daily Taprojections of CanESM2 RCP 2.6, RCP 4.5, and RCP 8.5 3. Results and Discussion 3.1. Short-term Performance Analyzing the extraction power from heat exchangers located at the three depth levels reveals a di erence 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 a ected by the seasonal uctuations at the surface. The annual variation of ground temperature pro le (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 uxes 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 e ect of surface conditions is more signi cant 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 10 Figure 5: Veri cation of the ground temperatures (June 2018 - July 2019) that the deep implementation requires longer drilling and a higher chance of reaching bedrock, which may result in a higher implementation cost. 3.2. 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 ux 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 rst ve 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 ve 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 in uenced 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 rst 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. 11 Figure 6: Extraction power per heat exchanger (2020-2030) [RCP 4.5] The balanced operation mode, on the other hand, shows a signi cant improvement in the energy extracted per block. It starts with a rapid boost in the rst 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 e ective 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. 3.3. Climate Resiliency With the looming impacts of climate change on built infrastructure and communities, it is important to highlight the in uence of climate change on the performance of GHP systems, i.e., resiliency. Increasing atmospheric temperature and solar radiation alter the ground temperature, which can a ect 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 signi cant 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 di erence 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 signi cantly a ected 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 ve times smaller, compared to seasonal balancing. 12 Figure 7: Calculated monthly ground temperature pro les (July 2018 - June 2019) 3.4. Carbon footprint and climate change mitigation GHP systems, though a source of thermal energy, require electricity to operate. The Coe- cient 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 Ocial Greenhouse Gas Inventory, the CO 2emission intensity of electricity generation is 130 g CO 2/kWh. However, it drastically varies between provinces due to di erent 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 in- tensity of electricity generation, the equivalent CO 2savings were calculated for the studied system Fig.11. Over a period of 50 years, the system can save more than 800 t CO 2eq, 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 2eq. Since a signi cant 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 signi cant, GHP systems can provide a stable yield throughout the cold season if the exchangers are implemented in a depth that is not a ected 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 signi cantly enhance the extraction capacity in the long term. 13 Figure 8: a) Extracted energy per heat exchanger at the three depth levels [RCP 4.5] b) MAGST [RCP 4.5] 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 prede ned 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 bene t 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 signi cantly 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 signi cantly vary between cases. Acknowledgments This research was supported by the New Frontiers in Research Fund - Exploration Grant [NFRF-2018-00966], and the University of Manitoba, Graduate Enhancement of Tri-Council Stipends (GETS) program. The authors would like to acknowledge KGS Group and Mr. Timi Ojo (Province of Manitoba) for providing the eld data used in this research and Mr. Rob Sinclair (KGS Group) for his technical assistance.