Geothermal pavements: field observations, numerical modelling and long-term performance

One of the key challenges of the twenty-first century is moving towards clean and renewable energy resources to both meet the ever-growing energy demand and reduce mankind's carbon dioxide footprint. Heating and cooling accounts for nearly 50% of global energy consumption, and most of the energy resources comes from burning fossil fuels (REN21, 2019). Ground source heat pump (GSHP) systems are considered to be one of the most efficient and environmentally benign solutions to be used in space heating and cooling (Brandl, 2006; Johnston et al., 2011; Kim et al., 2016). The typical coefficient of performance (COP) value of a GSHP system is around four to five, which indicates that every 1 kW electricity input can produce 4 to 5 kW heating/cooling energy (Narsilio et al., 2014). Generally, GSHP systems utilise the ground as a heat source in winter and a heat sink in summer by way of ground heat exchangers (GHEs). Based on the GHE configurations, GSHP systems can be categorised into horizontal (İnalli & Esen, 2004) and vertical (Hepbasli et al., 2003) systems. Restrictions to the wider utilisation of GSHP systems include the typically high capital costs, mainly due to the costly installation of GHEs, which may include the cost of drilling boreholes for vertical GSHP systems (Qi et al., 2019), or the need for a large installation area for horizontal GSHP systems.

Recent advances in energy geostructures have been increasingly gaining attention, owing to their strong potential for cost savings related to GSHP systems, as well as the adoption of GHEs in urban environments to solve the land availability issue. This is where the expertise of geotechnical engineers on heat transfer in soils (and rocks) can make a difference to the technology. Energy geostructures incorporate high-density polyethylene (HDPE) piping to underground structures that are primarily designed for stability, as well as to exchange heat. Hence, foundation piles (Brandl, 2006; Bourne-Webb et al., 2016; Makasis et al., 2018a), retaining walls (Barla et al., 2020; Makasis et al., 2020; Shafagh et al., 2020) and tunnel linings (Adam & Markiewicz, 2009; Bidarmaghz & Narsilio, 2018) can serve as GHEs with little additional cost, and achieve a dual purpose of structural stability and thermal provision. Although energy geostructures can reduce the installation cost significantly, the design and implementation of GHEs highly depends on the site layout and geometry of the underground structures. In addition, the heat transfer between GHEs can cause additional stress on structural stability, which needs to be considered carefully during design (Katzenbach et al., 2014; Batini et al., 2015). These limitations can raise complexity in designing geothermal systems that can provide the entirety of the desired thermal loads. Therefore, for some geothermal systems, it is more beneficial or necessary that they are used in combination with an auxiliary energy system to fulfil the thermal demand, thus forming a hybrid system.

The majority of existing energy geostructure research focuses on energy piles, given their similar geometry to traditional vertical GHEs such as boreholes. There are extensive studies focusing on the utilisation of analytical and numerical solutions to explore the thermal performance of energy piles (Abdelaziz et al., 2011; Loveridge & Powrie, 2013, Makasis et al., 2018b), developing further understanding on the response of piles to thermo-mechanical loads (Amatya et al., 2012; Batini et al., 2015; Gawecka et al., 2017) and proposing design options for improving the system performance (Alberdi-Pagola et al., 2020). A summary of recent analyses on energy geostructures, particularly on energy piles, can be found in Loveridge et al. (2020). It was found that research related to energy geostructures involving horizontal GHEs, such as geothermal pavements, has been largely overlooked.

The aim of the current research is to investigate the design of geothermal pavement systems, which unlike other energy geostructures, have not received much attention. With the rapid increase in human population leading to urbanisation and increased transportation, the idea of geothermal pavements becomes a promising solution and one that involves geotechnics professionals and pavement geotechnics expertise. Geothermal pavements take advantage of existing earthworks for road construction by placing the horizontal GHEs under the pavement surface at a relatively shallow depth (less than 1 m depth), thereby reducing capital costs for the GHE installation. The majority of existing studies focus on using geothermal pavement for road maintenance, such as preventing snow accumulation and controlling pavement surface temperatures (Eugster, 2007; Gao et al., 2010; Muñoz-Criollo et al., 2016; Zhao et al., 2020). As such, the option of considering geothermal pavement systems for thermal control of buildings remains largely unexplored in the available literature. In addition, a major difference lies in the much shallower placement of GHEs in pavements than in other typical geothermal horizontal systems. The impact of this feature has not been sufficiently studied to date.

