energies

Article
Modeling Artificial Ground Freezing for Construction
of Two Tunnels of a Metro Station in Napoli (Italy)
Alessandro Mauro 1 , Gennaro Normino 1, *, Filippo Cavuoto 2 , Pasquale Marotta 3
and Nicola Massarotti 1
1 Dipartimento di Ingegneria, Università degli Studi di Napoli “Parthenope”, Centro Direzionale, Isola C4,
80143 Napoli, Italy; [email protected] (A.M.); [email protected] (N.M.)
2 Studio Cavuoto, Via Benedetto Brin, 63/D2, 80142 Napoli, Italy; [email protected]
3 Consorzio di Ricerca per l’Ambiente i Veicoli l’Energia e i Biocombustibili (CRAVEB), Centro Direzionale, Is.
C4–80143 Napoli, Italy; [email protected]
* Correspondence: [email protected]; Tel.: +39-081-5476709




Received: 9 January 2020; Accepted: 3 March 2020; Published: 10 March 2020


Abstract: An artificial ground freezing (AGF) technique in the horizontal direction has been employed
in Naples (Italy), in order to ensure the stability and waterproofing of soil during the excavation of
two tunnels in a real underground station. The artificial freezing technique consists of letting a coolant
fluid, with a temperature lower than the surrounding ground, circulate inside probes positioned
along the perimeter of the gallery. In this paper, the authors propose an efficient numerical model to
analyze heat transfer during the whole excavation process for which this AGF technique was used.
The model takes into account the water phase change process, and has been employed to analyze
phenomena occurring in three cross sections of the galleries. The aim of the work is to analyze the
thermal behavior of the ground during the freezing phases, to optimize the freezing process, and to
evaluate the thickness of frozen wall obtained. The steps to realize the entire excavation of the tunnels,
and the evolution of the frozen wall during the working phases, have been considered. In particular,
the present model has allowed us to calculate the thickness of the frozen wall equal to 2.1 m after
fourteen days of nitrogen feeding.

Keywords: numerical modeling; heat transfer; artificial ground freezing; underground station; metro
in Napoli; GEO heating

1. Introduction
Artificial Ground Freezing (AGF) is a consolidation technique adopted in geotechnical engineering,
when underground excavations must be executed in granular soils, or below the groundwater level [1,2].
The realization of relevant underground structures in urban areas often involves the management of
constructive problems related to avoiding the presence of water in the excavation, especially if the soil
has poor geo-mechanical proprieties. The artificial ground freezing technique has been extensively
used in the last decades as an effective and powerful construction method, which provides ground
support, groundwater control, and structural underpinning during construction. However, the use
of this technology requires a good knowledge of frozen soil behavior and a robust numerical model
able to predict ground movements around the excavation. This is important, especially in densely
urbanized areas, where frost action is detrimental for surrounding structures.
The AGF method consists of letting a refrigerant circulate inside probes located along the perimeter
of the excavation, at a temperature significantly lower than the surrounding soil. The water in the
soil goes from liquid to solid phase, and it forms a block of frozen ground in the area surrounding
the probes.

Energies 2020, 13, 1272; doi:10.3390/en13051272 www.mdpi.com/journal/energies

Energies 2020, 13, 1272 2 of 24

The process is, in general, divided into two different phases: a first “freezing phase”, that ends
when the soil achieves the design temperature needed to start the excavation, and a “maintenance
phase”, which is characterized by heat absorption in order to keep the temperature constant during the
excavation [3].
The main advantages of this technique, among the available ground consolidation and
waterproofing technologies, are: (i) security and compatibility with the environment, since there is no
injection, and the dispersion of products in the ground. Water already present in the ground is, in fact,
frozen, using refrigerant fluids that are never directly in contact with the ground and groundwater,
avoiding contamination phenomena; (ii) applicability to any type of soil, from coarse to fine grain and
rocks [3,4].
Depending on the working fluid used, two types of methods can be identified: (i) the direct
method, which is based on the use of liquid nitrogen entering the probes at a temperature of −196 ◦ C
and released in the atmosphere in gaseous phase at a temperature between −80 ◦ C and −170 ◦ C; (ii) the
indirect method, which is based on the use of a mixture of water and calcium chloride (known as
brine), whose circulation temperature can vary between −25 ◦ C and −40 ◦ C. A combination of the two
previous methods is known as a mixed-method, which uses the direct method for the freezing phase,
and the indirect method for the maintenance phase.
Several numerical and experimental works analyzing the AGF technique are available in the
literature. Colombo [5] in the first part of his work invokes a well-known approximate approach for
the a priori evaluation of the parameters influencing the technique, such as the time required to reach
the target temperatures, or the heat flow rate needed by the plant. The results proposed by the author,
applied to Neapolitan tuff, were compared with those obtained from a series of numerical analyses
conducted by using the finite element method as a discretization technique, and with experimental data
measured on site during freezing operations carried out for the realization of the galleries for the stations
of Piazza Dante and Piazza Garibaldi of metro Line 1 in Napoli (Italy). Papakonstantinou et al. [6] first
analyzed the experimental data of monitored temperatures in the ground during the freezing process,
and then performed a numerical analysis through the FREEZE calculation code, a thermohydraulic
software developed at the ETH in Zurich. The authors found that the thermal conductivity of the soil
is an important parameter to be taken into account, and can be reasonably estimated by a posterior
numerical analysis, if it not known a priori. Subsequently, Pimental et al. [7] analyzed the results
of three applications of the AGF technique in urban underground construction projects, comparing
the experimental data with the thermohydraulic coupled code model FREEZE. The first case study
concerned the construction of a tunnel for the underground in Fürth (Germany) in soft ground with
significant infiltration flow. The second case study concerned a platform tunnel in a metro station in
Naples, and aimed at the determination of relevant thermal parameters through retrospective analysis,
and to compare the results obtained by using the forecasting model with on-site measures. In the
third case, regarding a tunnel under the river Limmat in Zurich, numerical simulations were used to
identify potential problems caused by geometrical irregularities in the well layout, in combination with
infiltration flow. Russo et al. [1] analyzed the experimental data collected during the execution of the
excavation with the AGF technique, and developed a numerical model to evaluate stress in the ground
during the freezing and defrosting of a frozen wall. The focus of the work was the settlement caused by
the tunnel excavation, and the use of the AGF technique to allow the safe digging of a service gallery
located half in the silty sand layer, and a half in the yellow tuff layer, below groundwater. The phases
of the tunnel construction were accompanied by the monitoring of the measurements and control
activities of the effects of the gallery excavation. The measurements collected during the construction
process, in fact, allowed us to monitor the freezing-thawing process, and the change in volume related
to the excavation, providing useful information for the future implementation of similar projects.
Finally, the analysis in the test procedure was conducted using a complete three-dimensional (3D)
model implemented in the DFM Flac3D package.
Energies 2020, 13, 1272 3 of 24

