Subsurface geothermal in combination with Ground Heat Exchanger (GHE) and heat pumps systems had proved over the past decade to be an effective way in space heating and cooling. GHE’s embedded in geostructures such as tunnels has become more and more attractive in the current state where we try to mitigate greenhouse gas effects. The scope of this paper is to provide a fully coupled numerical model to assess the different physical aspects of an energy tunnel such as the kinetic and thermal boundaries layers development, thermally induce stress, and the heat generation. To do so, a finite element model is used to assess the heat exchange between, the GHE, the concrete lining, the air within the tunnel, the surrounding ground, and groundwater flow. The numerical model is then confronted with data from real-scale experimentation that was performed in Torino in 2018. The predictions from the numerical model show a quite good match with the data from the experimental prototype. Later on, similar models are generated to assess the effect of the heat exchange on air velocity and temperature within the tunnel as well as the induced thermal stress. The numerical results highlight the formation of both the kinematic and thermal boundary layer and shows that those boundaries are affected by the GHE’s activation. For the kinematic aspect, the lower temperatures observed in the vicinity of the heat exchangers tend to lower the velocity magnitude especially close to the wall surface. This change in velocity may be linked to a buoyancy effect cause by the lower temperature of the air in this regions. Indeed natural convection occurs along the tunnel, cooler fluid is pushed down while the warmer fluid rise; locally a lower temperature tends to decrease the viscosity of the fluid and increase its volumic mass. It was observed that the activation of GHE’s enhances this phenomenon in their vicinity due to the temperature gradient that they produce. Moreover, GHE’s activation tends to increase the overall turbulence of the flow and thus the eddies’ formation. The locations of those disturbance are difficult to predict due to the fact that small disturbance may generate large vorticity within the flow. However, it is shown that after passing the GHE zone, the velocity profile tends to stabilize to the undisturbed case. The kinematic entry length seems to happen further down the stream in the case where GHE’s are activated. This is due to the temperature disturbance generated by GHE’s activation. This phenomenon is increased with respect to the number of GHE loops activated. Furthermore, some probing was performed in the vicinity of the heat exchanger to assess the thermal induce stress due to the heating operation. The recorded stress shows peaks at early stages due to the higher temperature gradient that is generated between the GHE’s and concrete at beginning. Finally, after some time, the temperature reaches a steady state hence the stabilized thermal response recorded.