The present work focuses on the thermal performance of geothermal pavements. Since GHEs play a crucial role in the heat exchange mechanism of a GSHP system, it is imperative to model accurately the behaviour of a shallow horizontal GHE. Although the thermal performance of traditional horizontal GHEs (buried at depths of 1–3 m) has been investigated by means of both experimental studies and analytical solutions (Gan, 2013; Jeon et al., 2018; Lamarche, 2019), there exists a paucity of research pertaining to very shallow horizontal GHEs (less than 1 m) and their modelling. To gain deeper insights into the heat transfer process, researchers are moving towards numerical solutions, owing to their wider applicability and fewer constraints in exploring geometry and configuration than laboratory or field-scale experiments permit (Selamat et al., 2016; Han et al., 2017; Kayaci & Demir, 2018; Yang et al., 2020). The soil thermal properties have the greatest impact on subsurface heat transfer; however, most of the developed numerical models assume constant soil properties throughout the modelling process for the purpose of simplification and saving computational time (Kayaci & Demir, 2018; Yang et al., 2020). In a recent study, Han et al. (2017) created a three-dimensional (3D) finite-element (FE) model to assess the thermal performance of various horizontal GSHP system configurations, as well as exploring the importance of using local geological data (measured ground temperatures and soil properties) as design inputs for improving modelling accuracy. In addition, regarding the proximity of the installation of GHEs to the ground surface, Li et al. (2019) studied the effects of different ground surface boundary conditions on the thermal performance of horizontal GHEs subjected to high building loads when buried at a shallow depth.

The literature review suggests that there are some numerical models available for shallow geothermal energy systems with research on traditional horizontal GHEs. However, models related to GHEs placed at a depth of less than 1 m, as is the case of geothermal pavement systems, are yet to be explored extensively. Although researchers have investigated geothermal pavements from a heat transfer perspective, their results were based on pilot or laboratory-based experimental data dealing with short-term performance (del-Tota-Maharaj et al., 2011; Castillo-García et al., 2013). To bridge these gaps, the aim of the current research work is to create a 3D FE numerical model of geothermal pavements based on non-uniform soil properties that utilises local geological data, with the goal of studying a geothermal pavement system as a means to harnessing sustainable energy for space heating and cooling. This model is successfully validated with full-scale experimental data collected at a site in Adelaide, South Australia (SA). In addition, the long-term performance of the system is simulated and evaluated using the validated model in conjunction with the thermal load of a typical residential building in Adelaide, SA. Apart from thermal performance study of GSHP systems, it is also important to assess the economic feasibility of designed systems. This is often required by engineers and clients and a pertinent question that remains largely unanswered in the literature; the authors intend to redress this situation herein to complement the technical study and for the benefit of the readership. The work by Nicholson et al. (2014) represents a notable exception, but with a focus on the planned geothermal tunnel of London's Crossrail project. In the work presented by Barla et al. (2016) on energy tunnels, apart from their main study on exploring the thermal efficiency and sustainability, they also carried out a cost analysis and comparison with traditional heating and cooling systems to highlight the economic benefit of energy tunnels. Franzius & Pralle (2011) investigated the thermal capacity and sustainability of implementing a heat-exchanger segment lining system for tunnel-boring machine-driven energy tunnels. Although Franzius & Pralle (2011) did not carry out an economic analysis in their paper, the economic viability and additional value for investment was emphasised. Cousin et al. (2019), through a more detailed analysis of the economic attractiveness (including levelised cost of electricity, net present value (NPV) and return on investment (ROI)) of energy segmental linings, concluded that every energy geostructure has to be carefully designed according to specific conditions, as a poor design can jeopardise both the energy efficiency and economic attractiveness of the application. Aditya et al. (2020) performed a comprehensive comparative lifetime cost analysis for GSHP systems against conventional heating and cooling systems for different climatic conditions. They concluded that the economic performance of GSHP systems depend on various parameters, and there is no ‘one size fits all’ rule. Life-cycle cost (LCC) analysis is a useful economic indicator for examining the application of GSHP systems (Zhu et al., 2012; Cui et al., 2018; Lu & Narsilio, 2019). Therefore, an LCC analysis is performed in this research to explore the financial benefits of such a system while drawing comparisons to other widely available conventional systems. This complements the study of heat transfer mechanisms and thermal yields of geothermal pavements.

This research develops a validated numerical model for a geothermal pavement system with the main aim of using this model to investigate the system's thermal performance. Detailed 3D FE numerical modelling approaches were utilised, incorporating the governing equations to simulate the thermal performance of the system under annual thermal loads of a typical residential dwelling in Adelaide, SA. Furthermore, the thermal performance of the system under more balanced loads was assessed by introducing a hybrid system approach, which incorporates the geothermal pavement system with an auxiliary heating/cooling system such as a reverse cycle air conditioner (RCAC) to provide the required thermal loads. Lastly, a LCC analysis was carried out to compare both capital and operating costs of the geothermal pavement system against a vertical GSHP system, a horizontal GSHP system, a RCAC system and a hybrid system.

The modelling of the geothermal pavement system was implemented in the FE package, Comsol Multiphysics, focusing on the GHEs embedded in the pavement. This 3D numerical FE model has been developed within the University of Melbourne, and validated herein with experimental data from a full-scale field study in Adelaide, SA (Motamedi et al., 2020, 2021) (further validations included in the Appendix). The governing equations for fluid flow and heat transfer are numerically coupled to capture the heat transfer between the carrier fluid in the pipes and surrounding ground.

