3. Results and discussion

Two dimensionless numbers, the Stefan number ( ${\textit{Ste}}$ ) and the Fourier number ( ${\textit{Fo}}$ ), are very important to characterise the thawing process (Naldi et al. Reference Naldi, Martino, Dongellini and Lorente2024). The non-monotonic dependence of core thawing time on pool temperature has been achieved in figure 2(c). Results for different pool temperatures ( ${\textit{Ste}}_{\textit{pool}}$ ) with a real density variation of water compared with an artificially linear density variation are presented in the figure. More details about the dimensionless numbers are provided in the supplementary material.

Figure 2. Dimensionless thawing time changing with pool temperature. (a) Average ice fraction changes with the dimensionless time ${\textit{Fo}}$ in three thawing stages which are distinguished by the red dash lines. The 45 $^\circ$ angle line identifies the bounds between stage 2 and 3. (b) Dominant zones of two different flow patterns in the pool: boundary zones and zigzag bulk zone away from the core and walls. (c) Full thawing dimensionless time ${\textit{Fo}}$ of inner core versus the dimensionless pool temperature ${\textit{Ste}}$ . The red solid circles are results with a real $\rho {-}T$ of water and the grey squares are calculated with an ‘artificial’ linear $\rho {-}T$ . (d) Temperature ( $T$ ) and velocity ( $U$ ) profiles combined with the fluid fraction profiles for three typical cases (dashed squares in panel (c): (i) ${\textit{Ste}}_{\textit{pool}}$ = 0.44; (ii) ${\textit{Ste}}_{\textit{pool}}$ = 0.56; (iii) ${\textit{Ste}}_{\textit{pool}}$ = 0.88. The black circle in panel (d) represents the shell of core.

To better describe the thawing process, three distinct thawing stages are recognised in figure 2(a). Stage 1 exhibits the rapid dynamics during initial thawing, forming an ice shell on the subcooled core ( $T_{\textit{core}}\lt T_m$ ) due to the steep radial temperature gradients. Although a net core thawing occurs, the ice-shell growth rate exceeds thawing, causing a slight increase in average ice fraction at the beginning. A critical pool temperature exists where the ice-shell formation is fully suppressed. Stage 2 subsequently develops with an accelerated thawing rate, characterised by a steep negative gradient. This enhancement results from an increased effective thermal conductivity within the ice–water mixed layer, sustained by maintained temperature gradients. The transition to stage 3 occurs at the demarcation point where the gradient equals -1 (45 $^\circ$ line). During stage 3, the thawing rates decline progressively as the growing water-layer thickness and diminishing ice coverage reduce both the thermal gradients and the conductivity, ultimately converging to classical Stefan problem behaviour (Hu & Argyropoulos Reference Hu and Argyropoulos1996). However, the thawing process changes with the pool temperatures. The pool-temperature variation is used to investigate the ‘thermal currents battle in ice-core thawing problem in a water pool’ to figure out the best thawing pool temperature and the potential mechanism.

An interesting non-monotonic thawing relationship has been observed, as shown by the red solid circles in figure 2(c), between the initial pool temperatures and the complete thawing time of the ice core. This implies that, interestingly, with fixed computational domain boundaries and ice-core parameters, a higher temperature of the water pool does not necessarily lead to a shorter thawing time, which is consistent with the old Chinese saying. The simulation results show that before the pool temperature reaches ${\textit{Ste}}_{\textit{pool}}$ = 0.50 ( ${\sim}10\,^\circ\mathrm{C}$ ), the complete thawing time decreases significantly and reaches a local minimum. Subsequently, as the pool temperature gradually increases, the thawing time rises and peaks at approximately ${\textit{Ste}}_{\textit{pool}}$ = 0.88 ( ${\sim}40\,^\circ\mathrm{C}$ ). The peak marks the end of the non-monotonic region. Beyond this point, the core thawing time decreases monotonically with rising pool temperatures, returning to an intuitive outcome.

In contrast, the simulation results with an ‘artificial’ linear $\rho {-}T$ relationship are plotted in the same figure, as shown as the grey squares in figure 2(c). The thawing time shows a monotonic dependence on the pool temperature as a result, which indicates that it is the density-temperature anomaly at ${\sim}4\,^\circ\mathrm{C}$ of water that triggers the non-monotonic thawing behaviour of ice in a water pool. Furthermore, the viscosity also follows a temperature-dependent variation consistent with the actual water. The NOB effect fundamentally drives the non-monotonic influence of pool temperature on complete core thawing time. Within this mechanism, the realistic density variation of water serves as the primary driver, while the temperature-dependent viscosity acts as a secondary factor that amplifies physical perturbations.

Visualisation of temperature and velocity fields reveals that there is a competition of two flow patterns noteworthy in the pool, as illustrated in figure 2(d). One is the chaotic flow in the bulk region away from the ice core and the boundary, as shown as around the zigzag line in figure 2(b), and the other appears as a ‘thermal jet’ leading to a large circular convective flow from the core bottom. Figure 2(d)(i)–(iii) shows such a flow pattern at different pool temperatures. The convective flow is near the boundaries of both core and walls of the pool, and the chaotic flow is in the bulk region between boundaries. The size of chaotic flow region will be compressed by the stronger convective flow around when the pool temperature increases (from case i to case iii). This means that a higher pool temperature may lead to a more stable flow field in the pool, which is consistent with the previous observation (Sugiyama et al. Reference Sugiyama, Calzavarini, Grossmann and Lohse2009). An interesting observation from figure 2(d) is that as ${\textit{Fo}}$ increases the flow in the pool seems evolving from chaotic to regular and then to re-chaotic flows in the non-monotonic region. The rich flow mechanism in this process remains undiscovered.

In the thawing process, as the inner structure of frozen core is considered as a porous medium, the flow inside the core is generally confined so that the heat transfer inside the core is only by conduction. Meanwhile, the heat transfer in the pool is majorly by convection since the radiation is negligible in this low-temperature scope. To better understand the transition mechanism of flow patterns during the thawing of ice core in a water pool, the instantaneous $Nu_b$ at the bottom surface of the ice core can be calculated by

(3.1) \begin{equation} Nu_{b}(t)\ =\ \frac {h(t)L}{k_{l}}, \end{equation}

where $h(t)$ is the instantaneous local convective heat transfer coefficient and expressed as $h(t) = q(t) / \Delta T(t)$ , $q(t)$ the instantaneous local heat flux, $\Delta T(t)$ the local temperature difference with reference boundary temperature, $k_l$ the thermal conductivity and $L$ the characteristic length which is the radius of ice core.

Figure 3 shows the evolution of dimensionless heat convection coefficient $Nu_b$ versus the dimensionless time ${\textit{Fo}}$ for a given pool temperature at ${\textit{Ste}} = 0.56$ . The simulation starts from a stationary state and needs some time for initialisation. Therefore, only the Nusselt number variation after the velocity-field initialisation ( ${\textit{Fo}} \gt 0.01$ ) is plotted in the figure. Combined with considerations of fluid field contours, as shown by the insets in figure 3, the flow and heat transfer during thawing experiences four clear regimes: temporary stable regime, oscillatory regime, transitional regime and fully chaotic regime. A similar stable and unstable convection also happens on the ice generation process (Wang et al. Reference Wang, Calzavarini, Sun and Toschi2021b ). The instantaneous Nusselt number together with the velocity field helps to quantify the bounds between regimes at ${\textit{Fo}}$ = 0.018, 0.04 and 0.064. First comes a temporary stable regime which has a clear ‘thermal jet’ with small oscillations. In this regime, the angle between ‘thermal jet’ and gravity is less than 10 $^\circ$ and $Nu$ slightly decreases with small fluctuations. After the bounds at ${\textit{Fo}}$ = 0.018 is the oscillatory regime, where $Nu$ experiences distinct periodic oscillations which correspond to a periodic swinging of the ‘thermal jet’ in the velocity field. In this regime, the Nusselt number is still generally below the temporary stable regime on the bounds ( ${\textit{Fo}}$ = 0.018, $Nu$ = 16.12) and the ‘thermal jet’ is still clear. As the ${\textit{Fo}}$ increases to ${\textit{Fo}} \approx \ 0.04$ , it experiences a transitional regime, $Nu$ begins to exhibit drastic oscillations and significantly exceeds the bounds $Nu$ of the stable regime which suggests a much stronger heat exchange. The ‘thermal jet’ in the velocity field is diffused and intermittently disappears, which suggests the flow exhibits a pronounced tendency towards chaos. Finally, after ${\textit{Fo}} \approx \ 0.064$ , the $Nu$ becomes stable and reaches the fully chaotic regime.

Figure 3. Variation of the instantaneous local $Nu_b$ with ${\textit{Fo}}$ at ${\textit{Ste}} = 0.56$ . The flow regime transitions from temporary stable, oscillatory, transitional flow and then to fully chaotic flow with ${\textit{Fo}}$ . The dashed lines are the bounds of regimes by comprehensively considering the variation of $Nu_b$ and the velocity field. Selected representative velocity fields are annotated above the curve to illustrate the evolving flow regime ( ${\textit{Fo}} = 0.016, 0.02, 0.054, 0.076$ ).

There is a clear competition between the chaotic flow and the circulation leads by ‘thermal jet’, which decides the flow regime in the pool region. Consequently, non-monotonic behaviour depends on the flow dynamics in the pool and can be expected to be influenced by the size of the pool. Simulations with pool sizes in figure 4(a) show different non-monotonic curves between the inner core thawing time and the pool temperature. The non-monotoinc trend decreases with the pool size to ‘narrower and shallower’. Generally, the full thawing needs shorter time in a larger pool after passing the minimal extreme. There is a migration of the buoyancy maximum position, which shifts from the top of the pool (near the thermal boundary) to the bottom of the cold core. For a pool temperature below 4 $^\circ$ C, wall heating induces density inversion, generating chaotic flows. Critical point analysis therefore initiates at 4 $^\circ$ C (where the thermal expansion coefficient $\beta$ becomes negative, rendering the $Ra$ to be negative and incompatible with stability criteria). As the pool temperature increases, the buoyancy maximum initially resides at the top of the pool, where boundary-imposed temperature (and density) differences peak. Upon further pool temperature elevation towards boundary conditions, boundary-driven density differences diminish relative to those generated near the density extreme (4 $^\circ$ C) within the core region. This shift drives downward migration of the buoyancy maximum to the lower core, forming a clear descending ‘thermal jet’.

Figure 4. Dimensional thawing time of ice core for different pool sizes. (a) Real time of full thawing of the inner core versus the pool temperatures for three different pool sizes with the same core size. The pool sizes are 9.60, 12.00 and 14.40 cm. The maximum buoyancy positions are marked with red circles for four typical cases respectively by insets: $T_{\textit{pool}}$ = (i) 5 $^\circ$ C, (ii)1 5 $^\circ$ C, (iii) 45 $^\circ$ C, (iv) 55 $^\circ$ C. (b) Velocity contours at $T_{\textit{pool}}$ = 15 $^\circ$ C for three different pool sizes at the starting (thawing ratio = 0.004) and final (thawing ratio = 0.9) stage of thawing. The grey part in the core is the remaining ice.

Figure 4(b) shows that the large pool size helps the maintain of ‘thermal jet’. During thawing, boundary-driven flow (analogous to Rayleigh–Bénard convection) exhibits chaotic dynamics, imparting disordered disturbances to the pool layer (evident in figure 4 b at select thawing fractions). Conversely, the thermal jet generating from the lower core promotes the formation of a stable, pool-scale circulatory flow. Thus, variations in pool temperature instigate competition between these two distinct flow regimes (chaotic boundary disturbances versus stable core circulation), resulting in diverse flow patterns. The transition zone corresponds to a state of dynamic equilibrium where the effects of both flow types are comparable, placing the system within a potential well. Outside the potential well, one flow regime predominates decisively. There are two main mechanisms that dominate the two patterns.

1. (i) The anomaly-triggered chaotic flow (ACF) mechanism. The chaotic flow is driven by the temperature difference ( $\Delta T$ ) between the boundaries and the cold core, analogous to multidirectional Rayleigh–Bénard convection (incorporating both horizontal and vertical thermal gradients). Chaotic flow state that is directly triggered and sustained by the density anomaly (maximum) of water. This multidirectional nature drastically reduces stability, implying no steady-state solution. Its intensity scales positively with $\Delta T$ , quantified by the dimensionless ${\textit{Ra}}_{\textit{ACF}}$ (anomaly chaotic $Ra$ ). Anomaly-triggered chaotic flow exerts persistent disturbances, driving the system towards chaos if without competing mechanisms.

2. (ii) The normal natural convection (NNC). Normal natural convection circulation originates from a 4–8 $^\circ$ C protective layer adjacent to the cold core wall (Léard et al. Reference Léard, Favier, Le Gal and Le Bars2020), characterised by: (i) self-restoring tendency due to the density extremum of water ( ${\sim}4\,^\circ\mathrm{C}$ ); (ii) intense disturbance dissipation under Prandtl number ( $Pr$ ) ${\gt } 1$ (viscous dominance); (iii) thin-layer structure (similar to thermal boundary layer), acting as a gravity-driven source via maximal density (Hu, Zhang & Xia Reference Hu, Zhang and Xia2024). It initiates downward flow that forms pool-scale circulation along walls, quantified by ${\textit{Ra}}_{\textit{NNC}}$ (normal natural convection $Ra$ ). The effect is strong at high $T_{\textit{pool}}$ .

The instability analysis may lead us to a threshold of flow pattern transition (Bergé & Dubois Reference Bergé and Dubois1984) since the non-monotonic thawing is caused by the competition of the two mechanisms dominated by different pool temperature. The analysis of competition between ACF and NNC may help to get the critical value by

(3.2) \begin{align} {\textit{Ra}}_{\textit{NNC}}&=\frac {(T_{\textit{pool}}-T_{m,{\textit{core}}})\beta _{w\textit{ater}}l_{\textit{eff}}^3\left |\boldsymbol{g}\right |}{\nu _l\alpha _l}, \end{align}

(3.3) \begin{align} {\textit{Ra}}_{\textit{ACF}}&=\frac {(T_{\textit{BC}}-T_{\textit{core}})\beta _{w\textit{ater}}l_p^3\left |\boldsymbol{g}\right |}{\nu _l\alpha _l}, \end{align}

where $T_{\textit{pool}}$ is the pool temperature, $T_{m,{\textit{core}}}$ the inner core thawing point, $T_{\textit{BC}}$ the temperature of the constant boundary condition, $\beta _{w\textit{ater}}$ the thermal expansion coefficient which may be a nonlinear function of temperature, $l_{\textit{eff}}$ the effective thermal boundary length of the NNC layer zone, $l_{p}$ the pool length, the $\boldsymbol{g}$ the gravity acceleration, $\nu _l$ the viscosity of the liquid phase and $\alpha _l$ the thermal diffusivity of the liquid phase. The thermal boundary length $\delta _{\textit{thermal}}$ and the effective length $l_{\textit{eff}}$ are

(3.4) \begin{align} &\,\,\delta _{\textit{thermal}} \sim x{{\textit{Gr}}_x}^{-\frac{1}{4}},\quad {\textit{Gr}}_x=\frac {{\textit{Ra}}_x}{\textit{Pr}}, \end{align}

(3.5) \begin{align}& l_{\textit{eff}}=\left (\frac {T_{\textit{pool}}-T_4}{T_{\textit{pool}}-T_{m,{\textit{core}}}}\right )l_p{\textit{Ra}}^{-\frac {1}{4}}{\textit{Pr}}^{\frac {1}{4}}. \end{align}

where $T_4$ is 4 $^\circ$ C for water.

Flow stability requires ${\textit{Ra}}_{\textit{NNC}}\geq \varLambda {\textit{Ra}}_{\textit{ACF}}$ , and $\varLambda$ is marked as stability coefficient. The stability results in

(3.6) \begin{equation} \varLambda \leq \frac {\left (T_{\textit{pool}}-T_{m,{\textit{core}}}\right )\left (\frac {T_{\textit{pool}}-T_4}{T_{\textit{pool}}-T_{m,{\textit{core}}}}\right )^3{\textit{Ra}}^{-\frac {3}{4}}{\textit{Pr}}^{\frac {3}{4}}}{T_{\textit{BC}}-T_{\textit{core}}}. \end{equation}

The stabilisation threshold value $\varLambda _{c}$ for (3.6) marks the dimensionless demarcation point for the domination of flow patterns. In figure 2, the critical point is around $T_{\textit{pool}}$ = 12.50 $^\circ$ C and then the corresponding threshold is determined as $\varLambda _{c}$ = $6.14\times {10}^{-7}$ .

The non-monotonicity of thawing originates from the density-temperature anomaly of water, which triggers two competing mechanisms in the pool zone: the anomaly-triggered chaotic flow (chaotic flow) and the normal natural convection (steady circulatory flow). Their comparable strength will induce a transitional region (manifested as a downward-curving potential well). The realistic $\nu {-}T$ relationship is not the key driving source yet plays a secondary but important and quantitative tuning role in the non-monotonic thawing phenomenon. As a dissipative term in the momentum equation, the fluid viscosity influences flow developments. If the viscosity was assumed constant, the viscous dissipation could be underestimated at a low temperature and overestimated at a high temperature. Consequently, the extreme points (the local minimum and maximum) on the non-monotonic curve may shift to compensate for this under- or over-estimated flow resistance. Accordingly, the width and depth of the heat potential well are also expected to change (wider and deeper) due to the altered balance between the competing flow mechanisms.

Although very interesting phenomenon and new mechanisms have been explored through numerical simulations in this work, the results are still constrained by some limitations. First, two-dimensional models are adopted for simplified simulations to capture the key mechanisms in a three-dimensional real system. This is a general simplification. The flow streamlines and vortices must be more complex in a real three-dimensional system and the critical values, including the optimal thawing time and stability threshold etc., may vary quantitatively due to three-dimensional effects. However, a recent study reported that for a spherical ice-core thawing problem, the difference between 2-D and 3-D modelling results is almost negligible (Yang et al. Reference Yang, Van Den Ham, Verzicco, Lohse and Huisman2024b ), which means that the conclusions and mechanisms of the non-monotonic thawing of ice core in a water pool will not be adapted by the 3-D effects. Second, some minor flows are not considered in this work and are believed not to change the phenomena and mechanisms. For example, the convection inside the core induced by temperature difference and concentration difference of sugar is ignored in this study because the inner porous structure of the core makes the flow rather weak. Once the inner convection is considered, the heat transfer in the core would be enhanced and may lead to a more pronounced non-monotonic trend. The small-scale Stefan flow at the solid–liquid interface is also not considered in our simulation. To capture this near-interface flow, one should adopt multi-block or adaptive mesh refinement techniques near the solid–liquid interface to obtain full resolution of pressure and consider pressure-sensitive parameter influences. Third, we consider a constant wall temperature of the pool in our study. In reality, the constant-temperature boundary condition may be rigidly valid only when the pool is large enough. Otherwise, the surrounding wall temperature of the pool would dynamically decrease as the ice melts and the real convection boundary condition outside the wall of pool would be considered.