Vitel et al. [8] developed a numerical model considering both the freezing tube and the surrounding
ground. The model is based upon the following principles: (i) heat conduction around the well is
solved by considering vertical heat transfer processes negligible compared to the trans-horizontal heat
transfer; (ii) heat transfer in the freezing probe is reduced to a one-dimensional (1D) calculation. In this
study, convection in the ground was not taken into account with respect to the heat conduction, and
therefore, the effects of groundwater flow were not considered. Vitel et al. [9] developed a numerical
thermohydraulic model in order to simulate artificial ground freezing by considering a saturated and
nondeformable porous medium under groundwater flow conditions. Marwan et al. [10] presented a
thermohydraulic finite element model integrated into an optimization algorithm, using the Ant Colony
Optimization (ACO). This technique allowed researchers to optimize the positions of the freezing
probes with respect to the groundwater flow. Kang et al. [11] combined a freezing method and a
New Tubular Roof (NTR) simulated by thermomechanical coupling analysis. The temperature range
obtained in the freezing process indicated that the thickness of the frozen wall grows of about 2.0 m
after 50 days of freezing. Moreover, the stability of the surrounding ground and the support structures
in the bench cutting phase were also studied. Panteleev et al. [12] focused on the development of a
monitoring system for the artificial ground freezing process for a vertical shaft. The temperature in
the wells was measured by using the fiber optic system Silixa, based on the Raman effect. Alzoubi
et al. [13] have evaluated the development of a frozen wall between two freezing probes, with and
without the presence of groundwater infiltration for 2D geometry, by was using ANSYS. Fan et al. [14]
show a case study concerning the monitoring of frozen wall formation during soil freezing using brine,
then developed a three-dimensional numerical model to analyze the temperature distribution. The
numerical simulation was conducted by using ADINA software.
Based on the analysis of the available literature, the interest of the research community on the
AGF technique is evident. Both experimental and numerical works can be found, however, more
research effort is needed to numerically analyze the evolution of the process for real cases, considering
the geometry of civil works and the development of the freezing probes in the ground. For these
reasons, the authors have developed an efficient transient numerical model to effectively analyze heat
transfer in the soil, and at the same time, save computing resources. The proposed approach is based
on the coupling of a heat transfer model between the freezing probes and the surrounding ground
with a heat transfer and phase change model of the soil, and for the first time in the literature, all of
the phases of AGF process have been reproduced. The model was used as a preliminary predictive
analysis for the construction of two tunnels in Napoli. The model has been validated against the data of
Colombo [5]. After validation, the numerical model has been employed to analyze a real case study of
two tunnels between Line 1 and Line 6 of the metro station in Piazza Municipio, Napoli, southern Italy.
The purpose of the work is to study in detail the heat transfer process during ground freezing for the
realization of the tunnels. The analysis is carried out by employing a FEM-based model, using Comsol
Multiphysics commercial software, to model artificial ground freezing during the whole processe.
The model developed in the present work allows us to simulate, for the first time in the literature,
a mixed-method used for the freezing process, from the first phase based on nitrogen feeding,
a maintenance phase, and a third phase that involves the use of brine. The maintenance phase is
necessary to avoid the freezing of brine in the probes. The novelty introduced by this work relies on
the development of a thermal analysis of the entire artificial ground freezing process, considering all
the phases and the influence of the process on the thermodynamic behavior of a second nearby tunnel
that was also subject to the AGF. Moreover, the excavation phase has been reproduced, by imposing
a convective heat transfer condition related to the presence of men and machines, while the second
tunnel was subject to AGF with a mixed method.
This paper is structured as follows: in Section 2, the characteristics of the AGF technique are
described, while in Section 3 the numerical model developed is presented. Section 4 reports the
validation carried out against literature data. The results of the parametric analysis performed after
model validation are reported in Section 5, while conclusions are drawn in Section 6 of the paper.
Energies 2020, 13, x FOR PEER REVIEW  4 of 26 
Energies 2020, 13, 1272 4 of 24

performed after model validation are reported in Section Five, while conclusions are drawn in the 
last section of the paper.   
2. Description of AGF Technique and Case Study
Artificial Ground Freezing (AGF) consists of freezing the ground by means of heat transfer, with a
2. Description of AGF technique and case study
refrigerant fluid circulating inside probes located along the perimeter of the excavation to be realized.
Artificial Ground Freezing (AGF) consists of freezing the ground by means of heat transfer, with 
In this way, the water contained in the soil undergoes a phase transition from liquid to solid, forming a
a  refrigerant  fluid  circulating  inside  probes  located  along  the  perimeter  of  the  excavation  to  be 
block of frozen ground called a “frozen wall” in the area surrounding the probes.
realized. In this way, the water contained in the soil undergoes a phase transition from liquid to solid, 
The mixed-method AGF process used for the Piazza Municipio galleries can be divided into
forming a block of frozen ground called a “frozen wall” in the area surrounding the probes. 
different phases:
The  mixed‐method  AGF  process  used  for  the  Piazza  Municipio  galleries  can  be  divided  into 
Phase 1-Nitrogen: The probes are fed with nitrogen at an inlet temperature of about −196 ◦ C,
different phases: 
while the expected outlet temperature is around −110 ◦ C. The duration of Phase 1 is related to the time
Phase 1‐Nitrogen: The probes are fed with nitrogen at an inlet temperature of about –196 °C, 
required for the formation of the minimum thickness of the frozen wall (1.5 m);
while the expected outlet temperature is around –110 °C. The duration of Phase 1 is related to the 
Phase 2-Waiting: At the end of nitrogen feeding, in order to obtain a temperature adequate for
time required for the formation of the minimum thickness of the frozen wall (1.5 m); 
brine feeding, and to avoid brine freezing inside the probes;
Phase 2‐Waiting: At the end of nitrogen feeding, in order to obtain a temperature adequate for 
Phase 3-Brine: Maintaining the ice thickness on the tunnel vault by feeding the probes with brine,
brine feeding, and to avoid brine freezing inside the probes; 
at a temperature of about −35 ◦ C. Phase 3 is used during the tunnel excavation, in order to maintain
Phase 3‐Brine: Maintaining the ice thickness on the tunnel vault by feeding the probes with brine, 
the soil temperature below water freezing over time, and the desired thickness of the frozen wall.
at a temperature of about –35 °C. Phase 3 is used during the tunnel excavation, in order to maintain 
For the realization of two tunnels between Line 1 and Line 6 of the underground station in Piazza
the soil temperature below water freezing over time, and the desired thickness of the frozen wall. 
Municipio in Napoli, the AGF technique with the mixed method has been used. The choice to use
For the realization of two tunnels between Line 1 and Line 6 of the underground station in Piazza 
this method is due to the possibility of combining the cryogenic power of nitrogen with the flexibility
Municipio in Napoli, the AGF technique with the mixed method has been used. The choice to use 
and safety of freezing with brine. It consists essentially of making complementary direct and indirect
this method is due to the possibility of combining the cryogenic power of nitrogen with the flexibility 
methods, using the same freezing probes.
and safety of freezing with brine. It consists essentially of making complementary direct and indirect 
An overview of the nitrogen feeding procedure and system used is reported in Figure 1, which
methods, using the same freezing probes. 
reports the truck and tanks at the construction yard, and the loading phase of liquid nitrogen into
An overview of the nitrogen feeding procedure and system used is reported in Figure 1, which 
the tanks.

(a)  (b)  (c) 

Figure 1. (a) Truck and tanks of liquid nitrogen; (b) loading of liquid nitrogen into the tanks; (c) view 
Figure 1. (a) Truck and tanks of liquid nitrogen; (b) loading of liquid nitrogen into the tanks; (c) view
of the construction site. 
of the construction site.

Figure 2 shows an overview of the construction site, where it is possible to see the nitrogen plant 
Figure 2 shows an overview of the construction site, where it is possible to see the nitrogen plant
and 
and the 
the two brine refrigeration 
two brine refrigeration units. 
units.The 
Thenitrogen feeding 
nitrogen feedingsystem 
systemworks 
worksby gravity. 
by gravity.The brine is 
The brine is
refrigerated by one unit, while the second is used as back up in case of the failure of the first one. 
refrigerated by one unit, while the second is used as back up in case of the failure of the first one.
The present case study involves the construction of two tunnels for the connection between Line
1 and Line 6 of the Metro station of Piazza Municipio in Napoli, southern Italy. The soil affected by
the excavation consists of a layer of pozzolana overlaying a bench of tuff. As shown in Figure 3a,
the horizontal distribution of the freezing probes is influenced by the actual development of the two
tunnels, that have a slight curvature. Instead, the freezing probes, for technological reasons, have a
straight distribution along their axis. Section A-A is located at 5.0 m from the Tunnel Boring Machine
(TBM) extraction well, section B-B is in a central position with respect to the tunnels, and section C-C is
located at 5.0 m from Line 6 Station well. Figure 3b shows a cross section of the case study with the
position of the freezing probes, while Figure 4 shows the axonometry of the two tunnels connecting
Line 1 and Line 6.
Energies 2020, 13, 1272 5 of 24
Energies 2020, 13, x FOR PEER REVIEW  5 of 26 

 
Figure 2. Overview of the construction site and view of the nitrogen and brine plant. 