The momentum and continuity equations (1) and (2) were used for modelling the incompressible fluid flow inside the HDPE pipes. Both heat conduction and convection are considered for modelling the heat transfer between carrier fluid and the ground. Conductive heat transfer mainly occurs in the solid materials, including the ground, GHE pavement structure backfills and the HDPE pipe walls, while convective heat transfer dominates in the carrier fluid within the HDPE pipes. In the absence of groundwater flow, the conductive heat transfer equation of the soil, HDPE pipe walls and GHEs was modelled by equation (3). In addition, the conductive–convection energy equations for heat transfer in the fluid flow are given in equations (4) and (5) (Lurie, 2008).

where $A$ is the inner cross-section of the HDPE pipe (m²); $ρw$ is the carrier fluid density (kg/m³); $v$ is the fluid velocity field (m/s) within the pipes; $t$ is time (s); $p$ is pressure (P_a); $fD$ represents the Darcy friction factor; $dh$ is the hydraulic diameter of the pipe (m); $ρm$ is the solid material density (kg/m³); $Cp,m$ is the specific heat capacity of a solid material (J/kg K); $Tm$ stands for the temperature of a solid material; $λm$ is the thermal conductivity of a solid material (W/m K); $Cp,w$ is the specific heat capacity of the fluid (J/kg K); $λw$ is the thermal conductivity of the fluid (W/m K); $Qwall$ is the external heat exchange rate through the pipe wall (W/m), and a function of the temperature of the pipe outer wall, $Tm,pipewall$ (K), and the temperature of the carrier fluid T (K).

The FE model, geometry and configuration can be seen in Fig. 1 and the relevant material properties are summarised in Table 1. The geothermal pavement system consists of a 0·05 m surface layer (asphalt), 0·30 m thick base layer with a mix of gravel (0·25 m) and fine sand (0·2 m), and the remaining layer is a highly plastic clay, which is the subgrade layer of the pavement. A single GHE circuit consisting of eight legs of HDPE pipes, and having an outer diameter of 25 mm (SDR11), is placed at a depth of 0·5 m below the surface and connected in series meandering in an area of approximately 4 m × 20 m. The total length of the HDPE pipes is 160·4 m with a spacing of 600 mm between the centres of the pipes. Moreover, the constant fluid flow rate through the inlet of the circuit is 12 l/min. Since this model is based on the conditions of a site in Adelaide, and given the shallow depth, groundwater flow is not considered in this research. The material characteristics were obtained following Australian Standard AS 1289 (Standards Australia, 2017), while the thermal conductivity of each layer was estimated by employing a successive forward simulation technique (Santamarina & Fratta, 2005). Key input parameters were successively varied in the numerical model and the resulting outlet fluid temperature, as well as ground temperatures compared against field measurements. The effective thermal conductivities of the pavement layers were determined as the thermal parameters that render the lowest error between simulated and measured data. Values obtained are consistent with those expected for the different materials and with thermal needle probe measurements (Manohar et al., 2000).

Appropriate assumptions, initial and boundary conditions are crucial to obtaining meaningful results. The geothermal pavement system including boundary conditions is shown in Fig. 1(a) and is summarised below.

• (a)

  The fluid flow within the pipes is assumed to be fully developed; therefore, the fluid flow and heat transfer occurring in the fluid are analysed as a one-dimensional problem.

• (b)

  The initial temperature of both the ground and the far-field boundaries are equal to the measured undisturbed ground temperature T for different depths z – that is, $Tfarfieldz$⁠, which can be seen in Table 2.

• (c)

  Thermal insulation that renders a zero net flux is prescribed to the far-field side boundaries of the numerical model. The distance between the edge of the GHEs and the outer boundary is determined as 10 m, making sure that a larger distance produces close to identical results, to ensure there is sufficient distance to avoid any boundary influences (similar results are achieved while assuming fixed temperatures over the simulation period at the boundaries). The temperature of the model bottom is set as a prescribed temperature on the boundary.

• (d)

  The surface boundary of the model is designed as a time-dependent ambient temperature $Tairt$⁠. The ambient air temperature data are the recorded air temperatures in Adelaide, SA.

• (e)

  The time-dependent inlet fluid temperature $Tint$ is determined through a function of the carrier fluid temperature at the outlet, $Toutt$ (obtained from solving the numerical model) and from the prescribed time-dependent annual thermal loads $QGHEt$ of a typical residential building in Adelaide. This simulates the operation of a heat pump, where $Toutt$ is the fluid temperature entering the heat pump, and $Tint$ is the temperature of the fluid leaving the heat pump after extracting/rejecting the required thermal demand.

  $Tint=Toutt−QGHEtm˙Cp,w$

  6

  where $Q$ is the provided thermal loads (W); $m˙$ is the mass flow rate (kg/s); and $Cp,w$ is the specific heat capacity of the fluid (J/(kgK)).

General GSHP systems can be defined as efficient and functional if the operating fluid temperature in the pipes is within a reasonable temperature range (typically between −5 and 40°C) for the life of the project (typically 20–25 years) (Makasis et al., 2018c). These temperatures are based on a number of parameters, most importantly, the thermal load distribution, which defines the thermal energy that the GSHP system needs to extract from/reject to the ground, the system geometry and the site conditions, among others. Since the GHEs of the geothermal pavement system are located at a shallower depth than traditional horizontal GSHP systems, these temperatures (and thus the overall system performance) can also be more sensitive to the ambient ground temperature. In this study, annual thermal loads of a typical residential building in Adelaide, SA, as shown in Fig. 2, are used with the FE method numerical model to determine the carrier fluid temperatures. In this study, the average fluid temperature $Tfluid$ (average from $Tin$ and $Tout$⁠) is used to quantify the circulated fluid temperature within the pipe.

