**4.2 Quantitative calculation of ice forming and melting**

Accurately, ice forming and melting in caves include at least three physical processes: water flow in porous media (limestone), natural convection of low viscous material (air), and water-ice phase change. The numerical simulation of each process involves complex mathematical method, especially the second physical process. Based on some assumptions, an equivalent method was used to deal with air natural convection [27].

## *4.2.1 Basic ideas of simulation*

Two heat transfer processes must be considered to interpret the existence of ice deposits in Ningwu ice cave, i.e., thermal conduction and convection. The water-ice phase change can reduce the rate of temperature change. So the phase change should be taken into account. The conducting process is governed by conduction equation, which is relatively easy to be computed, while for the convection process, the convection pattern of air and its thermal consequences are hard to determine exactly, because of the complicated geometry of the ice cave and complex varying boundary conditions. In consideration of this, a broadly used, simplified approach is applied in this study: evaluation of Nusselt number (*Nu*) and solving the conductive equation by introducing an equivalent thermal conductivity of the convecting air. For an upright circular tube, the physical relation between the temperature difference of the top and the bottom and *Nu* can be evaluated by adopting the experimental relation of natural convection. The enthalpy method can be adapted to compute the ice-water phase change.

In every computing time step of our modeling process, it is judged whether air convection occurs according to the temperature difference between the top and the bottom of the cave. If there is no convection, the simple conduction problem will be solved, while if convection occurs, an effective conductivity is used in the conduction equation.

### *4.2.2 Equivalent thermal conductivity*

The heat conduction equation and enthalpy method for dealing with water-ice phase change were detailed in our previous work [27]. Equivalent thermal conductivity method was a key point of computing and was shown as follows:

*Nu* is defined as the ratio of convection heat transfer to pure conduction heat transfer under the same conditions. The efficiency of energy transfer is *Nu* times greater than the conductive efficiency under the same conditions. Therefore, an equivalent thermal conductivity can be introduced, which is *Nu* times greater than the air thermal conductivity [28]. *Nu* is related to physical properties (e.g., viscosity and conductivity of air), to the temperature difference of air at the top and the bottom of the cave, and also to the

geometry of the cave.

Ningwu ice cave can be evaluated by an upright circular tube. For an upright circular tube, *Nu* can be calculated based on experimental fluid thermodynamics studies. Once Eq. (1) is satisfied [29, 30] (the case for Ningwu ice cave), the natural convection heat transfer relation [29, 31] can be expressed as Eq. (2):

$$d/h \ge \Im \mathfrak{F}/Gr^{1/4} \tag{1}$$

$$\text{Nu}\_m = \text{C} \left( Gr \cdot Pr \right)\_m^n \tag{2}$$

where *d* and *h* in Eqs. (1) and (2) are the diameter and the height of a circular tube respectively; *Num* is the Nusselt number; *Pr* is the Prandtl number and is dependent only on the material, e.g., *Pr* is 0.7 of air; and *C* and *n* are constants, the values of which are shown in **Table 1**.

*Gr* is the Grashof number:

$$Gr = \text{ g\u0\Delta T l}^3 / \text{o}^2 \tag{3}$$

**93**

*A Review of Chinese Ice Caves*

Eq. (2):

**Table 1.**

*DOI: http://dx.doi.org/10.5772/intechopen.89178*

Gr *number and constant for different flow types [28].*

*Gr* = 1.041 × 1014Δ*T* (4)

reach 1.041 × 1011. Based on **Table 1**, we infer that natural convection could occur and the fluid state of air is a turbulent flow. The relation between *Nu* to the temperature difference can be obtained when relevant parameters are substituted into

*Nu* = 11000 (0.0740Δ*T*)1/3 (5)

In other cases, Eq. (1) may be not satisfied; corresponding experimental physical

Although we took Ningwu ice cave as an example to show the equivalent thermal