The present case study involves the construction of two tunnels for the connection between Line 
1 and Line 6 of the Metro station of Piazza Municipio in Napoli, southern Italy. The soil affected by 
the excavation consists of a layer of pozzolana overlaying a bench of tuff. As shown in Figure 3(a), 
the horizontal distribution of the freezing probes is influenced by the actual development of the two 
tunnels, that have a slight curvature. Instead, the freezing probes, for technological reasons, have a 
straight distribution along their axis. Section A‐A is located at 5.0 m from the Tunnel Boring Machine 
(TBM) extraction well, section B‐B is in a central position with respect to the tunnels, and section C‐
C is located at 5.0 m from Line 6 Station well. Figure 3(b) shows a cross section of the case study with   
the  position  of  the  freezing  probes,  while  Figure  4  shows  the  axonometry  of  the  two  tunnels 
Figure 2. Overview of the construction site and view of the nitrogen and brine plant. 
Figure 2. Overview
connecting Line 1 and Line 6.    of the construction site and view of the nitrogen and brine plant.

The present case study involves the construction of two tunnels for the connection between Line 
1 and Line 6 of the Metro station of Piazza Municipio in Napoli, southern Italy. The soil affected by 
the excavation consists of a layer of pozzolana overlaying a bench of tuff. As shown in Figure 3(a), 
the horizontal distribution of the freezing probes is influenced by the actual development of the two 
tunnels, that have a slight curvature. Instead, the freezing probes, for technological reasons, have a 
straight distribution along their axis. Section A‐A is located at 5.0 m from the Tunnel Boring Machine 
(TBM) extraction well, section B‐B is in a central position with respect to the tunnels, and section C‐
C is located at 5.0 m from Line 6 Station well. Figure 3(b) shows a cross section of the case study with 
the  position  of  the  freezing  probes,  while  Figure  4  shows  the  axonometry  of  the  two  tunnels 
connecting Line 1 and Line 6.     
(a)  (b) 
Figure 3. (a) Key plan of the tunnels and freezing probes, with indication of the three cross sections 
Energies 2020, 13, x FOR PEER REVIEW 
Figure 3. (a) Key plan of the tunnels and freezing probes, with indication of the three cross sections6 of 26 
considered; (b) cross section of the case study with the position of the freezing probes. 
considered; (b) cross section of the case study with the position of the freezing probes.

 
(a)  (b) 
Figure 3. (a) Key plan of the tunnels and freezing probes, with indication of the three cross sections 
considered; (b) cross section of the case study with the position of the freezing probes. 

 
Figure 4. Axonometry of the two tunnels connecting Line 1 and Line 6.
Figure 4. Axonometry of the two tunnels connecting Line 1 and Line 6. 

In order to construct the two tunnels connecting Line 1 and Line 6, the excavation of the odd 
tunnel occurred before the one in the even tunnel. The odd tunnel excavation began after freezing 
the soil around it, by cutting the diaphragm of the Line 6 station well (see Figure 3), at the opposite 
point of nitrogen input into the probes, and continuing with the excavation one meter at the time, 
with the laying of steel ribs and spritz beton until the TBM extraction well was reached. Once the odd 
Energies 2020, 13, 1272   6 of 24
Figure 4. Axonometry of the two tunnels connecting Line 1 and Line 6. 

In order to construct the two tunnels connecting Line 1 and Line 6, the excavation of the odd
 
tunnel occurred before the one in the even tunnel. The odd tunnel excavation began after freezing the
In order to construct the two tunnels connecting Line 1 and Line 6, the excavation of the odd 
soil around it, by cutting the diaphragm of the Line 6 station well (see Figure 3), at the opposite point
tunnel occurred before the one in the even tunnel. The odd tunnel excavation began after freezing 
of nitrogen input into the probes, and continuing with the excavation one meter at the time, with the
the soil around it, by cutting the diaphragm of the Line 6 station well (see Figure 3), at the opposite 
laying of steel ribs and spritz beton until the TBM extraction well was reached. Once the odd tunnel
point of nitrogen input into the probes, and continuing with the excavation one meter at the time, 
excavation operations were completed, the even tunnel was frozen and constructed.
with the laying of steel ribs and spritz beton until the TBM extraction well was reached. Once the odd 
On the exterior side of the tunnels, two jet grouting walls reaching the depth of the tuff bench had
tunnel excavation operations were completed, the even tunnel was frozen and constructed. 
been employed with the purpose of containing fluids motion in the ground. The freezing of the tunnels
On the exterior side of the tunnels, two jet grouting walls reaching the depth of the tuff bench 
has been realized using 43 freezing probes (23 for the odd tunnel, and 20 for the even one) with a length
had been employed with the purpose of containing fluids motion in the ground. The freezing of the 
oftunnels has been realized using 43 freezing probes (23 for the odd tunnel, and 20 for the even one) 
about 40 m, arranged in an arch outside the excavation section with a constant wheelbase equal to
0.75 m.a  The
with  length  of  about 
installed 40  m, are
probes arranged 
made upin  an  arch concentric
of two outside  the tubes,
excavation  section 
as shown with  a 5.
in Figure constant 
The outer one
iswheelbase equal to 0.75 m. The installed probes are made up of two concentric tubes, as shown in 
made of steel and has a diameter equal to 76 mm, while the inner one is made of copper, with a
Figure 5. The outer one is made of steel and has a diameter equal to 76 mm, while the inner one is 
diameter of 28 mm.
made of copper, with a diameter of 28 mm. 

 
Figure 5. Freezing probe configuration and refrigerant direction: (left) scheme of freezing probe and
Figure 5. Freezing probe configuration and refrigerant direction: (left) scheme of freezing probe and 
temperature profile along the probes; (right) picture of the actual freezing probe.
temperature profile along the probes; (right) picture of the actual freezing probe. 

3. Mathematical Model of the Freezing Process

The model is based on 2D conductive heat transfer in the ground surrounding the probes, and takes
into account the phase change phenomenon of water. The mathematical model has been implemented
within the commercial software Comsol Multiphysics, based on finite element discretization technique,
and has been solved by using the MUltifrontal Massively Parallel sparse direct Solver. The ground
subdomain has a depth of 20 m and a length of 35 m, and can be considered sufficiently large to avoid
thermal interference with the external environment, and sufficiently deep to assume an undisturbed
soil temperature. Figure 6 shows the cross section of the computational domain considered in the
present analysis.
The assumptions underlying the present model are the following: (i) homogeneous and isotropic
materials in each layer of the computational domain; (ii) thermophysical properties of the soil varying
with temperature, between the frozen and unfrozen phases; (iii) for the whole volume of soil, phase
transition takes place at a temperature of 0 ◦ C within an interval of 1 ◦ C; (iv) the temperature of the
cooling fluid in the probes varies linearly along the axis; (v) heat transfer is purely conductive in the
soil, due to the limited convective motion of the water in the ground.
domain considered in the present analysis. 
The assumptions underlying the present model are the following: (i) homogeneous and isotropic 
materials in each layer of the computational domain; (ii) thermophysical properties of the soil varying 
with temperature, between the frozen and unfrozen phases; (iii) for the whole volume of soil, phase 
transition takes place at a temperature of 0 °C within an interval of 1 °C; (iv) the temperature of the 
Energies 2020, 13, 1272 7 of 24
cooling fluid in the probes varies linearly along the axis; (v) heat transfer is purely conductive in the 
soil, due to the limited convective motion of the water in the ground. 

 
6. Sketch
Figure 6.
Figure  Sketch ofof 
the computational
the  domain
computational  withwith 
domain  the two
the tunnels considered
two  tunnels  in the numerical
considered  analysis.
in  the  numerical 
analysis. 
3.1. Governing Equations
3.1. Governing equations 
The problem under investigation has been simulated by means of a dynamic model reproducing
the 2D conductive heat transfer in the ground [15,16], taking into account the phase change of water in
The problem under investigation has been simulated by means of a dynamic model reproducing 
the soil. The governing equation for heat transfer is reported as follows:
the 2D conductive heat transfer in the ground [15, 16], taking into account the phase change of water 
Transient conduction heat transfer
in the soil. The governing equation for heat transfer is reported as follows: 
Transient conduction heat transfer 
∂T ∂ ∂T
!
∂ ∂T
!
ρi cpi = ki + ki +Q i = 1, . . . . . . . . . ., n   (1)
𝜕𝑇 𝜕∂t 𝜕𝑇 ∂x ∂x 𝜕 ∂y ∂y
𝜕𝑇
𝜌𝑐 𝑘 𝑘 𝑄 𝑖 1, … … … . , 𝑛  (1) 
𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦
where ρi is the density of the materials constituting the subdomain (kg/m ), cpi is the specific heat
 
