Gr number and constant for different flow types .
Ice caves are a rare geological phenomenon. Ningwu ice cave, Zibaishan ice cave, and Wudalianchi ice cave are the most famous ice caves in China. We described each one in detail and carried out thermal-elastic modeling and heat conduction modeling to investigate the stability of the cave and the formation of ice deposit. In order to quantitatively study the mechanism of formation and preservation of the ice cave, we applied the FEM to model the heat exchange in the ice cave. The modeling results revealed that there is the seasonal asymmetric energy exchange. Heat energy is conducted inefficiently into the ice cave from outside and wall rock in spring, summer, and autumn. While in winter, heat energy is transferred very efficiently due to the air natural convection, thus cooling it down. We proposed that Ningwu ice cave and Zibaishan ice cave may be a self-regulating system, respectively. At Wudalianchi ice cave, airtight doors have been installed at these ice caves’ entrances. This actually prevents cooling in winter. We expect that no airtight door will be fixed at each ice cave’s entrance, and few people enter the ice cave before comprehensive and detailed studies, avoiding further affecting its natural conditions.
- ice cave
- air convection
- energy exchange
- water-ice phase change
- numerical modeling
Many caves developed in karstic regions and basaltic regions. But ice caves, permanent ice deposit preserved within caves, are a rare geological phenomenon. An ice cave is a type of natural cave that contains significant amounts of perennial ice. The most famous of ice caves are Eisriesenwelt ice cave, Austria [1, 2, 3]; Dobšinská ice cave, Slovakia [4, 5]; Scarisoara ice cave, Romania [6, 7]; and Monlesi ice cave, Switzerland [8, 9]. The Eisriesenwelt ice cave is the largest ice cave in the world. More than 10 ice caves have been found in China, and Ningwu, Zibaishan, and Wudalianchi ice caves are the most famous ice caves.
As long ago as 1861, studies of ice caves began . In the background of global climate change, seven international conferences on ice caves have been held in recent decades . Several articles have documented seasonal air temperature fluctuations of several degrees from cave systems [12, 13, 14]. Consequently, in order to approximate the impact of climatic conditions on cave environments, a better interpretation of subsurface heat transfers is necessary . In addition to this, ice caves are tourism resources. A better interpretation of subsurface heat transfers may assist in managing ice caves more scientifically.
Empirical calibrations have been carried out previously to evaluate the spatiotemporal distribution of cave temperature as a mathematical relation of the outside atmospheric conditions [15, 16]. In temperate karstic natural conditions, interpretation of the existence of subsurface ice deposits represents probably the most severe test for models of the magnitude and direction of heat and mass transfers induced by cave air circulation . In theory and application (mathematics and engineering), the finite element method (FEM) and the finite difference method (FDM) are classical methods for finding approximate solutions of partial differential equations, which governed physical processes (including heat transfer). The details of FEM and FDM can be found in many textbooks [17, 18]. These methods would be useful for ice cave studies.
Since Ningwu ice cave was found, ice cave studies began in 1998 in China. Ningwu ice cave in Shanxi Province has been broadly reported during the past decade [19, 20], but little was known about the physical processes controlling the formation and preservation of permanent subsurface ice deposits under temperate climate conditions . In addition to, basaltic cave stability is significant for ice preservation. FEM was applied to investigate the energy exchange of Ningwu ice cave and then quantitatively interpret the formation and preservation mechanism of the ice deposit. Thermal-elastic modeling was conducted to study lava tube stability during its cooling process.
2. Geological setting and characteristics of three ice caves in China
2.1 Ningwu ice cave
Ningwu ice cave (38°57′N, 112°10′E; 2121 m above sea level (asl)) is called “ten thousand years ice cave” by local people. It is located in the shaded slope of Guancen Mountain, Ningwu County, and Shanxi Province. The host rock consists of Ordovician Majiagou limestone, dolomitic limestone, argillaceous dolomite, and thin brecciated limestone and is densely fractured .
As a part of Ningwu National Geological Park, Ningwu ice cave is an important tourist attraction. Above 1500 visitors enter the cave per day from May to October. It has only a single entrance. People can walk into the inside of the cave by wooden spiral stairs. Ice covers the host rock almost completely. Ice stalactites and ice stalagmites can be found in all parts of the cave (Figure 1a).
In order to investigate the fine geometry of Ningwu ice cave, a classic geophysical exploration (using magnetotelluric measurement) was carried out and produced a high-resolution two-dimensional vertical cross-section of the ice cave  (Figure 1b). Figure 1b illustrates a bowling pin-shaped room, and the cave space is about 85 m depth. The widest part is in the middle of the cave with a width of 20 m.
The outside of Ningwu ice cave keeps a temperate climate. From June to September, the mean air temperature is about 14.6°C, and the average annual air temperature is 2.3°C . The nearest meteorological station to Ningwu ice cave is Wuzhai station (about 320 m lower than Ningwu ice cave), at which the daily temperature was measured continuously from 1957 to 2008. We averaged the observational air temperature at Wuzhai meteorological station to get the annual temperature and then derived the mean annual temperature at Wuzhai station. The difference between the average annual air temperature at Ningwu ice cave and that at Wuzhai station can be calculated. Finally, we obtained the annual temperature variation outside Ningwu ice cave through reducing the annual temperature at Wuzhai station by the difference.
2.2 Zibaishan ice cave
Zibaishan ice cave (33.7°N, 106.68°E; 2450 m asl) was found in August 2015, and it is located in the shaded slope of Zibai Mountain, Liuba County, Shanxi Province. Fieldwork was conducted in May 2016, and then the geometry and temperature profile of the ice cave was measured. The measurements are shown in Figure 2a. The cave space has a depth of about 90 m. The widest part is at the bottom with a width above 25 m.
Zibaishan ice cave, located in Zibaishan National Forest Park, is also a major tourist attraction. Most of the tourists only visited the opening of the ice cave, and a few experienced adventurers with professional climbing equipment could reach the interior of the ice cave. Consequently, Zibaishan ice cave maintained better natural environments (because of little artificial effects) than Ningwu ice cave. Ice deposits of the ice cave mainly included two parts: (1) ice cone which has a diameter of 25 m and a height of 16 m (Figure 2b) and (2) ice stalactites (Figure 2c). Ice deposits covered wall rocks of the ice cave partially. The host rock consists of limestone.
The outside of Zibaishan ice cave belongs to the warm temperate zone and humid monsoon section. The external annual mean air temperature is 11.5°C, and the frost-free period is 214 days per year . Generally, ice deposits could not be preserved under such a temperate environment. Therefore, there must be some special cooling mechanisms of Zibaishan ice cave.
2.3 Wudalianchi ice cave
There are two ice caves at Wudalianchi National Geological Park (48.647°, 126.25°, 400 m asl), Heilongjiang Province: one is Bailong (means white dragon; also named Dixiabinghe) ice cave, and another is Shuijinggong (means crystal palace) ice cave.
We conducted simple field work on Bailong ice cave in August 2012. The host rock is basalt. Dense fractures developed near the entrance. Basaltic pillars can be seen within the cave. The farthest distance we can reach is about 270 m from the entrance of the ice cave. Wang and Zhu measured the geometry and basaltic fracture of Bailong ice cave by applying geological radar (the antenna frequency is 40 MHz)  (Figure 3). From the radar profile (Figure 3b), we can see that seven large fractures with the same dip direction were developed above the ice cave.
Bailong ice cave was not a pure natural ice cave, but ice blocks were moved into the cave to decorate it, and thus its natural conditions were destroyed. The temperature in the ice cave was about −2°C in August 2012. Shuijinggong ice cave is much smaller than Bailong ice cave but with more man-made ice block which was used to decorate it. Therefore, it is difficult to obtain any natural information from Shuijinggong ice cave.
3. Formation of cave
Caves must be developed before the formation of their internal ice. Cave geometry mainly depended on the lithology of their wall rocks. Two types of lithology were involved in these three ice caves: limestone and basalt.
3.1 Limestone as the host rock
The host rock of Ningwu ice cave and Zibaishan ice cave mainly consists of limestone. Limestone is composed of calcium carbonate which is an unsteady material and could be readily dissolved. Limestone areas could become karst landform under wet and warm climate in the geological period. The limestone cave developed predominantly in the perpendicular direction, for example, Ningwu ice cave and Zibaishan ice cave. Generally, the karst region is not suitable for ice preservation. Ice preservation in karst cave may indicate the climate change from warm to cold globally or the presence of incomprehensible cooling mechanisms locally.
3.2 Basalt as the host rock
Basaltic caves (often as tubes) were formed in the process of the viscous lava flow, for example, Wudalianchi ice cave. Mechanical stability of lava tubes depends on their size (diameter), buried depth, and geometry in the cooling period. We conducted a series of two-dimensional thermo-elastic finite element calculations to estimate the effects of these three controlling factors .
Figure 4 shows the principle stress distributions around a lava tube of diameter 10 m at the last computing time step. It can be seen from Figure 4a that the largest maximum principal stress is at the top and bottom sides of the lava tube, and it is tensile stress. The left and right sides are also tensile stresses, but the magnitude is much smaller. Figure 4b is the minimum principal stress distribution. Obviously, the minimum principal stress of the upper and lower sides of the lava tube is tensile stress, while that of the left and right sides is compressive stress, and the magnitude of the stress around the lava tube is equal. These modeling results reveal that the most vulnerable place of the lava tube is on the top and bottom, and the rupture type is tension rupture.
We investigated three controlling factors (diameter, depth, and geometry) of lava tubes in the cooling period by FEM method, and the modeling results were illustrated in Figure 5. At the last time step, maximum principal stress of Figure 5a (normal lava tube) is smaller than that of Figure 5b (bigger lava tube), implying bigger lava tubes were easier ruptured. Similarly, lave tubes with larger buried depth would be easier ruptured than normal ones (Figure 5c). Elliptical lava tubes would be more stable than normal lava tubes (Figure 5d).
4. Formation of ice deposit
After caves were formed, ice may be generated gradually in a suitable climate. For basaltic cave (or lava tube), there may be a clear time boundary between cave formation and ice formation, because lava tubes are steady in water circulation. For limestone cave, there is no such clear time boundary due to limestone instability in wet environments during the geological period.
4.1 Qualitative analysis of cave cooling by air convection
There are variable hypotheses about the preservation mechanism of ice body in Ningwu ice cave. Chen  proposed that there is a “cold source” below Ningwu ice cave, which generates the negative geothermal anomaly and then preserves the ice body. Meng et al.  ascribed the existence of ice deposit to the combined effect of multiple factors, i.e., geographical location, the “icehouse effect,” the “chimney effect,” and the “thermal effect” produced by the ice deposit and the “millennial volcano.” Unfortunately, they did not supply more details on these factors. Gao et al.  considered two aspects, terrain and climate, and proposed that far more cold air than warm air entered the region, and thus the ice cave stayed cold over a year. The subsurface temperature usually rises with depth at a geothermal gradient of 1.0–3.0°C (100 m)−1 , which cannot support a permanent “cold source” underground. Furtherly, even if a cold region had somehow formed, it would be warmed up by the geothermal flux through the geological period, because the host rock successively transfer heat energy to the “cold source.” At another point, it is a temperate climate outside of the cave. It is difficult to maintain ice deposit in a warm climate without an efficient cooling mechanism. Therefore, we proposed that there must be a sustainable and efficient mechanism to remove the heat from underneath and ensure the existence of the ice deposit.
There is an annual cyclic variation on the air temperature outside Ningwu ice cave: it is warmer than the inside temperature in spring, summer, and autumn, but colder in winter. Because Ningwu ice cave has only a single opening at the top of the cave, cold air of the ice cave could be relatively heavy in spring, summer, and autumn and sinks into the bottom of the cave. Thus, it will not generate air natural thermal convection. In these three seasons, heat energy is transferred by conduction from outside down to the ice cave and from wall rock because of the terrestrial heat flux. Thermal conductivities are low for either wall rock (limestone) or air, and thus the conductive heat transfer efficiency is low. Consequently, the energy conducted to the inside of ice cave in these three seasons is quite limited. However, in winter, it is colder than the inside of the cave, and thus the air outside of ice cave could be heavier than the inside. Gravitational instability is generated, and air thermal convection could occur. The outside cold air promptly flows into the cave to cool it down, and it removes the heat energy out of the cave, which is conducted into the cave from the host rock and through the opening in spring, summer, and autumn. Convective heat transfer is much more efficient than conduction; therefore, the heat energy convected out of the cave in winter is enough to balance the heat conducted into the cave year-round. The annual heat budget of income and output is balanced, so the cave would be in a cyclic state with very small temperature fluctuations, and the average temperature is always lower than 0°C; thus ice deposits in the ice cave can be kept.
Studies in Zibaishan and Wudalianchi ice caves were still at a relatively low level by far. We proposed that similar air convection cooling could occur in these two ice caves.
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 .
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 . 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 . 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):
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:
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−2, β = 3.67 × 10−3 k−1, l = 80 m, and υ = 13.30 × 10−6 m2 s−1 and are substituted into Eq. (3) to obtain
According to Eq. (4), once the temperature difference is only 10−3°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):
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 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 , −4.0, and −6.0°C . 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-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-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.
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 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.
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 also suggests that scientific management is significant for sustainable application of natural resources. Otherwise, well-meaning acts such as installing a door to completely seal the entrance will actually destroy the ice cave in a few decades.
We expect that no airtight door will be fixed at Zibaishan ice cave entrance and few people enter the ice cave before comprehensive and detailed studies, avoiding affecting its natural conditions.
This chapter has described three famous ice caves in China, i.e., Ningwu ice cave, Zibaishan ice cave, and Wudalianchi ice cave. We reviewed quantitative analysis of the formation and preservation mechanism of an ice deposit in Ningwu ice cave, a quasi-static ice cave. The systematic FEM computing leads to the conclusion below: the air natural convection in winter is the crucial controlling factor for forming and sustaining the ice deposit in the cave and can cool the ice cave efficiently. Heat transfer by conduction in spring, summer, and autumn is very limited to warm up the cave. Water-ice phase change plays an important role in summer and can prevent melting of ice. We proposed that Ningwu ice cave and Zibaishan ice cave may be a self-regulating system, respectively, controlled by air natural convection between inside and outside of the ice caves. We expect that no airtight door will be fixed at each ice cave entrance and few people enter the ice cave before comprehensive and detailed studies, avoiding further affecting its natural conditions.
We appreciate Prof. Lada Bozic’s invitation. Many thanks to Prof. Huai Zhang, Yuanze Zhou, and Guiping Zhao for their kind help in field work. We thank Jiaqi Sun who supplied the pictures of Zibaishan ice cave. This research is supported by the National Natural Science Foundation of China (NSFC), project 41590860 and U1839207.