There is an annually periodic change on the air temperature outside of Ningwu ice cave. Therefore, the air temperature inside the ice cave would show a periodic fluctuation corresponding to the heat conduction and convective processes. The evolution of the air temperature at the bottom of the ice cave is shown in **Figure 6**. **Figure 6a** represents the air temperature evolution in Ningwu ice cave during its initial 16 years of the formation process of ice deposit. It can be seen that the air temperature in the cave rises in three seasons (spring, summer, and autumn) and drops in winter only, showing annually periodic variation. The air temperature of Ningwu ice cave drops promptly in winter while rises slowly in other three seasons (spring, summer, and autumn), because the efficiency of heat transfer in these three seasons is much lower than that in winter. If the water-ice phase change process was considered (black line), the dropped rate of air temperature in summer is smaller than that without the water-ice phase change (red line), because the latent heat of the water-ice phase change is required to melt ice near the cave entrance, thus delaying the conduction of heat energy to the bottom of the cave, while the convective cooling in winter is so efficient that the temperature difference is minimized. The temperature decreases below 0°C all year-round after winter cooling for about 5 years. That means the permanent cave ice deposit can be formed

**Figure 6b** represents the annual cave air temperature periodic fluctuations for the case where the heat transfer process has lasted two centuries, long enough to have evolved to a quasi-stable cyclic state. The amplitude of the cave air temperature fluctuation is about 1.0°C (from −3.9 to −2.9°C). Ningwu ice cave has been

conductivity method, the method could be used directly to Zibaishan ice cave.

°C, *Gr* can

According to Eq. (4), once the temperature difference is only 10<sup>−</sup><sup>3</sup>

relations can also be found in this literature [30, 32–33].

*4.2.3 Simulated results of ice deposit forming and melting*

Because Zibaishan ice cave is a typical upright circular tube.

and then maintained after winter cooling for about 5 years.

where *g* is the acceleration of gravity, *β* is the coefficient of cubical expansion, *ΔT* is a temperature difference, *l* is a characteristic length, and υ is the coefficient of kinematic viscosity. The values are *g* = 9.8 m s<sup>−</sup><sup>2</sup> , *β* = 3.67 × 10<sup>−</sup><sup>3</sup> k<sup>−</sup><sup>1</sup> , *l* = 80 m, and *υ* = 13.30 × 10<sup>−</sup><sup>6</sup> m2 s<sup>−</sup><sup>1</sup> and are substituted into Eq. (3) to obtain


#### **Table 1.**

*Earth Crust*

to compute the ice-water phase change.

*4.2.2 Equivalent thermal conductivity*

bottom of the cave, and also to the geometry of the cave.

values of which are shown in **Table 1**. *Gr* is the Grashof number:

> m2 s<sup>−</sup><sup>1</sup>

*Gr* = *g*β*ΔT l*

of kinematic viscosity. The values are *g* = 9.8 m s<sup>−</sup><sup>2</sup>

phase change can reduce the rate of temperature change. So the phase change should be taken into account. The conducting process is governed by conduction equation, which is relatively easy to be computed, while for the convection process, the convection pattern of air and its thermal consequences are hard to determine exactly, because of the complicated geometry of the ice cave and complex varying boundary conditions. In consideration of this, a broadly used, simplified approach is applied in this study: evaluation of Nusselt number (*Nu*) and solving the conductive equation by introducing an equivalent thermal conductivity of the convecting air. For an upright circular tube, the physical relation between the temperature difference of the top and the bottom and *Nu* can be evaluated by adopting the experimental relation of natural convection. The enthalpy method can be adapted

In every computing time step of our modeling process, it is judged whether air convection occurs according to the temperature difference between the top and the bottom of the cave. If there is no convection, the simple conduction problem will be solved, while if convection occurs, an effective conductivity is used in the conduction equation.

The heat conduction equation and enthalpy method for dealing with water-ice phase change were detailed in our previous work [27]. Equivalent thermal conduc-

*Nu* is defined as the ratio of convection heat transfer to pure conduction heat transfer under the same conditions. The efficiency of energy transfer is *Nu* times greater than the conductive efficiency under the same conditions. Therefore, an equivalent thermal conductivity can be introduced, which is *Nu* times greater than the air thermal conductivity [28]. *Nu* is related to physical properties (e.g., viscosity and conductivity of air), to the temperature difference of air at the top and the