3
where 
capacityρi  is  the  density 
(J/kg·K), of  the 
ki is the materials constituting 
thermal conductivity (W/m·K), the  subdomain  (kg/m3), cpi(K),
T is the temperature   is  the 
andspecific 
finally, heat 
Q is
capacity (J/kg∙K), k
the heat generatedi is the thermal conductivity (W/m∙K), T is the temperature (K), and finally, Q is 
or absorbed per unit of volume (W/m3 ). This last term on the right-hand side
the  heat 
allows usgenerated 
to model theor  absorbed 
latent heatper  unit  of  volume 
of solidification, (W/m
such
3).  This  last  term  on  the  right‐hand  side 
as the heat absorbed or released at a constant
allows us to model the latent heat of solidification, such as the heat absorbed or released at a constant 
temperature during the phase change of water in the soil. In fact, this phenomenon is characterized by
temperature during the phase change of water in the soil. In fact, this phenomenon is characterized 
a significant variation of the thermal diffusion coefficient and the specific heat of the saturated soil, in
by a significant variation of the thermal diffusion coefficient and the specific heat of the saturated 
addition to the absorption of melting latent heat. The time-step employed in the simulations is equal
soil, in addition to the absorption of melting latent heat. The time‐step employed in the simulations 
to twelve hours.
is equal to twelve hours.   
 3.2. Phase Change in the Soil
 
The formulation used in the present work provides the latent heat as an additional term in the
heat capacity.
Instead of adding the latent heat L in the energy balance equation exactly when the material
reaches its phase change temperature, Tpc, it is assumed that the transformation occurs in a temperature
interval between Tpc − ∆T/2 and Tpc + ∆T/2. ∆T is the temperature interval which occurs within
the phase change of water. In this interval, the material phase is modeled by a smooth function, ϑ,
representing the fraction of phase change during transition, which is equal to 1 below Tpc − ∆T/ 2 and
to 0 above Tpc + ∆T ⁄ 2. The density, ρ, and the specific enthalpy, h, of the ground are then calculated as:

ρ = ϑρphase1 + (1 − ϑ)ρphase2 (2)

h = ϑρphase1 hphase1 + (1 − ϑ)ρphase2 hphase2 (3)

Energies 2020, 13, 1272 8 of 24

where phase1 and phase2 indicate the characteristics of the material during the different phases of water
within the soil. The specific heat at constant pressure can be defined as:

∂h
!
cp = (4)
∂T p

that becomes, with the product derivatives:

1    ∂αm
cp = ϑ1 ρphase1 cp,phase1 + ϑ2 ρphase2 cp,phase2 + hphase2 − hphase1 (5)
ρ ∂T

where ϑ1 and ϑ2 are, respectively, equal to ϑ and 1 − ϑ. The term αm is defined as:

1 ϑ2 ρphase2 − ϑ1 ρphase1
αm = (6)
2 ρ

and it is assumed equal to −1/2, before the phase change process, and 1/2 at the end of the process.
Therefore, the specific heat during the phase change phenomenon is given by the sum of two
terms, one proportional to the equivalent thermal capacity Ceq :

1 
ceq = ϑ1 ρ f 1 Cp, f 1 + ϑ2 ρ f 2 cp, f 2 (7)
ρ

and the other proportional to the latent heat CL :

  dαm dαm
cL ( T ) = h f 2 − h f 1 =L (8)
dT dT
so that the total heat per unit of volume released during the phase change process is equal to the latent
heat of solidification:
Z T f 1 + ∆T Z T f 1 + ∆T
2 2 dα
m
Q= CL (T )dT = L dT (9)
T ∆T T ∆T
dT
f 2+ f 2+
2 2

Finally, the apparent thermal capacity Cp used in the heat conservation equation, is given by:

1 
cp = ϑ1 ρphase1 cp,phase1 + ϑ2 ρphase2 cp,phase2 + cL (10)
ρ

The effective thermal conductivity of the portion of soil affected by the phase change is expressed as:

k = ϑ1 kphase1 + ϑ2 kphase2 (11)

while the effective density is calculated as:

ρ = ϑ1 ρphase1 + ϑ2 ρphase2 (12)

Finally, continuity of heat flux is assumed on internal interfaces between the materials. To solve
the equation of transient heat conduction, appropriate values must be assigned to the coefficients ρi ,
cpi , ki . The values used in this work have been derived from the literature (Papakonstantinou et al. [6]
and Rocca [3]). The thermal characteristics as mineral density, dry density, porosity, wet density of the
soil layers and jet-grouting, and the thermal characteristics dependent on the frozen and unfrozen
phase, as thermal conductivity and heat capacity, are reported in Table 1.
Energies 2020, 13, 1272 9 of 24

Table 1. Characteristics of the soil layers and jet-grouting ([3,6]).

Property Tuff Pozzolana Jet-Grouting

Porosity, n 0.5 0.51
Mineral density, ρs (kg/m3 ) 2713 2392 3000
Dry density, ρd (kg/m3 ) 1223 1172
Wet density, ρwet (kg/m3 ) 1733 1682
Unfrozen/frozen ground properties
(at 16 ◦ C/−50 ◦ C)
Thermal conductivity, k (W/mK) 1.48/3.14 1.28/2.61 1.40
Heat capacity, cv (kJ/m3 K) 3120/1990 3150/2709 900

3.3. Initial and Boundary Conditions

The initial condition in the whole domain is:

T (x, y, 0) = T0 = 16 ◦ C ∀(x, y)  Ω (13)

where Ω is the computational domain for each of the three sections considered in this work and
reported in Figure 3.
The following Dirichlet condition is imposed on the external surface of each probe during Phase 1
and Phase 3 of the AGF process:

T (x, y, ϑ) = T (ϑ) ∀(x, y)  Ωprobes (14)

The boundary conditions employed in the present model refer to the soil and probes domain
and are sketched in Figures 7 and 8. In particular, Figure 7 refers to the domain considered before
excavation of the first tunnel, while Figure 8 refers to the domain considered after excavation of the
first tunnel. The temperature of the top, bottom and lateral surfaces of the soil has been assumed to be
constant during the analysis, equal to the average yearly temperature of the site under investigation,
Energies 2020, 13, x FOR PEER REVIEW  10 of 26 
T0 = 16 ◦ C.
As previously specified, due to the complexity of nitrogen phase change phenomena occurring
During “Phase 3‐Brine”, the temperature of the probe boundary has been imposed equal to the 
inside the freezing probes, a linear temperature profile has been assumed for the refrigerant fluid
temperature of the brine, –33 °C, in all the sections of the excavation. 
between the inlet and outlet sections of the probes (refer to Figure 5).
 

 
Figure 7. Computational domain and boundary conditions before the excavation.
Figure 7. Computational domain and boundary conditions before the excavation. 
  