This section introduces the results of the model validation and presents analysis results of evaluating the thermal performance of the geothermal pavement system. The model validation is discussed in the next section entitled ‘Model validation’, and then the thermal performance of different model configurations, including a stand-alone geothermal pavement system design and a hybrid system, are presented in the section entitled ‘Investigation on the thermal performance of the geothermal pavement systems’. Furthermore, a sensitivity analysis exploring the effect of changing the surface boundary condition on system performance is provided in the section entitled ‘Sensitivity analysis on surface boundary condition’. To determine the feasibility and functionality of a geothermal pavement system, it is also important to evaluate its economic implications; therefore, the section entitled ‘Cost comparison with various heating and cooling system configurations’ presents a preliminary cost analysis on capital and operating costs. The geothermal pavement system is compared with conventional space heating and cooling systems, including a traditional vertical GSHP system, a horizontal GSHP system and a RCAC system.

The numerical model must first be validated before being utilised to simulate the system's annual thermal performance. To gain confidence in the model created, it should be compared against experimental data for validation. Also, to reinforce the assumption of incorporating $Tairt$ for the surface boundary condition, a numerical model which assigns a general energy balance equation to the ground surface is discussed and compared. In addition, more detailed validation is carried out by modelling the ground soil temperature at three locations where thermistor sensors have been installed. In this numerical model, the mesh uses free tetrahedral elements, and the final mesh size used herein is determined to balance the computational time and accuracy, with finer mesh refinement within the region near the structural elements. The FE mesh of the model consists of 525 453 domain elements. The mesh element size was carefully chosen following validation, ensuring that a finer mesh provides close to identical results.

A geothermal pavement system field study was conducted in 2019 at a carpark in Adelaide, SA. Temperate climate conditions prevail in this part of the country. To explore the thermal benefits of geothermal pavement systems, this field study carried out a thermal response test (TRT) for 1·5 days with a constant heat power of 4·5 kW and a flow rate of 16·8 l/min. A single GHE circuit is approximately 160 m long (between inlet and outlet, plus header piping). The geometry and instrumentation of the field set-up are illustrated in Fig. 3 and material details were further discussed earlier in the section entitled ‘Model geometries and model parameters’. In addition to the measured inlet and outlet temperatures from the TRT unit, thermistors were installed along the length of the pipes to the soil to help gain understanding of the heat transfer. Further details of this unique experimental set-up were reported by Motamedi et al. (2020, 2021).

Selection of the surface boundary condition is arguably one of the most debatable assumptions while modelling for horizontal GHEs. Two surface boundary condition choices are explored and discussed in this section. Boundary condition 1 (BC1) assumes ground surface temperature as a function of time (equal to $Tair$⁠), which has been detailed in the section ‘Model assumptions, initial and boundary conditions’, whereas boundary condition 2 (BC2) introduces a surface boundary condition where the heat flux due to varying weather conditions is modelled (Chiasson, 2010)

where n is the normal vector on the boundary; $λm$ is the material thermal conductivity (⁠$W/mK$⁠); $T$ represents the material temperature (K); $h$ is the heat transfer coefficient due to wind; $Tair$ is the ambient air temperature (K); $q′′solar$ stands for the solar irradiance on the surface (⁠$W/m2$⁠); $ε$ represents the ground surface thermal emissivity, which is dimensionless; $σ$ is the Stefan–Boltzmann constant (⁠$5⋅67×10−8W/m2K4$⁠); and $Tsky$ is the sky temperature (K). These ambient thermal properties were obtained from the processing of measured data from the American Society of Heating, Refrigeration and Air-Conditioning Engineers (ASHRAE) weather data viewer (ASHRAE, 2017) and from the Australian Bureau of Meteorology (particularly solar radiation) for Adelaide, SA (Bureau of Meteorology, 2019). Note that parameters including $Tair$⁠, $q′′solar$ solar radiation on the horizontal, wind speed, $Tsky$ and $T$ are time dependent.

The ground surface boundary condition in BC2 considers the solar radiation entering the pavement surface along with the net radiation from the surface to the surroundings and convection heat transfer caused by the wind. The comparisons between modelled and measured outlet fluid temperature for both boundary conditions are shown in Fig. 4. On the one hand, both numerical models show good agreement with measured outlet fluid temperature. On the other hand, the numerical model with the proposed BC1 shows a slightly better overall simulation performance with a lower root mean square error (RMSE) value of 0·36°C than the model with BC2. While these observations are valid for the TRT duration, they indicate that applying a time-dependent $Tair$ to the ground surface boundary is a simple and appropriate assumption and is capable of accurately simulating the thermal behaviour of the geothermal pavement system. It should be noted that wind speed fluctuations during testing were relatively minor. The selection of boundary conditions is worth investigating in detail and will be contrasted against longer term full-scale experimental data in future research.