Ningwu ice cave can be evaluated by an upright circular tube. For an upright circular tube, *Nu* can be calculated based on experimental fluid thermodynamics studies. Once Eq. (1) is satisfied [29, 30] (the case for Ningwu ice cave), the natural

*d*/*h* ≥ 35/ *Gr*1/4 (1)

*Num* = *C* (*Gr* ⋅ *Pr*)*<sup>m</sup>*

where *d* and *h* in Eqs. (1) and (2) are the diameter and the height of a circular tube respectively; *Num* is the Nusselt number; *Pr* is the Prandtl number and is dependent only on the material, e.g., *Pr* is 0.7 of air; and *C* and *n* are constants, the

where *g* is the acceleration of gravity, *β* is the coefficient of cubical expansion, *ΔT* is a temperature difference, *l* is a characteristic length, and υ is the coefficient

and are substituted into Eq. (3) to obtain

3

*<sup>n</sup>* (2)

/ υ<sup>2</sup> (3)

k<sup>−</sup><sup>1</sup>

, *l* = 80 m, and

, *β* = 3.67 × 10<sup>−</sup><sup>3</sup>

tivity method was a key point of computing and was shown as follows:

convection heat transfer relation [29, 31] can be expressed as Eq. (2):

**92**

*υ* = 13.30 × 10<sup>−</sup><sup>6</sup>

Gr *number and constant for different flow types [28].*

$$Gr = \mathbf{1.041} \times \mathbf{10}^{14} \,\Delta T \tag{4}$$

According to Eq. (4), once the temperature difference is only 10<sup>−</sup><sup>3</sup> °C, *Gr* can reach 1.041 × 1011. Based on **Table 1**, we infer that natural convection could occur and the fluid state of air is a turbulent flow. The relation between *Nu* to the temperature difference can be obtained when relevant parameters are substituted into Eq. (2):

$$\text{Nu} = \text{11000} \text{ (0.0740 \Delta T)}^{1/3} \tag{5}$$

In other cases, Eq. (1) may be not satisfied; corresponding experimental physical relations can also be found in this literature [30, 32–33].

Although we took Ningwu ice cave as an example to show the equivalent thermal conductivity method, the method could be used directly to Zibaishan ice cave. Because Zibaishan ice cave is a typical upright circular tube.

#### *4.2.3 Simulated results of ice deposit forming and melting*

There is an annually periodic change on the air temperature outside of Ningwu ice cave. Therefore, the air temperature inside the ice cave would show a periodic fluctuation corresponding to the heat conduction and convective processes. The evolution of the air temperature at the bottom of the ice cave is shown in **Figure 6**.

**Figure 6a** represents the air temperature evolution in Ningwu ice cave during its initial 16 years of the formation process of ice deposit. It can be seen that the air temperature in the cave rises in three seasons (spring, summer, and autumn) and drops in winter only, showing annually periodic variation. The air temperature of Ningwu ice cave drops promptly in winter while rises slowly in other three seasons (spring, summer, and autumn), because the efficiency of heat transfer in these three seasons is much lower than that in winter. If the water-ice phase change process was considered (black line), the dropped rate of air temperature in summer is smaller than that without the water-ice phase change (red line), because the latent heat of the water-ice phase change is required to melt ice near the cave entrance, thus delaying the conduction of heat energy to the bottom of the cave, while the convective cooling in winter is so efficient that the temperature difference is minimized. The temperature decreases below 0°C all year-round after winter cooling for about 5 years. That means the permanent cave ice deposit can be formed and then maintained after winter cooling for about 5 years.

**Figure 6b** represents the annual cave air temperature periodic fluctuations for the case where the heat transfer process has lasted two centuries, long enough to have evolved to a quasi-stable cyclic state. The amplitude of the cave air temperature fluctuation is about 1.0°C (from −3.9 to −2.9°C). Ningwu ice cave has been

**Figure 6.**

*Ice deposit forming process of Ningwu ice cave. (a) Initial stage of the formation process and (b) quasi-stable stage.*

opened to tourists, and about 1500 tourists visited the cave per day from May to October. Consequently, the cave air temperature has been disturbed. The lowest air temperature inside the cave was −1.5°C by our in situ measurement on 5 June 2012. Through the record in the literature, the actual measured inside air of the ice cave ranges between −1.0 [20], −4.0, and −6.0°C [19]. The inconsistent measurements may be contributed to variable measuring methods, measuring time and measuring positions. Similar to what is illustrated in **Figure 6a**, the cave air temperature shows annually periodic fluctuation. The overall rising rate of cave air temperature is lower than its dropping rate, resulting from variable heat transfer efficiency of conduction and convection. The variation on cave air temperature of the model with the water-ice phase change considered (black line) is almost the same as that without the water-ice phase change considered (red line). The reason is that, although the water-ice phase change was considered at the computation, the ice deposit temperature could be always maintained below 0°C when it reaches a quasistable cyclic state, and thus no water-ice phase change actually occurred.

The ice deposit in Ningwu ice cave would be melted if there is no air heat convection in winter. Taken the temperature distribution at 220 years (**Figure 6b**) as initial temperature, the evolution of temperature distribution can be computed with or without considering the water-ice phase change effects. The relatively modeling results are presented in **Figure 7** as a black line and a red line, respectively. They are the same when the air temperature is below the water-ice phase change temperature. However, the ice deposit body takes much longer to be thawed when the water-ice latent heat of melting is taken into account than when it is not considered. It takes 37 years to melt all of ice deposit with the water-ice phase change, while Ningwu ice deposit may be melted completely within 23 years without the water-ice phase change.

#### **4.3 Dynamic balance of water-ice and energy**

Water and ice of Ningwu ice cave may be in a dynamic equilibrium state. Ice stalactites and ice stalagmites (**Figure 1a**) can be found everywhere inside of Ningwu ice cave. The observation indicates that water infiltrates into Ningwu ice cave through fractures of limestone throughout the year and then forms ice deposit. Meanwhile, ice deposit at the bottom of Ningwu ice cave could be melted under geothermal flow, and then the produced water infiltrates into places beneath Ningwu ice cave. There are no direct observations to support this inference. This hypothesis may be correct, because of the conservation of water.

**95**

again.

**Figure 7.**

**5. Discussion**

*A Review of Chinese Ice Caves*

*DOI: http://dx.doi.org/10.5772/intechopen.89178*

predominantly transferred by conduction; if *Ra* exceeds the critical value, heat energy is mainly transferred by the air natural convection. We proposed that the evolution of ice deposit in Ningwu ice cave is a typical self-regulating process. When too much ice deposit accumulates in the cave, then the cavity will become small gradually. Thus, *Ra* and *Nu* will be reduced, meaning the freezing efficiency decreases, resulting in some of the cave ice deposit melting. Ice deposit melting will lead the cavity to become large, and thus the corresponding *Ra* and *Nu* will be risen. This implies that the freezing efficiency will increase and ice deposit in Ningwu ice cave will grow

*Illustration of internal temperature evolution of ice cave while ice deposit is melting.*

We proposed that Zibaishan ice cave may be also a self-regulating system, controlled by air natural convection between inside and outside of the ice cave.

In spring, summer, and autumn, heat energy was transferred into Ningwu ice cave through air and wall rock by thermal conduction, resulting in increasing the air temperature of the ice cave limitedly, because the efficiency of conduction is relatively low. But in winter, the air temperature of the ice cave decreases promptly, caused by air natural thermal convection. The water-ice phase change buffers the energy exchange. Considering these physical mechanisms, the modeling results present that (1) starting from a normal ground temperature distribution, a year-round ice deposit will be produced in the cave within a decade, about 5 years (**Figure 6a**), and the air temperature of the cave will drop gradually for more than a century, and also that (2) the air temperature of the ice cave will finally reach a quasi-stable cyclic state and will vary within a certain range (less than 1.0°C, from −3.9 to −2.9°C). At this stage, the annual total heat energy conducted into the cave and the heat energy removed from the cave by air natural convection are balanced. Installing an airtight door at a cave entrance is what the Wudalianchi National Geological Park has done to "prevent" the ice deposit from melting every night during the tourist season and the entire winter. When the cave airtight door is closed, in fact, the air convection in winter is blocked. Consequently, cold air cannot enter the cave and cannot remove heat energy from the cave. Accumulation of heat energy conducted by air and rock will eventually cause the ice deposit to melt in Wudalianchi cave. Our modeling results present that it takes less than 40 years to completely melt the whole ice deposit in the cave. This means that Ningwu ice cave is probably not undergoing melting of the ice deposit at present. This study

The Rayleigh number (*Ra*, a dimensionless number) is related to the buoyancy-driven flow. When *Ra* is below a critical value of that fluid, heat energy is *A Review of Chinese Ice Caves DOI: http://dx.doi.org/10.5772/intechopen.89178*

*Earth Crust*

**Figure 6.**

*stage.*

opened to tourists, and about 1500 tourists visited the cave per day from May to October. Consequently, the cave air temperature has been disturbed. The lowest air temperature inside the cave was −1.5°C by our in situ measurement on 5 June 2012. Through the record in the literature, the actual measured inside air of the ice cave ranges between −1.0 [20], −4.0, and −6.0°C [19]. The inconsistent measurements may be contributed to variable measuring methods, measuring time and measuring positions. Similar to what is illustrated in **Figure 6a**, the cave air temperature shows annually periodic fluctuation. The overall rising rate of cave air temperature is lower than its dropping rate, resulting from variable heat transfer efficiency of conduction and convection. The variation on cave air temperature of the model with the water-ice phase change considered (black line) is almost the same as that without the water-ice phase change considered (red line). The reason is that, although the water-ice phase change was considered at the computation, the ice deposit temperature could be always maintained below 0°C when it reaches a quasi-

*Ice deposit forming process of Ningwu ice cave. (a) Initial stage of the formation process and (b) quasi-stable* 

stable cyclic state, and thus no water-ice phase change actually occurred. The ice deposit in Ningwu ice cave would be melted if there is no air heat convection in winter. Taken the temperature distribution at 220 years (**Figure 6b**) as initial temperature, the evolution of temperature distribution can be computed with or without considering the water-ice phase change effects. The relatively modeling results are presented in **Figure 7** as a black line and a red line, respectively. They are the same when the air temperature is below the water-ice phase change temperature. However, the ice deposit body takes much longer to be thawed when the water-ice latent heat of melting is taken into account than when it is not considered. It takes 37 years to melt all of ice deposit with the water-ice phase change, while Ningwu ice deposit may be melted completely within 23 years without the

Water and ice of Ningwu ice cave may be in a dynamic equilibrium state. Ice stalactites and ice stalagmites (**Figure 1a**) can be found everywhere inside of Ningwu ice cave. The observation indicates that water infiltrates into Ningwu ice cave through fractures of limestone throughout the year and then forms ice deposit. Meanwhile, ice deposit at the bottom of Ningwu ice cave could be melted under geothermal flow, and then the produced water infiltrates into places beneath Ningwu ice cave. There are no direct observations to support this inference. This

The Rayleigh number (*Ra*, a dimensionless number) is related to the buoyancy-driven flow. When *Ra* is below a critical value of that fluid, heat energy is

hypothesis may be correct, because of the conservation of water.

**94**

water-ice phase change.

**4.3 Dynamic balance of water-ice and energy**

**Figure 7.** *Illustration of internal temperature evolution of ice cave while ice deposit is melting.*

predominantly transferred by conduction; if *Ra* exceeds the critical value, heat energy is mainly transferred by the air natural convection. We proposed that the evolution of ice deposit in Ningwu ice cave is a typical self-regulating process. When too much ice deposit accumulates in the cave, then the cavity will become small gradually. Thus, *Ra* and *Nu* will be reduced, meaning the freezing efficiency decreases, resulting in some of the cave ice deposit melting. Ice deposit melting will lead the cavity to become large, and thus the corresponding *Ra* and *Nu* will be risen. This implies that the freezing efficiency will increase and ice deposit in Ningwu ice cave will grow again.

We proposed that Zibaishan ice cave may be also a self-regulating system, controlled by air natural convection between inside and outside of the ice cave.