Energies 2020, 13, 1272 10 of 24
Figure 7. Computational domain and boundary conditions before the excavation. 

 
Figure 8. Computational domain and boundary conditions after the excavation of the odd tunnel.
Figure 8. Computational domain and boundary conditions after the excavation of the odd tunnel. 
The temperature boundary conditions applied on the external perimeter of each probe depend on
In order to simulate the excavation of the even tunnel, the same steps employed for that of the 
the phase of the freezing process. During “Phase 1-Nitrogen”, the temperature has a linear variation
odd  tunnel  have  been  considered.  During  these  phases,  the  excavated  odd  tunnel  (left)  has  been 
along the axis, from −196 ◦ C to −110 ◦ C, as shown in Figure 5. During “Phase 2-Waiting”, the adiabatic
reproduced  by  eliminating  the  corresponding  domain  of  soil  and  applying  a  proper  boundary 
condition, ∇T·n = 0, has been imposed on the probe boundary.
condition (refer to Figure 8). This condition takes into account the presence of men, vehicles and air 
During “Phase 3-Brine”, the temperature of the probe boundary has been imposed equal to the
circulation in the excavated tunnel, and is represented by convective heat transfer on the walls of the 
temperature of the brine, −33 ◦ C, in all the sections of the excavation.
odd tunnel: 
In order to simulate the excavation of the even tunnel, the same steps employed for that of the odd
tunnel have been considered. During 𝑘∇these
𝑇 ∙ 𝑛 phases,
ℎ 𝑇 the𝑇 excavated
  (16) 
odd tunnel (left) has been reproduced
by eliminating the correspondingℎ 𝑥, domain
𝑦, 𝜗 ℎof soil15and
𝑊/𝑚applying
𝐾  a proper boundary condition (refer
(17)  to
Figure 𝑇  This
where  8). condition takes into account the presence of men, vehicles and air circulation in the
is equal to 30 °C. 
excavated
  tunnel, and is represented by convective heat transfer on the walls of the odd tunnel:

− k∇T·n = h(T − T∞ ) (15)

h(x, y, ϑ) = h∞ = 15 W/m2 K (16)

where T∞ is equal to 30 ◦ C.

3.4. Mesh Sensitivity Analysis

A mesh sensitivity analysis has been carried out in order to obtain grid-independent numerical
results. A domain of 20 × 25 m2 has been considered. All the computational grids are made by
triangular quadratic elements, and are refined near the freezing probes (Figure 9).
Table 2 reports the details of the eight grids considered, together with a summary of the main
numerical results. In particular, considering the nitrogen activation (phase 1), the days required for the
formation of the frozen wall at 1.5 m have been calculated and reported in Table 2, together with the
computing time necessary to reach the convergence. The nitrogen freezing phase is stopped when
the desired design value of the frozen wall thickness is reached. The evolution of the frozen wall
can be monitored, controlling the temperature of the ground at 0.5 m from the probe axis, which
generally must be around −10 ◦ C. Therefore, in the analysis, the temperature at this point has been
taken into account.
On the basis of the present sensitivity analysis, the grid employed for the calculations is the one
with 103,482 elements (letter f in the table), since the difference between the results obtained by using
this grid, and the ones obtained by using the most effective adaptive grid, is around 1%.
 Energies 2020, 13, x FOR PEER REVIEW  11 of 26 
Energies 2020, 13, 1272 11 of 24

3.4. Mesh sensitivity analysis 

   
(a)  (b) 

   
(c)  (d) 

   
(e)  (f) 

   
(g)  (h) 
Figure 9. Meshes used for the grid independent study (a) 5339 elements; (b) 9596 elements;
(c) 42,671 elements; (d) 73,804 elements; (e) 82,046 elements; (f) 103,482 elements; (g) 113,342 elements;
(h) 177,964 elements.
Energies 2020, 13, 1272 12 of 24

Table 2. Details of the meshes and obtained results for the grid-independent study.

Number of Time for Frozen Wall

Temperature at 0.5 m Computing Time
Mesh Code Elements Formation at 0.6 m
(◦ C) (min)
(-) (d)
a 5339 9 −13.1 18
b 9596 9 −14.6 33
c 42,671 8 −12.3 399
d 73,804 8 −12.0 1231
e 82,046 8 −11.3 1820
f 103,482 8 −11.4 4213
g 113,342 8 −11.6 8122
h 177,964 8 −11.6 25,015

4. Model Validation
The present model has been validated against the numerical data reported by Colombo [5], which
is validated with on-field data. In that study, the software ABAQUS was used to solve the numerical
model, and the computational domain was meshed by using DC2D4 elements. The geometry of the
tunnel was symmetric. The present study reproduces the case study of Colombo [5], based on a 2D
model, by taking into account a computational domain representative of a portion of land equal to
10 × 20 m. A mesh consisting of 33,605 triangular elements has been considered, and the thermal
characteristics of the materials are those reported in Colombo [5], in particular: volumetric heat capacity
of 1910 kJ/m3 K and 3100 kJ/m3 K for solid and liquid phases, respectively; thermal conductivity of
3.07 W/mK and 1.48 W/mK for the solid and liquid phases, respectively; volumetric latent heat of
179,280 kJ/m3 and a saturated tuff density of 1550 kg/m3 . The initial temperature has been assumed
equal to 18 ◦ C, and the analysis has been carried out, imposing a linear variation of temperature down
to −33 ◦ C for the first day, on the nodes representing the perimeter of the freezing probe. In order
to compare the results obtained from the present FEM analysis with those reported in Colombo [5],
the authors have considered two points located on a line orthogonal to the junction between the probes,
at 0.50 m and 0.90 m. Figure 10 shows the conditions that determine the propagation velocity of the
freezing front in the tuff, evaluated for two points located at 0.5 m and 0.9 m from the freezing probe.
It is evident that the temperature gradually decreases over time, and that a good agreement between
the present numerical results and those available in the literature [5] is observed.
Energies 2020, 13, x FOR PEER REVIEW 13 of 26

21 days
Figure 1. Mesh and development of the frozen wall for model validation.

1 day 5 days

Figure 10. Cont.

 21 days
Figure 1. Mesh and development of the frozen wall for model validation.
Energies 2020, 13, 1272 13 of 24

1 day 5 days

10. Comparison
FigureFigure between
10. Comparison field
between fieldemission
emission microscopy (FEM)
microscopy (FEM) analysis
analysis and and Colombo
Colombo [5]. [5].

5. Results and Discussion

5.1. Odd Tunnel Freezing

5. Results and discussion
The numerical analysis has been developed considering the probes connected in parallel, i.e., all with
5.1. Odd tunnel freezing
the same temperature distribution. The boundary conditions on the surface of probes, described in
The numerical analysis has been developed considering the probes connected in parallel, i.e. all
Section 3, are summarized in Table 3, together with the calculated thickness of frozen wall after 14 days.
with the same temperature distribution. The boundary conditions on the surface of probes, described
in Section 3, are summarized in Table 3, together with the calculated thickness of frozen wall after 14
Table 3. Boundary conditions imposed on the surface of the probes at the considered sections and
days.
for the different freezing phases, and frozen wall thickness calculated after 14 days of AGF with
liquid nitrogen.
Table 3. Boundary conditions imposed on the surface of the probes at the considered sections and for
the different freezing phases, and frozen wall thickness calculated after 14 days of AGF with liquid
Phase 1 Phase 2 Phase 3
nitrogen.
Probe Probe Probe Thickness of
z Duration Phase 1 Duration Phase 2 Duration Phase 3
Section Temperature Temperature Temperature frozen wall
(m) (days) (days) (days)
(◦ C) (◦ C) (◦ C) (m)
A-A 5 14 −180 1 Adiabatic 15 −33 1.50
B-B 20 14 −160 1 Adiabatic 15 −33 1.70
C-C 35 14 −120 1 Adiabatic 15 −33 1.90

The results refer to the freezing process in the odd tunnel (left), for which it has been assumed
duration of the direct freezing phase equal to 14 days. The waiting phase, before switching to
the indirect method phase with brine, has been imposed equal to 1 day, enough to have suitable
temperatures for the brine intake in the pipe system, without freezing it. The third phase of brine
retention has been considered to last for 15 days, for a total of 30 days of the whole freezing process in
the odd gallery.
Figure 11 reports the temperature field calculated after 14 days of nitrogen feeding in the three
sections considered for the odd tunnel (refer to Figure 3), and it is possible to clearly see the formation
of the frozen wall, since the figure reports only the temperature values below −2 ◦ C.
In section C-C, the frozen wall reaches the desired thickness of 1.5 m in less time than the other
two sections, due to the temperature distribution in the probes. The thickness of 1.5 m ensures the
static stability of the ground, according to the design specifications. After 14 days of freezing with
nitrogen, the frozen wall reaches the thickness reported in Table 3 for the three sections.
 indirect  method  phase  with  brine,  has  been  imposed  equal  to  1  day,  enough  to  have  suitable 