The (arguably small) discrepancies in the results could be attributed to the fact that the sensors used to record the inlet and outlet fluid temperatures were attached to the outside of the pipes rather in direct contact with the water. Although well thermally insulated from the ambient, these attached sensors may be more susceptible to the ambient air temperatures and solar radiation than if placed inside the pipes. For example, fluctuations in the measured $Tout$ that can be seen around the 15th hour of the test occurred in the middle of the day when the air temperature was higher compared to the rest of the day. The effects of these are seen on the upper right part of Fig. 4(b). Moreover, inconsistencies observed at the initial stage of the test (i.e. lower temperatures, lower left part of Fig. 4(b)) could be caused by modelling initialisation errors or discrepancies in the measured initial ground temperature data (Table 2), which could influence the initial modelling results.

For a more complete validation, the temperature distribution in the pavement system is computed with the model described in the earlier section entitled ‘Methodology’ with T_air(t) as surface boundary, and compared against field measurements. Fig. 5 shows a comparison between the modelled and measured soil temperatures at the location of thermistors 1, 2 and 4 (see Fig. 3). The comparison for sensor 3 is not shown in Fig. 5, because the ground temperatures observed in T3 and T4 are very similar. The modelled temperatures at both T1 and T2 have captured the measured results very well. The sensor at location T4 (and similarly T3), however, shows a less good agreement. Even though the largest variation between measured and modelled at location T4 is less than 1°C, the overall trend of the modelled soil temperature for the T4 sensor is consistently lower than the measured data with lower fluctuations. Since the GHE of the geothermal pavement system is located at a shallow depth of 0·5 m, the soil temperatures can be strongly influenced by the surface temperature, and unlike locations T1 and T2 the effect of the GHE operation is less dominant at this depth. Note that the RMSE is in the order of 0·5°C, the same as the accuracy of the thermistors used in the experiment.

In summary, the 3D numerical model is capable of simulating the performance of the geothermal pavement system while also showing a very acceptable agreement with measured data, not only in terms of fluid temperature, but also capturing the soil temperature distribution well.

In this section, the validated numerical model is employed to evaluate the long-term thermal performance of the geothermal pavement system, with a case study centred on a typical residential building in the temperate climate conditions of Adelaide (SA) (see earlier section ‘Thermal performance evaluation’). In a traditional vertical GSHP system, unbalanced thermal loads lead to thermal accumulation effects on the ground (Bidarmaghz et al., 2016; Li et al., 2018). This is generally not the case for horizontal systems. To ensure there is no thermal drift or accumulation in the pavement package, a 5 year thermal load case (yearly repetition of the thermal loads in Fig. 2) is simulated. Based on the results shown in Fig. 6, it is evident that there are no thermal accumulation effects on the horizontal GHE performance, since the $Tfluid$ shows a yearly cyclic performance. Note that, given the annual net balanced (ground and fluid) temperature, it is unlikely that the starting time of operation of the (horizontal) systems (e.g. in the middle of winter or in the middle of summer) would have a long-term impact on thermal performance. Based on the results from Fig. 6, an annual (1 year) simulation of a typical residential dwelling in Adelaide, SA is used in this research to represent the performance throughout the system's design life. Two geothermal system configurations are analysed: a standalone geothermal pavement system, the details of which are discussed in the next section ‘Thermal performance of the geothermal pavement system’; and a hybrid system (using the geothermal pavement to provide the base thermal loads (6 kW) and an RCAC system to support shortfalls). Results are analysed in the section entitled ‘Thermal performance of the hybrid geothermal pavement system’.

As mentioned in the section entitled ‘Thermal performance evaluation’, the geothermal pavement system design is considered acceptable if the average fluid temperatures $Tfluid$ over the simulation are within a reasonable operating range. An average $ΔT=Tin−Tout$ of about 1·5°C is observed from the modelling results. Fig. 7 presents the simulated $Tfluid$ of multiple GHE circuits under annual thermal demands. According to Fig. 7, the $Tfluid$ of two identical GHE circuits reach temperatures outside the typical acceptable GSHP operating range. Therefore, in this case, it was found that four identical GHE circuits can meet the thermal loads resulting in a highest $Tfluid$ value of 39·8°C and lowest of 0°C while maintaining a suitable $ΔT$ between inlet and outlet for optimal system operation.

To obtain a better understanding of the thermal performance of the system, Fig. 8 shows the temperature distribution corresponding to the hottest and coldest day during the simulation period, 1 February (building cooling required) and 16 July (building heating required), at the cross-section A–A indicated in Fig. 1(b). The higher differential temperature and rate of heat transfer occurs in both heating and cooling modes on the right-most pipes. This is a subtlety captured by the detailed model and is expected since the inlet is placed on the right side of the GHE in this section and, depending on whether the system is operating in heating or cooling mode, the entering water extracts or rejects the maximum amount of heat due to the temperature differential between the fluid in the pipes and the surrounding ground. As the carrier fluid travels through the pipes, the rate of heat transfer decreases moving left. It is also noticed, in Fig. 8(a), that the temperature of ground surface is similar to that of the vicinity of the pipes (they fall approximately within the same temperature contour interval). This can be explained by considering the leakage of heat from the region surrounding the pipe towards the surface. Furthermore, this also indicates that the shallower pipe depths of geothermal pavements may degrade the heat extraction capacity of the surrounding ground during extremely hot (or cold) days.

In the common scenario of a large residential development with multiple dwellings and a finite pavement area available, it may not be possible to assign four GHE circuits per dwelling. In this case, a hybrid system may be considered and result in a more economical solution and better operation of the system. Here one must find the maximum amount of thermal energy that the geothermal pavement can provide (per dwelling). With this aim, an iterative process varying the peak heating and cooling capacity of a GSHP such that the resulting average fluid temperature is found to be within operational ranges is undertaken. As a result of these multiple successive forward simulations, the hybrid system considered in this study uses the geothermal pavement system to supply a ‘base load’ of thermal demands that are lower than 6 kW (for both heating and cooling) or approximately 40% of the peak heat demand and 33% of the peak cooling demand, while the exceeding thermal loads are supplied by an auxiliary RCAC system. The resulting $Tfluid$ values for the geothermal pavement system comprising a hybrid design are displayed in Fig. 9, noting that an anti-freeze solution rather than plain water may be required – in this case, 15% of propylene glycol is recommended by International Ground Source Heat Pump Association (IGSHPA) guidelines (Remund & Carda, 2009). Given the ‘peaky’ thermal load pattern found in Adelaide, despite the modest base load (from 33 to 40%), the hybrid system can provide 78% of the annual thermal load using the geothermal pavement system, requiring the use of the auxiliary system for 137 days per year.

The location of geothermal pavement systems designed for both standalone and hybrid system configurations are presented in Fig. 10. Compared to the four-GHE circuit design as discussed in the section ‘Thermal performance of the geothermal pavement system’, the hybrid system requires only two GHE circuits, which would significantly reduce the capital costs as well as the required construction area. Therefore, instead of serving a single building, the same area could be utilised to supply the demand of two separate buildings with the help of an auxiliary system. Note that with this concept design, not all residential buildings in a new neighbourhood would be able to be conditioned with a GSHP system (unless additional GHEs are built on front or back yards). The cost comparison for these systems can be found in the section ‘Sensitivity analysis on surface boundary condition’.

Owing to the very shallow buried depth of the GHEs, it is important to explore the influence of the ambient air temperature on the thermal performance of the designed pavement system. A sensitivity analysis is carried out in an attempt to understand the impact of varying air temperature on the $Tfluid$ of the proposed geothermal pavement designs. All parameters remain unchanged except $Tairt$⁠, which is modified by simply adding values in a range of −4 to 4°C with an increment of two. The resulting relationships between maximum and minimum $Tfluid$ to $Tairt$ for both system design configurations are presented in Fig. 11. The maximum and minimum $Tfluid$ (⁠$Tmax$ and $Tmin$⁠) shows a linear direct proportionality to changes in $Tairt$⁠, with correlation coefficient R² values for all linear regressions above 0·99. Moreover, the slopes for both standalone and hybrid geothermal designs are greater than 0·9, which indicates that a change of 1°C of $Tairt$ would result in more than 0·9°C change in the peak $Tfluid$ values expected.

To understand the cost efficiency and not only the technical but also the financial feasibility of geothermal pavement systems, a comparison of life-cycle analysis is undertaken for a traditional vertical GSHP system, a horizontal GSHP system and an RCAC system. All the systems are assumed to be operated for a design life of 20 years. In addition, based on Lu et al. (2017) and Aditya & Narsilio (2020), the lifespans of a GSHP and a RCAC unit are assumed to be 20 and 10 years, respectively. Thus, every tenth year, a new RCAC unit must be installed to replace the previous one, which incurs additional capital costs. The net present cost (NPC) method is used in this study as the indicator for all the heating and cooling systems considered. The calculation of NPC can be seen in equation (8)

where $Ct$ is the net cash outflow during a single period $t$ in years, which can be either capital cost or an operating cost. The discount rate $r$ is taken as 3·5% (see Table 3) for the present cost at time $t$ over a 20-year design life of the system.