temperatures for the brine intake in the pipe system, without freezing it. The third phase of brine 
retention has been considered to last for 15 days, for a total of 30 days of the whole freezing process 
in the odd gallery. 
Figure 11 reports the temperature field calculated after 14 days of nitrogen feeding in the three 
Energies 2020, 13, 1272 14 of 24
sections considered for the odd tunnel (refer to Figure 3), and it is possible to clearly see the formation 
of the frozen wall, since the figure reports only the temperature values below –2 °C. 

 
(a) Section A‐A.  (b) Detail 

 
(c) Section B‐B. 

Energies 2020, 13, x FOR PEER REVIEW  15 of 26 

(d) Section C‐C. 
 

. 

(e) Detail of the mesh employed in the simulations. 
Figure 11. from (a) to (d): Temperature field below –2 °C after 14 days of AGF with liquid nitrogen in 
Figure 11. From (a–d): Temperature field below −2 ◦ C after 14 days of AGF with liquid nitrogen in the
the three sections of the odd tunnel reported in Figure 3; (e) detail of the employed mesh. 
three sections of the odd tunnel reported in Figure 3; (e) detail of the employed mesh.

 
In section C‐C, the frozen wall reaches the desired thickness of 1.5 m in less time than the other 
two sections, due to the temperature distribution in the probes. The thickness of 1.5 m ensures the 
static stability of the ground, according to the design specifications. After 14 days of freezing with 
nitrogen, the frozen wall reaches the thickness reported in Table 3 for the three sections. 
 the three sections of the odd tunnel reported in Figure 3; (e) detail of the employed mesh. 

 
In section C‐C, the frozen wall reaches the desired thickness of 1.5 m in less time than the other 
two sections, due to the temperature distribution in the probes. The thickness of 1.5 m ensures the 
Energies 2020, 13, 1272 15 of 24
static stability of the ground, according to the design specifications. After 14 days of freezing with 
nitrogen, the frozen wall reaches the thickness reported in Table 3 for the three sections. 
Figure 12(a) shows the frozen wall thickness reached after 14 days of nitrogen freezing phase, 
Figure 12a shows the frozen wall thickness reached after 14 days of nitrogen freezing phase,
by showing the isotherms at –2 °C, and reports a detail of the frozen wall for section A‐A. Figure 
by showing the isotherms at −2 ◦ C, and reports a detail of the frozen wall for section A-A. Figures 13b
13(b)  and  Figure  12(c)  show  the  frozen  wall  thickness  reached  after  14  days  of  nitrogen  freezing 
and 12c show the frozen wall thickness reached after 14 days of nitrogen freezing phase, by showing
phase, by showing the isotherms at –2 °C, for sections B‐B and C‐C. Moreover, the figure presents a 
the isotherms at −2 ◦ C, for sections B-B and C-C. Moreover, the figure presents a vertical segment
vertical segment across probe n. 13, that is considered useful for the analysis of the temperature field 
across probe n. 13,  that is considered useful for the analysis of the temperature field in the frozen wall.
in the frozen wall. 

 
Energies 2020, 13, x FOR PEER REVIEW  16 of 26 
(a)   

 
(b)   

 
(c)   

Figure 12. Thickness of the frozen wall represented through isotherm at –2 °C, for section A‐A (a), 
Figure 12. Thickness of the frozen wall represented through isotherm at −2 ◦ C, for section A-A (a),
section B‐B (b); section C‐C (c). 
section B-B (b); section C-C (c).

Figure  13  shows  the  temperature  profile  calculated  on  a  vertical  segment  of  1.5  meters  from 
probe 13 (point a) to point b for the three sections considered, as shown in the previous figure. From 
the analysis of the figure, it is evident that at the end of the first phase, the frozen wall has reached 
the minimum required thickness of 1.5m.   
Energies 2020, 13, 1272 16 of 24

Figure 13 shows the temperature profile calculated on a vertical segment of 1.5  m from probe
13 (point a) to point b for the three sections considered, (c)    as shown in the previous figure. From the
analysis of the figure, it is evident that at the end of the first phase, the frozen wall has reached the
Figure 12. Thickness of the frozen wall represented through isotherm at –2 °C, for section A‐A (a), 
minimumsection B‐B (b); section C‐C (c). 
required thickness of 1.5 m.
Figure 14 shows the temperature profile and thickness of the frozen wall with respect to the probe
axis, afterFigure 
14 days13 of
shows  the  temperature 
freezing, for sectionsprofile  calculated 
A-A, B-B on  a 
and C-C. Thevertical 
graph segment 
showsof  1.5 the
that meters  from 
frozen wall at
probe 13 (point a) to point b for the three sections considered, as shown in the previous figure. From 
section A-A after 14 days of activation with nitrogen has a thickness of 1.5 m for section A-A, 1.7 m for
the analysis of the figure, it is evident that at the end of the first phase, the frozen wall has reached 
section B-B, and 1.9 m for section C-C.
the minimum required thickness of 1.5m.   

 
Figure 13. Temperature profile after 14 days of artificial ground freezing (AGF) with liquid nitrogen
Figure 13. Temperature profile after 14 days of artificial ground freezing (AGF) with liquid nitrogen 
along the segment a-b. along the segment a‐b. 

For the odd14 

Figure  tunnel
shows (left), further segments
the  temperature  are considered,
profile  and  in order
thickness  of  the  frozen towall 
verify the
with  thickness
respect  to  the of the
ice wall at the end of Phase 1. Figure 15a shows the cut lines taken into account, while Figure 15b
probe axis, after 14 days of freezing, for sections A‐A, B‐B and C‐C. The graph shows that the frozen 
shows wall at section A‐A after 14 days of activation with nitrogen has a thickness of 1.5 m for section A‐A, 
the temperature profiles calculated after 14 days for section C-C.
1.7 m for section B‐B, and 1.9 m for section C‐C. 
Figure 16a presents the temperature profiles at the sections A-A, B-B, and C-C. From the analysis
of this figure, it is possible to notice that the section located 5 m from the TBM well (A-A) is the most
disadvantaged for the freezing process. This is due to the nitrogen freezing temperature imposed on the
probe perimeter in this section (−120 ◦ C), which is higher than that in the other sections. Consequently,
more time is needed for the formation of the frozen wall in this section, defining the total duration
of the freezing Phase 1 (nitrogen activation). After almost 5 days of activation with brine, constant
temperature values can be reached, but always below the safety value of −5 ◦ C.
Moreover, it is possible to notice that the ground reaches the temperature of 0 ◦ C after about
ten days of liquid nitrogen feeding. This finding is in good agreement with the results reported by
Colombo [17] and Manassero [18]. Moreover, from the analysis of Figure 16a, it is possible to notice
that the defrosting phase for the odd tunnel starts after 60 days and that after 90 days, the temperature
rises up to 0 ◦ C. These temperature values are important to evaluate the ground displacement.
Figure 16b shows the temperature profile for the freezing process of the even tunnel. The beginning
of the freezing with nitrogen starts at the end of the freezing of the odd tunnel.
Table 4 shows the thickness of the frozen wall after 14 days of artificial ground freezing with
nitrogen and is always larger than the design value of 1.5 m for both tunnels.
 Energies 2020, 13, x FOR PEER REVIEW  17 of 26 

wall at section A‐A after 14 days of activation with nitrogen has a thickness of 1.5 m for section A‐A, 
Energies 2020, 13, 1272 17 of 24
1.7 m for section B‐B, and 1.9 m for section C‐C. 
 

 
(a) 

 
Energies 2020, 13, x FOR PEER REVIEW  (b)  18 of 26 

 
(c) 
Figure 14. Temperature profile after 14 days of AGF with liquid nitrogen on segment a‐b at Section 
Figure 14. Temperature profile after 14 days of AGF with liquid nitrogen on segment a-b at Section
A‐A (a); section B‐B (b); section C‐C (c). 
A-A (a); section B-B (b); section C-C (c).

 
For the odd tunnel (left), further segments are considered, in order to verify the thickness of the 
ice wall at the end of Phase 1. Figure 15(a) shows the cut lines taken into account, while Figure 15(b) 
shows the temperature profiles calculated after 14 days for section C‐C. 
Energies 2020, 13, 1272 18 of 24

Table 4. Thickness of the frozen wall after 14 days of AGF with nitrogen liquid.

Odd Tunnel z Thickness Frozen wall (m)

Section (m) a-a’ b-b’ c-c’ d-d’ e-e’ f-f’ g-g’
C-C 35 2.0 2.1 2.1 2.1 2.0 2.1 2.0
Even Tunnel z Thickness Frozen wall (m)
Section (m) a-a’ b-b’ c-c’
A-A 5 1.8 1.6 1.6
B-B 20 1.9 1.9 1.9
C-C 35 2.1 2.1 2.1
Energies 2020, 13, x FOR PEER REVIEW  19 of 26 

 
(a) 

 
(b) 
Figure 15. (a) Cutline for the section C‐C; (b) Temperature trend after 14 days for section C‐C and 
Figure 15. (a) Cutline for the section C-C; (b) Temperature trend after 14 days for section C-C and
frozen wall was formed for all cut lines considered. 
frozen wall was formed for all cut lines considered.
Figure  16(a)  presents  the  temperature  profiles  at  the  sections  A‐A,  B‐B,  and  C‐C.  From  the 
analysis of this figure, it is possible to notice that the section located 5 m from the TBM well (A‐A) is 
the most disadvantaged for the freezing process. This is due to the nitrogen freezing temperature 
imposed  on  the  probe  perimeter  in  this  section  (–120  °C),  which  is  higher  than  that  in  the  other 
sections.  Consequently,  more  time  is  needed  for  the  formation  of  the  frozen  wall  in  this  section, 
defining  the  total  duration  of  the  freezing  Phase  1  (nitrogen  activation).  After  almost  5  days  of 
activation with brine, constant temperature values can be reached, but always below the safety value 
 days  of  liquid  nitrogen  feeding.  This  finding  is  in  good  agreement  with  the  results  reported  by 
Moreover, it is possible to notice that the ground reaches the temperature of 0 °C after about ten 
Colombo [17] and Manassero [18]. Moreover, from the analysis of Figure 16(a), it is possible to notice 
days  of  liquid  nitrogen  feeding.  This  finding  is  in  good  agreement  with  the  results  reported  by 
that the defrosting phase for the odd tunnel starts after 60 days and that after 90 days, the temperature 
Colombo [17] and Manassero [18]. Moreover, from the analysis of Figure 16(a), it is possible to notice 
rises up to 0 °C. These temperature values are important to evaluate the ground displacement. 
that the defrosting phase for the odd tunnel starts after 60 days and that after 90 days, the temperature 
Figure  16(b)  shows  the  temperature  profile  for  the  freezing  process  of  the  even  tunnel.  The 
rises up to 0 °C. These temperature values are important to evaluate the ground displacement. 
Energies 2020, 13, 1272
beginning of the freezing with nitrogen starts at the end of the freezing of the odd tunnel.  19 of 24
Figure  16(b)  shows  the  temperature  profile  for  the  freezing  process  of  the  even  tunnel.  The 
beginning of the freezing with nitrogen starts at the end of the freezing of the odd tunnel. 

Freezing odd tunnel Freezing even tunnel

 
Freezing odd tunnel Freezing even tunnel
(a)   

  (a) 

Freezing even tunnel  

Freezing
(b)  even tunnel  

Figure 16. (a) Temperature profiles at the control (b) 

point on the top of the odd tunnel for the three
sections; (b) Temperature profiles at the control point on the top of the even tunnel for the three sections.

5.2. Even Tunnel Freezing

At the end of the 60 days, after the freezing and excavation activities of the odd tunnel, freezing
is repeated with the same conditions for the even tunnel, disabling completely the supply of brine
for the maintenance phase of the odd tunnel. The temperature applied on the freezing probes for the
three freezing phases is reported in Table 3. The duration of each freezing phase of the even tunnel is
equal to fourteen days for Phase 1 with nitrogen, one day for Phase 2, and forty days for Phase 3 (brine
maintenance), for a total of sixty days. The results obtained at the end of the first phase with nitrogen
for the even tunnel have been reported in Figure 17.
is repeated with the same conditions for the even tunnel, disabling completely the supply of brine for 
the maintenance phase of the odd tunnel. The temperature applied on the freezing probes for the 
three freezing phases is reported in Table 3. The duration of each freezing phase of the even tunnel 
is equal to fourteen days for Phase 1 with nitrogen, one day for Phase 2, and forty days for Phase 3 
(brine maintenance), for a total of sixty days. The results obtained at the end of the first phase with 
Energies 2020, 13, 1272 20 of 24
nitrogen for the even tunnel have been reported in Figure 17. 
 
 
 
 
 
 
 
 
 
 
 
 
  Energies 2020, 13, x FOR PEER REVIEW  22 of 26 
(a) Section A‐A 
 
 

 
(b) Section B‐B 

 
(c) Section C‐C 

Figure 17. Frozen wall at 0 °C after 14 days. 
Figure 17. Frozen wall at 0 ◦ C after 14 days.

The results of artificial freezing of soil with the use of liquid nitrogen, after fourteen days from 
The results of artificial freezing of soil with the use of liquid nitrogen, after fourteen days from the
the beginning of the phase, are presented in Figure 18. In particular, the thickness of the frozen wall 
beginning of the phase, are presented in Figure 18. In particular, the thickness of the frozen wall in the
in the three sections is reported, identified by the isotherm at –2 °C. The minimum thickness of 1.5 m 
three sections is reported, identified by the isotherm at −2 ◦ C. The minimum thickness of 1.5 m is used
as is used as a design target necessary to consider Phase 1 completed. 
a design target necessary to consider Phase 1 completed.
Energies 2020, 13, 1272 21 of 24
Energies 2020, 13, x FOR PEER REVIEW 23 of 26

(a)

(b)

(c)

Figure 18. Cont.

Energies 2020, 13, 1272 22 of 24
Energies 2020, 13, x FOR PEER REVIEW  24 of 26 

 
(d)   
Figure 18. The frozen wall after 14 days of nitrogen feeding for section A‐A (a); for section B‐B (b); for 
Figure 18. The frozen wall after 14 days of nitrogen feeding for section A-A (a); for section B-B (b); for
section C‐C (c); Temperature profile for section C‐C (d). 
section C-C (c); Temperature profile for section C-C (d).

Figure 18d shows the temperature profile at three cut lines of section C‐C during the freezing 
Figure 18d shows the temperature profile at three cut lines of section C-C during the freezing
process of the the 
process  of  eveneven  tunnel. 
tunnel. FromFrom  the  analysis 
the analysis of thisof  this  figure, 
figure, it  is  possible 
it is possible to observeto that
observe  that  the 
the temperature
values ◦
temperature values below ‐2 °C involve a thickness of ground around 2.1 m, in all the three cut lines 
below −2 C involve a thickness of ground around 2.1 m, in all the three cut lines considered,
considered, as reported in Table 4. This is a relevant result, proving that the minimum ice wall limit 
as reported in Table 4. This is a relevant result, proving that the minimum ice wall limit (1.5 m) is
(1.5 m) is respected in the whole ice vault. 
respected in the whole ice vault.