It is essential to investigate the economic implication of the geothermal pavement system, especially with regard to the capital costs, which can often be the deciding factor on whether a system design will become reality. The operating cost is also important to consider, mainly by examining the system COP. For example, different design configurations can result in different $Tfluid$ ranges, which could influence heat pump efficiency and thus the operating costs. Besides, both capital and operating costs are the key components for calculating the cash outflow in equation (8). The assumptions adopted for this cost analysis can be found in Table 3 and are based on the literature, the current economic status of Australia and the research group's experience with system installations. The Adelaide thermal load distribution is considered as the design thermal load for determining the pipe length of traditional GSHP (vertical and horizontal trench) systems. The designed borehole length of a vertical GSHP system is 330 m (based on a separate numerical modelling not presented in this work), while using the model developed in this study, for a GHE placed at 2·0 m depth horizontal GSHP system, the design pipe length is 600 m. The cost of a 20 kW heat pump is assumed to be Australian (AUD) $13 300 for the stand-alone geothermal system designs, and the cost of a 6 kW GSHP (used in hybrid system design) is AUD $6000 (note: 1 AUD = US$0·78). The average drilling cost for the borehole is considered to be AUD $130 per metre (Lu et al., 2017). The installation and material cost of pipes are provided in the assumptions as well. Since the operating temperatures of fluid within pipes may fall below zero, the cost of adding antifreeze solution is also considered in this research. Regarding the annual running cost, the COP values of each system are determined as a basis for the comparison. The daily COP values of the geothermal pavement system and vertical and horizontal GSHP systems were computed based on the water temperatures entering the GSHP, along with the technical specifications provided by Climate Master (Climate Master TC Series 50 Hz TCH/V072 – HFC410A). In addition, it should be noted that while calculating the COP of the geothermal pavements within the hybrid system, a smaller heat pump was used (Climate Master^TM TC Series 50 Hz TCH/V030 – HFC410A), and the COP of RCAC was obtained from Qi et al. (2019). The operating cost of each system was computed based on the system's daily COP values, primarily a function of the T_out(t). The annual averaged COP values are listed in Table 4, and the analysed NPC over the 20 year lifespan for each system is presented in Fig. 12.

The total NPC of the designed geothermal pavement system resulted in AUD $69 268, which is found to be the most cost-effective among all the systems considered. From the observation in Fig. 12, geothermal pavement system has the lowest capital cost, as it takes advantage of the construction of pavement, and hence does not require any additional earthworks, such as drilling boreholes or trenching and excavation. Although the operating cost of geothermal pavement system is higher compared to traditional vertical and horizontal GSHP systems, the much lower capital cost gives pavement systems great economic potential. Unlike the geothermal systems, a RCAC system relies on ambient air as the exchange medium rather than the fluid in the GHEs, which has more stable temperatures, thus providing more efficient space heating and cooling. Therefore, RCAC systems have the lowest annual average COP values among all the systems, amounting to the highest running cost. In addition, the design life of an RCAC is only 10 years, which results in an additional replacement cost in the tenth year of the design life period. Thus, the RCAC system is the most expensive system and requires 52% extra cost compared to the geothermal pavement system. Although the vertical GSHP system appears to be the most economical in terms of operational cost, the requirement for expensive drilling results in high capital cost, which is a major drawback for that option. The overall LCC for the vertical GSHP system is just slightly lower than for the RCAC system, while exceeding the pavement system by 50%. As for the horizontal GSHP system, the greater GHE depth (2 m) saves 6% in design pipe length but only reduces the running cost by 1%. However, as mentioned above, there are additional trenching costs associated with it. The hybrid system (geothermal pavement with two GHEs and auxiliary) also provides a viable alternative. As shown in Fig. 12, the capital cost of the hybrid system is slightly more than that of the standalone geothermal pavement system due to the additional installation and replacement cost of the RCAC system. The advantage of considering the hybrid system is that this system only requires two GHE circuits, which would be beneficial if the available road space is limited.

Apart from the NPV cost analysis, it is also interesting to look at the ROI for these structures, herein considering a design period of 20 years. ROI is calculated by dividing the profit gained on an investment by the cost of that investment. Since the RCAC system results in the highest NPV and is also a common system used in residential houses, the RCAC system would be assumed as the initial thermal control system in the house. Therefore, the profits gained by the other proposed systems in Fig. 12 are the total cost savings compared with the RCAC system. The computed ROI values for each system follow the same trend as Fig. 12, where geothermal pavement has the highest ROI value of 52%, which is 29% higher than the second-best performing system, the hybrid system (23%). As regards the vertical GSHP system, its ROI value is only 1%, and the ROI of the horizontal GSHP system is 14%. It is worth noting here that, if a longer design period is taken into account, these values are expected to increase. Among the designed systems, the geothermal pavement design yields the highest ROI value, by a significant margin, which suggests its higher economical attractiveness.

Energy geostructures have recently gained increasing attention for providing space heating and cooling, alongside contributing towards continued efforts in relation to climate change. Among the most popular energy geostructures are energy piles, energy tunnels and energy retaining walls, but limited information exists relating to geothermal pavements. This research has utilised numerical modelling and available field data to explore the potential and the long-term thermal performance of geothermal pavement systems. A FE numerical model has been developed and subsequently validated with experimental data. Annual thermal demand of a typical residential building in Adelaide, SA was utilised to investigate the long-term performance of the system. This study has focused on the performance of the following two system design configurations: one uses a standalone geothermal pavement system to supply the required thermal loads, while the other system consists of a combination of a geothermal pavement system to provide a base load and an auxiliary RCAC system to supply the balance of the thermal demand. In addition to the thermal performance, life-cycle analysis for different suitable heating and cooling systems has been presented, to better understand the potential and cost implication of the geothermal pavement system. The key findings are summarised below.

• (a)

  A developed FE modelling methodology has been presented and validated against two different cases, one found in the literature and one from a field TRT test undertaken in Adelaide, SA.