6. 6. Conclusions 
Conclusions
This work presents an analysis of the heat transfer phenomena occurring during the artificial 
This work presents an analysis of the heat transfer phenomena occurring during the artificial
ground freezing (AGF) process on the actual case study of the tunnels between Line 1 and Line 6 of 
ground freezing (AGF) process on the actual case study of the tunnels between Line 1 and Line 6 of the
the Underground station in Piazza Municipio in Napoli, southern Italy. An efficient numerical model, 
Underground station in Piazza Municipio in Napoli, southern Italy. An efficient numerical model,
based on conductive heat transfer and water phase change, has been developed and validated against 
based on conductive heat transfer and water phase change, has been developed and validated against
the data available in the literature.   
the data available in the literature.
The present model has us allowed to simulate, for the first time in the literature, a mixed‐method 
The present model has us allowed to simulate, for the first time in the literature, a mixed-method
used for the freezing process from the first phase, based on the use of nitrogen, a maintenance phase 
used for the freezing process from the first phase, based on the use of nitrogen, a maintenance phase
that allows us to raise the temperature of the probes, and a third phase that involves the use of brine. 
that allows us to raise the temperature of the probes, and a third phase that involves the use of brine.
The numerical model has allowed us to reproduce the evolution of the temperature field during the 
The numerical
whole  model
excavation  has allowed
process.  us tohas 
The  model  reproduce
been  used theto 
evolution
illustrate of theheat 
the  temperature field during
transfer  phenomena 
theassociated with the phase change and the influence of latent heat, the influence of phase change, and 
whole excavation process. The model has been used to illustrate the heat transfer phenomena
associated with the phase change and the influence of latent heat, the influence of phase change,
temperature variation. 
and temperature variation.
The present modeling activity allows us to identify the possible solutions for reducing the time 
The present modeling activity allows us to identify the possible solutions for reducing the time
required for the completion of the excavation activities and the freezing of the galleries. In addition 
required for the
to  properly  completion
planning  of the excavation
the  nitrogen  activities
feeding  phase,  and the freezing
this  approach  of the
allows  the  galleries.
analysis  In addition
of  alternative 
to solutions 
properly to  accelerate 
planning thethe  freezing feeding
nitrogen of  the  soil, 
phase,such  as approach
this increasing allows
the  number  of  probes, 
the analysis or  using 
of alternative
different configurations. 
solutions to accelerate the freezing of the soil, such as increasing the number of probes, or using
differentThe results show that the time needed to complete the freezing process and excavation of the 
configurations.
two tunnels, consists of 120 days. Instead, 14 days are required to obtain an ice vault with a thickness 
The results show that the time needed to complete the freezing process and excavation of the two
larger than the minimum value required for safety reasons (1.5 m). In particular, the section that first 
tunnels, consists of 120 days. Instead, 14 days are required to obtain an ice vault with a thickness
reaches the formation of the frozen wall is the one at 5 m from the well of Line 6, and from this section, 
larger than the minimum value required for safety reasons (1.5 m). In particular, the section that first
the excavation operations can begin. Moreover, the defrosting phase for the odd tunnel starts after 60 
reaches the formation of the frozen wall is the one at 5 m from the well of Line 6, and from this section,
days, and after 90 days, the temperature rises up to 0 °C.   
the excavation operations can begin. Moreover, the defrosting phase for the odd tunnel starts after
60 days, and after 90 days, the temperature rises up to 0 ◦ C.
Energies 2020, 13, 1272 23 of 24

The authors believe that the present model is useful to optimize the AGF technique, and to better
understand the heat transfer phenomena occurring between probes and ground.

Author Contributions: Conceptualization, N.M. and A.M.; methodology, Eng. G.N.; software, N.M., A.M., G.N.
and P.M.; formal analysis, G.N. and P.M.; resources, N.M.; data curation, G.N.; writing—original draft preparation,
G.N.; writing—review and editing, N.M. and A.M.; supervision, N.M. and A.M.; project administration, N.M.
and F.C.; funding acquisition, N.M. and F.C. All authors have read and agreed to the published version of
the manuscript.
Funding: This research was funded by a research agreement with Ansaldo STS | A Hitachi Group Company.
Acknowledgments: The authors thank the support of Ansaldo STS | A Hitachi Group Company, in particular
Eng. Giuseppe Molisso, MetroTec Scarl, and Trevi SpA.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Russo, G.; Corbo, A.; Cavuoto, F.; Autuori, S. Artificial Ground Freezing to excavate a tunnel in sandy soil.
Measurements and back analysis. Tunn. Undergr. Sp. Technol. 2015, 50, 226–238. [CrossRef]
2. Smith, N. Artificial ground freezing. In Temporary Works: Principles of Design and Construction; Grant, M.,
Pallet, F.P., Eds.; ICE Publishing: London, UK, 2012; pp. 109–116.
3. Rocca, O. Congelamento Artificiale del Terreno; Hevelius Edizioni: Napoli, Italy, 2011.
4. Andersland, O.B.; Ladanyi, B. Frozen Ground Engineering; Wiley: Hoboken, NJ, USA, 2004.
5. Colombo, G. Il congelamento artificiale del terreno negli scavi della metropolitana di Napoli: Valutazioni
teoriche e risultati sperimentali. 2010.
6. Papakonstantinou, S.; Pimentel, E.; Anagnostou, G. Analysis of artificial ground freezing in the Pari-Duomo
platform tunnel of the Naples metro. Numer. Methods Geotech. Eng. 2000, 281–284. [CrossRef]
7. Pimentel, E.; Papakonstantinou, S.; Anagnostou, G. Numerical interpretation of temperature distributions
from three ground freezing applications in urban tunnelling. Tunn. Undergr. Space Technol. 2012, 28, 57–69.
[CrossRef]
8. Vitel, M.; Rouabhi, A.; Tijani, M.; Guerin, F. Modeling heat transfer between a freeze pipe and the surrounding
ground during artificial ground freezing activities. Comput. Geotech. 2015, 63, 99–111. [CrossRef]
9. Vitel, M.; Rouabhi, A.; Tijani, M.; Guérin, F. Modeling heat and mass transfer during ground freezing
subjected to high seepage velocities. Comput. Geotech. 2016, 73, 1–15. [CrossRef]
10. Marwan, A.; Zhou, M.-M.; Abdelrehim, M.Z.; Meschke, G. Optimization of artificial ground freezing in
tunneling in the presence of seepage flow. Comput. Geotech. 2016, 75, 112–125. [CrossRef]
11. Kang, Y.; Liu, Q.; Cheng, Y.; Liu, X. Combined freeze-sealing and New Tubular Roof construction methods
for seaside urban tunnel in soft ground. Tunn. Undergr. Space Technol. 2016, 58, 1–10. [CrossRef]
12. Panteleev, I.; Kostina, A.; Zhelnin, M.; Plekhov, A.; Levin, L. Intellectual monitoring of artificial ground
freezing in the fluid-saturated rock mass. Procedia Struct. Integr. 2017, 5, 492–499. [CrossRef]
13. Alzoubi, M.A.; Madiseh, A.; Hassani, F.P.; Sasmito, A.P. Heat transfer analysis in artificial ground freezing
under high seepage: Validation and heatlines visualization. Int. J. Therm. Sci. 2019, 139, 232–245. [CrossRef]
14. Fan, W.; Yang, P. Ground temperature characteristics during artificial freezing around a subway cross passage.
Transp. Geotech. 2019, 20, 100250. [CrossRef]
15. Colangelo, F.; De Luca, G.; Ferone, C.; Mauro, A. Experimental and numerical analysis of thermal and
hygrometric characteristics of building structures employing recycled plastic aggregates and geopolymer
concrete. Energies 2013, 6, 6077–6101. [CrossRef]
16. Carotenuto, A.; Marotta, P.; Massarotti, N.; Mauro, A.; Normino, G. Energy piles for ground source heat
pump applications: Comparison of heat transfer performance for different design and operating parameters.
Appl. Therm. Eng. 2017, 124, 1492–1504. [CrossRef]
Energies 2020, 13, 1272 24 of 24

17. Colombo, G.; Lunardi, P.; Cavagna, B.; Cassani, G.; Manassero, V. The artificial ground freezing technique
application for the Naples underground. In Proceedings of the World Tunnel Congress, Agra, India,
22–24 September 2008; pp. 1–15.
18. Manassero, V.; Di Salvo, G.; Giannelli, F.; Colombo, G. A Combination of artificial ground freezing and
grouting for the excavation of a large size tunnel below groundwater. In Proceedings of the International
Conference on Case Histories in Geotechnical Engineering, Arlington, VA, USA, 11–16 August 2008; pp. 1–13.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).