• (b)

  Under design thermal loads of a typical dwelling in Adelaide, four GHE circuits (having total HDPE pipe length of 640 m) of the geothermal pavement system can fulfil the thermal demands.

• (c)

  Adopting a hybrid system can reduce the required GHE circuits from four to two, thus reducing the capital costs related to material and installation as well as future maintenance costs of the system. A base load of just about 35% of the peak demand is to be provided by the geothermal system, but providing thermal comfort for 78% of the time, while the balance of the peak load is met by the auxiliary system.

• (d)

  Geothermal pavement systems produce higher cost-saving than traditional geothermal systems. Geothermal pavement systems have the lowest LCC, because this system design does not require any additional excavation or drilling, unlike traditional GSHP systems. The annual running cost of the geothermal pavement system is not the lowest, as vertical GSHP systems have the highest COP throughout the year. However, the vertical GSHP system amounts for the most expensive capital costs and leads to a high NPC value. In addition, a RCAC system is the most costly design because it requires replacement of the RCAC unit, as well as having high running costs compared to the geothermal systems in this study.

• (e)

  Both the hybrid and geothermal pavement system configurations show a great potential in terms of economic implications and thermal performance. Although the hybrid system only needs two instead of four GHE circuits, the overall cost of a hybrid system design configuration is slightly higher than that of a standalone geothermal pavement system. Therefore, if there is a limitation on pavement length or land availability, the hybrid system would be a better design option.

Overall, in this paper a comprehensive study has been conducted, exploring the thermal performance and feasibility of geothermal pavement systems, and results indicate that it is worthwhile to investigate these systems in detail. Further studies could explore the influence of a ground temperature variation with seasonality, since the HDPE pipes are located at 0·5 m depth, and therefore are placed shallower than traditional GSHP systems; in addition, an even more thorough life-cycle analysis could be conducted, including environmental aspects and carbon dioxide emission calculations.

Funding from the Australian Research Council (ARC) (project number LP170100072) and the University of Melbourne is much appreciated. Moreover, the fifth and sixth authors would also like to acknowledge the support from the National Science and Technology Development Agency (NSTDA), Thailand, under Chair Professor program (P-19-52303). Assistance from the City of Mitcham (SA), Dr Mahdi Disfani and Mr Ramin Raeesi in the fieldwork is much appreciated.

Apart from the validation with experimental field observations, this research also presents a validation against published results which can corroborate the goodness of the employed modelling approach. The 3D numerical model was validated with experimental data from a full-scale horizontal GSHP system, gathered at the Yildiz Technical University, Istanbul, Turkey. Details about the experimental setup and the data can be found in (Kayaci & Demir, 2018). In the referred study, the GHEs were buried in the soil at a depth of 1·8 m and consisted of three pipes connected in parallel, each having a length of 40 m and spaced 33 m apart. While the GHE circuits designed by (Kayaci & Demir, 2018) were configured in parallel. Fig. 13 presents the top view of the model geometry and model validation results. Since the experimental inlet fluid temperature to the GHE is an input to the model, the comparison is based on the outlet fluid temperatures, which are the output of the model. Other initial and boundary conditions used in this model includes total flow rate $q$  = 21·3 l/min, far-field temperature $Tfarfield=13°C$⁠, and T_air = 9·3°C kept as an average value throughout the duration of the test. As shown in Fig. 13, the modelled outlet fluid temperature shows good agreement with the experimental data. All the peaks are simulated accurately with the maximum difference between experimental and modelled data being less than 1°C. Overall, the numerical model shows a very good match with the experimental data from the literature and it provides more confidence to further validate the model with the geothermal pavement field observations.

$A$

inner cross-section area of a pipe (m²)

$Cp,m$

specific heat capacity of soil material (J/(kg K))

$Cp,w$

specific heat capacity of carrier fluid (J/(kg K))

$Ct$

net cash outflow during a single period t in years (AUD)

$dh$

hydraulic diameter of pipe (m)

$fd$

Darcy friction factor

$h$

heat transfer coefficient due to wind

$m˙$

mass flow rate (kg/s)

$n$

normal vector on the boundary

$p$

pressure (Pa)

$QGHE$

external heat exchange rate through pipe walls (W/m)

$Qwall$

external heat exchange rate through pipe walls (W/m)

$q′′solar$

solar irradiance on the surface (W/m²)

$r$

discount rate

$T$

material temperature (K)

$Tairt$

time-dependent ambient temperature (°C)

$Tfarfield$

undisturbed ground temperature (°C)

$Tfluid$

average fluid temperature (°C)

$Tinl$

inlet pipe temperature (°C)

$Tm,pipewall$

temperature of solid material at the pipe outer wall (K)

$Tmax$

maximum average fluid temperature (°C)

$Tmin$

minimum average fluid temperature (°C)

$Tout$

outlet pipe temperature (°C)

$Tsky$

sky temperature (K)

$t$

time (s)

$λm$

thermal conductivity of solid material (W/m K)

$λw$

thermal conductivity of carrier fluid (W/m K)

$v$

fluid velocity field (m/s)

$ρm$

density of solid material (kg/m³)

$ρw$

density of carrier fluid (kg/m³)