Melting of floating ice cylinders in fresh and saline environments

Abstract

Motivated by the need for a better understanding of the melting and stability of floating ice bodies, we experimentally investigated the melting of floating ice cylinders. Experiments were carried out in a tank, with ice cylinders with radii between 5 and 12 cm, floating horizontally with their axis perpendicular to gravity. The water in the tank was at room temperature, with salinities ranging from 0 to 35 g l⁻¹. These conditions correspond to Rayleigh numbers in the range 10$^5\lesssim$ Ra $\lesssim$ 10$^9$. The relative density and thus the floating behaviour was varied by employing ice made of H$_2$O–D$_2$O mixtures. In addition, we explored a two-layer stable stratification. We studied the morphological evolution of the cross-section of the cylinders and interpreted our observations in the context of their interaction with the convective flow. The cylinders only capsize in fresh water but not when the ambient is saline. This behaviour can be explained by the balance between the torques exerted by buoyancy and drag, which change as the cylinder melts and rotates. We modelled the oscillatory motion of the cylinders after a capsize as a damped nonlinear oscillator. The downward plume of the ice cylinders follows the expected scalings for a line-source plume. The plume’s Reynolds number scales with Rayleigh number in two regimes, namely Re $\propto$ Ra$^{1/2}$ for Ra $\lt \mathcal{O}(10^7)$ and Re $\propto$ Ra$^{1/3}$ for Ra $\gt \mathcal{O}(10^7)$, and the heat transfer (non-dimensional as Nusselt number) scales as Nu $\propto$ Ra$^{1/3}$. Although the addition of salt substantially alters the solutal, thermal and momentum boundary layers, these scaling relations hold irrespectively of the initial size or the water salinity. While important differences exist between our experiments and real icebergs, our results can qualitatively be connected to natural phenomena occurring in fjords and around isolated icebergs, especially with regard to the melting and capsizing behaviour in stratified waters.

1. Introduction

The lack of understanding of ice-sheet melting dynamics is the largest source of uncertainty in the projections of sea-level rise (Robel, Seroussi & Roe 2019). (Throughout the Introduction, the terms sea ice, ice-sheet, icebergs, etc. are used as in the reference they correspond to. The physics and interaction mechanisms of each are different, and the terms are not interchangeable.) The complexity of the phase-change physics and the uncertainties in the parametrisations (e.g. for albedo changes and for the interactions of the ice with ocean currents) implies that most of today’s climate models poorly model the sea ice variations (Stroeve et al. 2007; Malyarenko et al. 2020). The different contributors to sea level rise (among which, thermal expansion and ice sheets) are commonly assumed to be independent of each other, and this is deemed to be the reason for an underestimation of the uncertainty in sea-level projections (Le Bars 2018). One of the processes that connect the ocean with ice structures is calving at the Greenland Ice Sheet and Antarctic Ice Sheet, where ice chunks up to hundreds of metres of length (Orlowski 2012) detach from the edge of a glacier, becoming icebergs. On top of being a major component in both the Antarctic Ice Sheet and Greenland Ice Sheet mass balance (Bigg et al. 2014; Cenedese & Straneo 2023), icebergs and their freshwater fluxes have been proved to be a contributor to the freshwater balance in the Southern Ocean (Silva, Bigg & Nicholls 2006). An increase in the freshwater input from icebergs and meltwater discharge can have strong consequences on the global climate (Schloesser et al. 2019).

Oceanographic surveys are as a precaution kept far from icebergs, due to the unpredictability of capsizing events, thus limiting the data collection in close proximity of icebergs (Yankovsky & Yashayaev 2014). Correspondingly, limited research has been conducted on freely floating melting objects. A historically significant branch of studies developed in the second-half of the 20th century, when the interest in harvesting icebergs for freshwater motivated research on their melting during towing, see e.g. Hult & Ostrander (1973) and Russell-Head (1980). More recently, growing interest in iceberg research led to experimental and numerical investigations on the capsizing of icebergs. Burton et al. (2012) and Bonnet et al. (2020) model natural icebergs as plastic rectangular cuboids, with well-defined aspect ratios, but excluded melting. Their research focuses on the capsizing events: How it is affected by the aspect ratio, how energy is released by a calving event and how the capsizing triggers mixing and flow around the iceberg?

The essential physics of an isolated iceberg melting in the ocean is that of a Stefan problem (Rubinstein 1971), where a cold ice mass lays on top of a warm, saline water reservoir. Three fields control the melt rate of the ice: temperature, salinity and velocity of the fluid. First and foremost, the melting temperature of the water is a function of the salinity $T_{\textit{melt}}=T_{\textit{melt}}(S)$ and decreases with increasing salinity. Secondly, the Lewis number (Le $=\kappa _S/\kappa _T$ , the ratio between saline diffusivity $\kappa _S$ and thermal diffusivity $\kappa _T$ ) controls the relative thickness of the thermal and saline boundary layers. Lastly, the density of water is itself a function of temperature and salinity, $\rho _\infty =\rho _\infty (T,S)$ . The interactions between the melting and the water density are mutual, with the temperature and salinity fields that are affected by the melting ice, and the resulting density differences that set the water in motion, advecting water to the ice surface, whose temperature and salinity affect the melting itself.

Additional complexity arises if the free surface around the ice is considered, and the ice is unconstrained. In fact, in this configuration, the ice floats expose only a fraction of their volumes to air, and the floating behaviour is controlled by the ratio between ice and water densities, $\varLambda =\rho _{\textit{ice}}/\rho _\infty$ . The lower heat transfer rate with air compared with water causes a significant difference in the melt rate between the submerged ice and the ice exposed to air. Furthermore, the ice’s motion can be both in response to its own morphological evolution and to the liquid’s motion, with the ice exposing different parts of itself to the liquid, depending on its orientation.

Hitherto, most experimental laboratory-scale research has focused on the melting of submerged ice. The studies on boundary layers of submerged slabs are extensive (Josberger & Martin 1981; Carey & Gebhart 1982), with recent particular interest in the appearance of scallops, both in melting and dissolution problems (Davies Wykes et al. 2018; Cohen et al. 2020). At the same time, melting in symmetrical domains has received interest due to the relative simplicity of theoretical description of the heat transfer. For this reason, Fukusako et al. (1992) and Yamada et al. (1997) have investigated the melting of a fully submerged fixed horizontal ice cylinder, for different ambient temperature and salinity conditions. They grew ice around a cooled pipe in a water bath and then let it melt, imaging the cross-section with a single lens reflex (SLR) camera and visualising the flow using streakline photography. They concluded that the flow profile in the vicinity of the cylinder was laminar and bidirectional in the lower portion of the cylinder, and turbulent unidirectional in the upper portion of the cylinder, with the flow enhancing the melt rate.

Another line of research that has run in parallel to the aforementioned ones, is the study of plumes around a heated submerged object. The most relevant references here are the work of Grafsrønningen et al. (2011) and Grafsrønningen & Jensen (2012, 2017), who extensively investigated the plumes above a horizontal heated cylinder, both numerically and experimentally. They control the Rayleigh number by changing the temperature difference between the cylinder and the fluid, allowing them to explore Rayleigh numbers from $2.05 \times 10^{7}$ to $7.94 \times 10^{7}$ .

To the best of our knowledge, the only study that integrated the analysis of a floating melting object with a theoretical description of the heat transfer is by Hosseini & Rahaeifard (2009). The authors put an ice cylinder to float in a tank, and measured its cross-sectional width constraining it between vertical wires and measuring the distance between the wires. They find that the Rayleigh number (in their case, Ra $=\mathcal{O}(10^7)$ ) scales with the Nusselt number with a 1/3 power law, for experiments conducted in freshwater. We extend their results by analysing images of the cross-section, thus being able to deliver a more complete morphological analysis.

In several natural circumstances, icebergs (or ice chunks, of all sizes) melt in stratified waters. The most studied example is the one of Greenlandic fjords, where icebergs calve from a glacier and accumulate in the fjord, where they melt, releasing cold fresh water, which accumulates on top of of the warmer saline water. The draft of icebergs is various, and with it the fraction of mass immersed in each layer (Jackson, Straneo & Sutherland 2014; FitzMaurice et al. 2016). A similar stratification mechanism applies to open-water icebergs, as described in existing studies on plumes around icebergs (e.g. Helly et al. 2011, Yankovsky & Yashayaev 2014). As described in Yankovsky & Yashayaev (2014), the meltwater from an iceberg accumulates around the iceberg itself, and the iceberg will find itself in a stably stratified two-layer environment, where the less saline and thus lighter top layer has been produced by the accumulation of the iceberg’s meltwater itself.

Recently, the problem of melting has been studied experimentally in several geometries and flow configurations. In FitzMaurice, Cenedese & Straneo (2017) the authors used ice suspended in a flume to study the effect of oceanic currents on the melting, and highlighted the existence of two melting regimes, one with the melt plume attached and one with the melt plume detached. A similar set-up was used in Hester et al. (2021), where the authors compare the melt rate in the different faces of a cuboid iceberg, and show the dependency of the melt rate variations between the faces on the water velocity. The authors of Meroni et al. (2019) carried out experiments with ice melting in a rotating water tank, and found a relation between the Rossby number (Ro, non-dimensionalised rotation speed) and the melt rate. The authors of McCutchan, Meyer & Johnson (2024) report the results for an ice ball that melts submerged in an enclosed tank, with varying water temperatures and turbulence conditions, and show, for quiescent fresh water, the appearance of a well-defined upward or downward plume, depending on the water temperature. In Xu et al. (2024) the authors study the melting behaviour of a submerged ice cylinder whose axis is parallel to gravity, and report on the heat transfers and the scallops that appear on the surface of the cylinder itself. Lastly, Waasdorp et al. (2024) elaborates on the problem of submerged melting of an object composed of a material (olive oil) which is immiscible in water, and derive theoretical descriptions for the boundary layer characteristics. The problem of floating melting was explored by Johnson et al. (2023) in a gallery of fluid motion where they show that floating ice cylinders repeatedly capsize while melting.

In this paper we will focus on three aspects of the floating melting problem. First, how the size of ice and the water’s salinity affect the morphology of the ice. The melt characteristics of the ice are a direct consequence of the local heat transfer with the water, which is a function of salinity, temperature, and velocity fields.

Second, how the ice exchanges heat and momentum with the surroundings, in terms of (non-dimensional) heat transfers and plume characteristics. On these topics, we seek for scaling relations that summarise the physics.

Last, how stable the ice cylinders are against rotations around their axes. Both uneven melting and flow asymmetries can give rise to changes in the rotational stability of the object, which, in turn, causes oscillatory, rotational motion.

This article is organised as follows. In § 2 we provide our experimental methods. In § 3 we develop theoretical interpretations to the morphological evolution, convective heat transfers, plume, rotational stability and rotations. In § 4 we present and analyse our experimental results and put them in the theoretical framework. In § 5 we extend our findings to the more general two-layers stable stratifications. Lastly, in § 6 conclusions and an outlook are provided.

2. Experimental apparatus and procedures

2.1. Experimental set-up

All the experiments were carried out in a glass tank with outer dimensions W $\times$ L $_{\textit{tank}}$ $\times$ H = 80 cm $\times$ 40 cm $\times$ 50 cm and 8 mm-thick walls. The tank was filled with water to a level of 38.3 cm, for a total volume of 115 l. The ice cylinders were put to float parallel to the short-side of the tank, centred with respect to the long side. To avoid the cylinder drifting around in the tank, vertical nylon wires bound the cylinder laterally, see figure 1. The wires also prevent rotations around the vertical (as seen by Dorbolo et al. 2016 and Schellenberg, Newton & Hunt 2023), which would have been detrimental to our backlight illumination set-up. Note that we found that neither the thermal and dynamical interaction between the cylinder and the wire is not important for the melting dynamics, nor that the friction between the ice and the wires prevents possible rotational motion of the cylinders around their axes. During the placement of the cylinders at the beginning of every experiment, care was taken to not excessively disturb the water. Nonetheless, standing and travelling waves were always present for the first few seconds of the experiments. To avoid any parallax problem in the imaging of the near subsurface region for subsequent experiments, an equivalent amount of water was removed from the tank after the melting of every cylinder. Similarly, all the water was mixed after every experiment, to avoid unwanted stratification or inhomogeneities in temperature or salinity. We froze ice cylinders of different diameters: 5, 8.1 and 12 cm, with a length of approximately 30 cm, resulting in a total volume of 0.5 l, 1.5 l and 3.4 l, respectively. A summary of the explored parameters is shown in table 1. The presence of a tiny number of bubbles in the ice (resulting from dissolved air in the water) seemed not to influence the melting process, in line with the findings of Wengrove et al. (2023). Despite this, we used deionised water (Milli-Q), which we put under vacuum with constant stirring for a sufficiently long time before freezing to remove the majority of the dissolved gas. The small (5 cm diameter) and large (12 cm diameter) cylinders were made using aluminium moulds, and the middle-sized (8.1 cm diameter) cylinders in PVC moulds. The freezer was kept at a temperature of ${-}5\,{}^\circ \textrm{C}$ $\pm$ 0.1 K, to ensure a slow freezing process. When frozen, the cylinders were quickly removed from the freezer, and from the moulds, and put to thermalise in another freezer at ${-}16\,{}^\circ \textrm{C}$ . After a day (the time scale for thermal diffusion in the ice ( $r_{{cyl}}^2/\kappa _T$ ) is of the order of 1 hr), they were removed and put to float in the tank.

Figure 1. Sketches of the experimental set-up. An ice cylinder (initial temperature $T_{\textit{ice}}$ , initial radius $R_i$ ) is left free to float on the surface of water (initial temperature $T_\infty$ , initial salinity $S_\infty$ , initial density $ \rho _\infty =\rho (T,S)$ ) in an open tank. (a) Vertical nylon wires create a ‘cage’ around the ice, with the ice being in contact with the wires only occasionally. The dashed circle sketches the equivalent radius of the ice, calculated as the radius of a disk that has the same area as the cross-section of the ice (with a radius $r$ and a perimeter $P=2\unicode{x03C0} r$ ). (b) For contour tracing, we use background illumination from a linearly polarised light-emitting diode (LED) light source from a panel at the back of the ice. A linear polarising filter, with the polarisation direction orthogonal to the light source, is screwed onto the objective of the camera. This configuration of polarisation only lets light through that passes through a birefringent material, which crystalline ice is. The side (non-polarised) panel shines light that illuminates the opaque parts of the ice cylinder, scattering non-polarised light, which is also captured by the sensor of the camera. For flow field measurements, we use planar particle image velocimetry (PIV) with the laser plane being vertical and perpendicular to the cylinder’s axis. The length of the cylinders is indicated with $L$ , and the lateral area as $A_{\textit{lat}}=LP$ .

Table 1. Values of initial Rayleigh numbers (calculated according to (3.7)) for the explored cases of salinity and diameter, and the cases with deuterium oxide. Initial temperature of the water $T_\infty \approx$ ${+}20\,{}^\circ \textrm{C}$ , initial temperature of the ice $T_{\textit{ice}}\approx$ ${-}16\,{}^\circ \textrm{C}$ .

The bulk water temperature was measured at the beginning and end of every melting experiment with a K-type thermocouple (temperature sensitivity 0.1 K). To measure the homogeneity, the temperature was probed in several locations of the tank, and the absolute variation was considered as the error. Note that the typical error was of the order of 0.2 K, for a typical temperature of $19\,{}^\circ \textrm{C}$ .

To vary the relative density between ice and water ( $\varLambda = \rho _{\textit{ice}}/\rho _\infty$ ) we modified our experiment in two ways. To explore ice that would float less, we performed experiments with deuterium oxide (D $_{\text{2}}$ O, colloquially known as heavy water). For those, a mixture of 81 wt % D $_{\text{2}}$ O and 19 wt % H $_{\text{2}}$ O was frozen. The ice of this mixture closely matches the density of the surrounding water, $({\rho _\infty -\rho _{\textit{ice}}}/{\rho _\infty })\approx 1\,\%$ , this made the ice float; just touching the surface, without breaking surface tension.

Similarly, to investigate ice melting in denser water, we added sodium chloride (NaCl) to the bulk water. This not only varies the ratio $\rho _{\textit{ice}}/\rho _\infty$ , but introduces added complexity as the fluid’s density depends on both temperature and salinity ( $\rho _\infty =\rho _\infty (T,S)$ ).

To study the ice in absence of rotational motion, a 3D-printed plastic holder was frozen in one end of the cylinder. The holder can slide along a vertical guide, allowing for vertical translation of the ice without rotation. The size and density of the holder made its impact on the buoyancy of the cylinder negligible.

2.2. Measurements of shape

In order to capture the shape and orientation evolution during the melting, a background illumination set-up was arranged around the tank, as seen in figure 1. Two LED-light panels shine light from the side and back of the cylinder. The one in the back projects linearly polarised light. A DSLR camera (Nikon D850, 45.7 Mpixel, pixel size $4.35\,\unicode{x03BC} \textrm {m}$ ) with a 300 mm objective (Nikon AF-S Nikkor, resulting in a resolution of approximately $50\,\unicode{x03BC} \textrm {m}\,\textrm {pixel}^{-1}$ ) records the shadow of the cylinder. A linear polarising filter is screwed on the objective of the camera, with the axis perpendicular to the one of the LED panel. This set-up is most effective when high contrast is needed, but has not been used for all the experiments. The effect is based on the birefringent properties of clear (crystalline) ice, which changes the polarisation direction of the incoming polarised light, thus letting some component through the two perpendicular polarising filters. On the other hand, the sharp index of refraction changes in opaque ice (with lattice defects and microbubbles), leading to scattered light. Hence, the posterior (polarised) LED panel illuminates clear parts of the ice, while the lateral panel illuminates the opaque parts. The optical axis of the camera was aligned to be perpendicular to the sides of the tank, and the cylinder aligned with the vertical wires to be in line with the optical axis.

Images were recorded at intervals of 60 s for the cylinders of 5 and 8.1 cm diameter, and 180 s for the cylinders of 12 cm diameter, corresponding to the different melting rates. Contours were determined manually, with the aid of an assisted drawing tool (MATLAB drawassisted). Each melting experiment produces 30–40 contours. The contours retain all the information about the two-dimensional shape of the ice cylinder, and the cross-sectional area and equivalent radius can be computed. Investigation on the minor lengthwise variations of the morphology of the cylinder and end effects are out of scope for this work.

2.3. Measurements of rotation dynamics

To investigate the rotational dynamics more precisely, videos of the ice rotations were recorded using the same lighting set-up as described above, with a resolution of $3840 \times 2160$ pixels at a framerate of 30 f.p.s. Of these, every frame was analysed as described in § 2.2.

2.4. Fluid velocity measurements

To investigate the convective flow arising from a melting cylinder we use a PIV set-up, see figure 1. The beam of a Litron LDY-300 Nd:YLF (532 nm) pulsed laser was widened into a vertical sheet with cylindrical optics (LaVision), and crossed the cylinder in a spanwise section, roughly at 1/3 of its length. The water was seeded with rhodamine B-coated poly(methyl methacrylate) (PMMA) particles with a diameter of 1–20 $\unicode{x03BC} \textrm {m}$ from Dantec (PMMA-RhB-10). A Photron SA-X2 camera imaged the flow with a field of view of approximately 15 cm $\times$ 15 cm around the cylinder, at a framerate of 50 Hz and resolution of 1 Mpixel. The synchronisation and triggering of the laser and camera were controlled through LaVision’s DaVis software. The recording time was chosen to be 40 s, resulting in a total of 2000 frames. The camera buffer and data transmission allowed to make a recording roughly every 5 m. The PIV data were analysed using DaVis. Frames were undersampled at one every five (effective 10 Hz). The cross-correlation method was set to multipass (from $32 \times 32$ pixel window to 16 $\times$ 16 pixel window with 50 % overlap), which corresponds to a final window size of 1.22 mm. Note that our PIV measurements were only done in non-stratified water. Furthermore, the regions of interest for our experiments are always outside the thermal BL around the melting ice, such that refractive index matching techniques are not required.

3. Theoretical background

3.1. Morphological evolution: qualitative description

In this subsection we examine the evolution of the morphology and how it is affected by salinity in the background fluid. In the case of floating ice, in contact with air and (saline) water, the heat transfer mechanisms are radiation, conduction and convection. We find radiation to be negligible ( $\sigma (T_{\textit{room}}^4-T_{\textit{ice}}^4)S_{\textit{exposed}}\approx$ 1 W, while convective flows, estimated from $ ({\mathcal{L}\rho _{\textit{ice}}V_{{cyl}}})/{t_{\textit{melt}}}$ , are $\mathcal{O}(10^2$ W).) compared with the other mechanisms. Measured at $20\,{}^\circ \textrm{C}$ , the thermal diffusivity of air is more than 100 times higher than that of water, and the air’s thermal conductivity is roughly 25 times smaller than that of water. These differences imply that the conductive term in the heat balance is dominated by transfers with the water. On top of that, the convective motions of the liquid triggered by buoyancy differences reduce the extent of the thermal boundary layers, further increasing the heat transfer rate.

In the freshwater case, both the meltwater and the cooled down surrounding water are denser than $\rho _\infty$ , hence they sink, and while doing so they exchange heat with the cylinder, getting colder as they sink. The local melt rate depends on the temperature gradient between the ice and the water, and consequently decreases with depth.

Figure 2. Sketch of the boundary layers around an ice cylinder melting in fresh, saline and stratified water. In freshwater (a), a single direction of flow is developed. Both the cold meltwater and the cooled down surrounding water are denser, and sink, forming a plume below the cylinder. In (b i) we see how a sufficiently high salinity separates the flow in two directions, which carry two different density anomalies: the cold, fresh meltwater becomes lighter than the surrounding saline water, and rises, accumulating under the free surface; on the contrary, the cooled down, surrounding saline water becomes denser and sinks, forming a plume. The variation in intensity of the two flows along the cylinder implies that parts of the flow are bidirectional and other are unidirectional, see (b ii). In the case of two-layered (stratified) water (c), the meltwater will accumulate across the density step, either being heavier than the top-layer water or lighter than the dense bottom-layer.

As for the saline case (for the salinities that we explored), the meltwater is less dense than $\rho _\infty$ (despite being much colder, the density is dominated by salinity), and tends to rise. The melting of the cylinder releases fresh meltwater, which further decreases the relative density of this innermost fluid. At the same time, the much wider extension of the thermal boundary layer compared with the solutal boundary layer ( $\delta _{\textit{thermal}}/\delta _{\textit{solutal}}\approx$ 20, as the Lewis number Le $=\kappa _T/\kappa _S \approx 400$ and $\delta _{\textit{thermal}}/\delta _{\textit{solutal}}\approx \sqrt {{Le}}$ ) implies that the bulk water decreases its temperature, and hence increases its density, see figure 2(b i). Just like in the freshwater case, the water continues to get cooled during its sinking motion, thus increasing its relative density. For most of the cylinder’s latitudes, the flow in the vicinity of the ice surface is bidirectional, but there is a point where the positive buoyancy of the meltwater overcomes the negative buoyancy of the cooled saline water, and the flow becomes unidirectionally rising, see figure 2(b ii). Lastly, for the relative high salinities analysed here, the fresh meltwater accumulates below the surface, creating a cold, insulating layer which inhibits the melting.

3.2. Heat transfer mechanisms

A number of parameters play a role in the process of melting of a floating object, both for the fluid and for the melting solid. We use non-dimensional numbers to effectively compare experiments carried out under different configurations. The relevant non-dimensional numbers are the Rayleigh number, describing the intensity of the thermal driving, and the Nusselt number, describing the heat flux as compared with the conductive case.

From the heat balance equation for a melting horizontal right prism, the Stefan condition reads

(3.1) \begin{equation} \frac {1}{P}\rho _{\textit{ice}} \frac {\partial V}{\partial t}\mathcal{L} = \bigg \langle Lk_{\textit{ice}}\frac {\partial T}{\partial n}\biggr \rvert _{-} - Lk_{w\kern-0.3pt\textit{ater}}\frac {\partial T}{\partial n}\biggr \rvert _{+}\bigg \rangle _P\,. \end{equation}

Here $V$ is the volume of the prism, $L$ its length, $P$ its cross-sectional perimeter, $\mathcal{L}$ the latent heat, $k=\kappa _T\rho c_p$ the thermal conductivity of the material, $c_p$ its specific heat, $\kappa _T$ its thermal diffusivity, $ ({\partial }/{\partial n})\big \rvert _{\pm }$ the spatial derivative in the direction normal to the ice surface, pointing into the water, evaluated at the interface, either in the ice (−) or in the water (+). The symbols $\langle \boldsymbol{\cdot }\rangle _P$ indicate the average over the perimeter. The expression equates the heat flux from the water with the sum of those needed to melt the ice and that is diffused in the ice. The contributions of the ends are neglected (which we think is justified, given the aspect ratio of our cylinders).

Given the convective motions that arise in the vicinity of our cylinders, we can express the terms on the right-hand side of (3.1) using an averaged convective heat transfer coefficient $h$ ,

(3.2) \begin{align} \bigg \langle Lk_{\textit{ice}}\frac {\partial T}{\partial n}\biggr \rvert _{-} - Lk_{w\kern-0.3pt\textit{ater}}\frac {\partial T}{\partial n}\biggr \rvert _{+}\bigg \rangle _P &= \frac { h(T_\infty -T_{\textit{melt}})A_{\textit{lat}}}{P}, \end{align}

(3.3) \begin{align} A_{\textit{lat}}&=LP, \end{align}

where $A_{\textit{lat}}$ is the lateral surface of the ice. We assume here that the liquidus temperature $T_{\textit{melt}}$ is equal to $0\,{}^\circ \textrm{C}$ , hence $T_\infty -T_{\textit{melt}}\approx T_\infty$ , which is a good approximation only for the salinity and temperature ranges explored in this work, and not for geophysical applications.

Under the assumption that the melting occurs in the radial direction, the volume evolves as follows:

(3.4) \begin{equation} \frac {\partial V}{\partial t} = L\frac {\partial A}{\partial t} \end{equation}

where $A$ is the cross-sectional area. Putting everything together, we obtain

(3.5) \begin{equation} \rho _{\textit{ice}} \frac {\partial A}{\partial t} \mathcal{L} = h T_\infty P. \end{equation}

If we define a typical length scale for a melting right prism of base area $A$ as the square root of $A$ , and using (3.5) then the Nusselt number can be defined as

(3.6) \begin{equation} \textit{Nu}=\frac {h\sqrt {A}}{k_{w\kern-0.3pt\textit{ater}}}= \frac {\rho _{\textit{ice}} \sqrt {A} \frac {\partial A}{\partial t}\mathcal{L}}{k_{w\kern-0.3pt\textit{ater}} T_\infty P}\,\,\,. \end{equation}

Under these assumptions for the length scales, the expression for the Rayleigh number is

(3.7) \begin{equation} {Ra} = \frac {g\Delta \rho \big(\sqrt {A}\big)^3}{\kappa _T\nu \bar {\rho }} \end{equation}

with $\Delta \rho /\bar {\rho }$ the density difference between the meltwater and the far-field water, normalised by the average between the two densities $\bar {\rho }$ , $g$ the gravitational acceleration and $\nu$ the kinematic viscosity of water.

The expressions for Nu and Ra derived so far hold for any right prism. However, in the case of a right prism with a polygonal base with sharp angles, corrections due to the Gibbs–Thomson effect could be needed.

The densities of the ambient water are calculated with the formula provided in Millero & Huang (2009), which accounts for the temperature and salinity of the water. For the heavy water cases, we estimated the density of our frozen mixture of H $_2$ O and D $_2$ O as a linear combination of the two densities. However, we could not find a valid measure of the density of heavy water ice for temperatures below $0\,{}^\circ \textrm{C}$ , hence we assumed the ice to be at $0\,{}^\circ \textrm{C}$ . (Our best estimate is that the temperature dependency of the density of heavy water D $_2$ O behaves like that of ‘normal’ water H $_2$ O, with a $\approx$ 0.3 % difference between the density at ${-}16\,{}^\circ \textrm{C}$ as compared with that at $0\,{}^\circ \textrm{C}$ .)

3.3. Line-source plumes

As sketched in figure 2, density differences generate downward flows, both in the fresh and saline water experiments. Such flows interact with the quiescent water in the tank, and are expected to show similar dynamics as line-source plumes. Hereafter we adapt classic plume theory (Straneo & Cenedese (2015), and references therein) for fully developed, turbulent line-source plumes, to obtain simple expressions for the plume width and velocity, as a function of depth.

The theory assumes that a localised, two-dimensional source of mass and negative buoyancy generates a region of downward motion in a fluid, where the width $b(y)$ and velocity $w(y)$ can be calculated from the entrainment coefficient $\alpha$ and the source’s volume and buoyancy fluxes ( $Q$ and $B$ , respectively), both per unit length. (The length of our cylinders can be assumed constant as the axial melt rate is much slower than both the typical plume velocity and the radial melt rate.) Under the assumption that the ambient water density is homogeneous, the plume buoyancy flux $B$ is predicted to be constant with depth. For the plume width $b$ , the theory predicts a linear relation with depth $y$ , with a slope equal to the entrainment coefficient $\alpha$ , thus $b=\alpha y$ . In our case, this expression will have with a virtual origin correction $\tilde {y}$ ( $\tilde {y}$ such that $w(\tilde {y})=0$ (Cenedese & Linden 2014)) from our cylinders having a non-zero horizontal extent. For the plume velocity, the theory predicts the expression

(3.8) \begin{equation} w = \left (\frac {B}{2\alpha } \right )^{\frac {1}{3}} \end{equation}

that is constant in depth under the assumptions stated before. Additionally, given the definition and theoretical prediction for the reduced gravity in the plume $g'(y)$ ,

(3.9) \begin{equation} g'(y):=g\frac {\rho _p(y)-\rho _\infty }{\rho _0}=\left (\frac {B}{2\alpha }\right )^{\frac {2}{3}}\frac {1}{y} \end{equation}

where $\rho _p(y)$ is the density in the plume, $\rho _\infty$ the ambient density and $\rho _0$ a reference density (taken here as the mean of $\rho _p$ and $\rho _\infty$ ), the plume velocity can then be rewritten as

(3.10) \begin{equation} w = \sqrt {yg'(y)}=\sqrt {yg\frac {\rho _p(y)-\rho _\infty }{\rho _0}} \,\, , \end{equation}

highlighting the relation between the velocity and the plume density.

To compare the velocity of the plumes generated by experiments with different Rayleigh number (3.7), we define the Reynolds number as

(3.11) \begin{equation} \textit{Re} = \frac {w\sqrt {A}}{\nu } \end{equation}

with $w$ the typical plume velocity, and $A$ the cross-sectional area. Note that the Rayleigh number, the driving parameter of the problem, is based on the cylinder, while the Reynolds number, the response parameter of the problem, is based on the velocity of the plume.

Figure 3. Gravitational and buoyancy effects on a floating melting cylinder (D $_{\textit{equivalent}}$ = 5.9 cm, D $_{\textit{initial}}$ = 8.1 cm, T $_\infty$ = $19.3\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, S = 0.0 g l⁻¹, Ra $_{\textit{initial}}$ $\approx$ $4 \times 10^{7}$ , t $\approx$ 10 min). The dashed horizontal line is the water level. The blue and orange dots (visible also in the insets) indicate the COM and COB, respectively. In the inset, the distance $d$ is the horizontal distance between the two centres. The inset shows that in an equilibrium position the COM and COB are vertically aligned (up to experimental accuracy), but the COB gets displaced if the shape is rotated around the COM. Note that the horizontal scale is 10x stretched to highlight the precision of the experiment and to visualise the horizontal distances easier.

3.4. Conditions for stability of cylinders

At melting temperature, the density of liquid water is 9 % higher than that of crystalline ice, so whenever left free to float, ice exposes approximately 9 % of its volume to air (the literal tip of the iceberg). The higher melt rate of the immersed ice compared with the ice exposed to air (see § 3.1), makes the ice bottom-heavy (with gravity effectively acting upwards). In the case of an ice cylinder, the stability dynamics will manifest itself through rotations about the cylinder’s axis of symmetry. The cylinder (or any irregular right prismatic shape) will rotate around its axis until a stable equilibrium is found.

Any floating object is subject to gravity and buoyancy, which can be considered as acting on two ideal points, the centre of mass (COM) and centre of buoyancy (COB). If the object is in an equilibrium position, the COM and COB are vertically aligned. When a torque is applied, the induced rotation horizontally displaces the COB, which generates a torque that can stabilise or destabilise the object, see the sketch in figure 3. For rotations happening around the COM, the arm of gravity’s torque is zero. We call $d$ the horizontal displacement between COM and COB ( $d=x_{\textit{COM}}-x_{\textit{COB}}$ ), see figures 3 and 4. Given our signs convention, in a stable equilibrium, an increase in the rotation angle will produce a negative $d$ , with the corresponding buoyancy torque that will act in the opposite direction as the angle increase, thus stabilising the shape. On the contrary, in an unstable equilibrium, an increase in the rotation angle will produce a positive $d$ , with the buoyancy torque acting in the same direction as the angle increase, thus destabilising the shape.

It is important to notice that in order to destabilise a shape, the angular displacement needs to reach at least the closest unstable equilibrium point, otherwise the restoring force will push back the shape to where it was.

Two mechanisms control the stability of the ice against rotations: the magnitude of the displacement $d$ (which is responsible for the magnitude of the torque) and the angular distance to the closest unstable equilibrium point. This argument is expanded in Appendix B.

With the aim of finding a quantitative measure of the stability which accounts for both these processes, we defined the stability index (SI) as the angular distance to the closest unstable equilibrium point multiplied by the maximum horizontal distance between the COM and the COB throughout a full rotation of the shape, as follows:

(3.12) \begin{equation} {{SI}(\theta ) = (\theta _{\textit{closest}\,\textit{unstable}}-\theta )\, \max _\theta (d(\theta ))\, g\,\rho _{\textit{ice}}\, A} \end{equation}

with $A$ the cross-sectional area. The choice of using $\max _\theta (d(\theta ))$ arises from the fact that in both a stable and an unstable equilibrium point, the COM and COB are vertically aligned, so $d(\theta )=0$ , and this would contrast with one being more stable than the other.

Note that we do not non-dimensionalise the stability index for two reasons: first, applying a uniform scaling to the cylinders changes the magnitude of the forces and torques involved; and second, the stability that the cylinders reach just before a rotation cannot be assumed universal, as the magnitude of the (random) perturbations that trigger a rotation is not. Hence, to maintain physically meaningful dimensions, we multiply this by the buoyancy force per unit length. Consequently, our stability index has the dimensions of a torque per unit length ( $[N]$ ).

Figure 4. Horizontal displacement $d$ of COM and COB ( $d$ $=\text{COM}-\text{COB}$ , along the horizontal axis) as a function of the rotation angle $\theta$ for the shape in figure 3. The points where the curve crosses zero are equilibrium points. If the first derivative of the curve is positive (negative) in an equilibrium point, the buoyancy-induced torque is acting in the same (opposite) direction of the angle perturbation, hence the equilibrium point is unstable (stable). Stable equilibrium points are indicated with green triangles, while unstable with red circles.

3.5. Rotation dynamics

We now develop a dynamical model to describe the rotation dynamics of a capsizing event.

We modelled the rotational motion as a damped, rotating nonlinear oscillator. Newton’s second law for a rotating object states

(3.13) \begin{equation} {I \ddot {\theta } \hat {\boldsymbol{e}}_z= \sum \boldsymbol{\tau } = \sum \boldsymbol{a} \times \boldsymbol{F},} \end{equation}

with $I$ being the moment of inertia of the object, and $\sum \boldsymbol{\tau }$ the sum of the torques exerted on the object ( $\boldsymbol{a}$ being the arm of each force $\boldsymbol{F}$ ). It is assumed that the mass of the displaced fluid due to the rotational motion is negligible compared with the mass of the cylinder, hence no added mass is considered in our model. Consequently, under the assumption of cylindrical cross-section the moment of inertia is simply $I=({1}/{2})mr^2$ , with $r$ the equivalent radius. The deviation from the actual moment of inertia (calculated from the cross-section) is of the order of 1 %. The cylinder is rotating around its COM under the effect of the buoyancy force $\boldsymbol{F}_{\!B}$ , that acts at a distance $d$ from the COM (note that $d = d(\theta )$ , see figure 3). Therefore, $\boldsymbol{F}_{\!B} = -\rho _\infty \boldsymbol{g} V_{\textit{sub}}$ , with $V_{\textit{sub}}=(\rho _{\textit{ice}}/\rho _\infty )V_{{cyl}}$ . The corresponding torque will be $\boldsymbol{\tau }_{\!B} = \rho _\infty V_{\textit{sub}} \boldsymbol{d}(\theta ) \times \boldsymbol{g}$ . Lastly, the cylinder is subject to inertial drag, and as such to a drag force that acts with an arm equal to the radius of the cylinder, with a resulting torque of $\boldsymbol{\tau }_D = -({1}/{2})A_{\textit{lat}}\rho _\infty C_D r^2\dot {\theta }|\dot {\theta }| \boldsymbol{r} \times \hat {\boldsymbol{e}}_\theta$ , with $A_{\textit{lat}}$ being the lateral surface of the cylinder and $C_D$ a drag coefficient.

The resulting ordinary differential equation for the damped nonlinear oscillator thus reads

(3.14) \begin{align} \frac {1}{2}\:\varLambda \:V_{{cyl}}\:r^2\:\ddot {\theta }&=-g\:d(\theta )\:\varLambda \:V_{{cyl}}-\frac {1}{2}\:A_{\textit{lat}}\:C_D\:r^3\:\dot {\theta }\:|\dot {\theta }|\,, \end{align}

(3.15) \begin{align} \ddot {\theta } &= -\gamma \:d(\theta ) - \frac {2\:C_D}{\varLambda } \:\dot {\theta }\:|\dot {\theta }| , \end{align}

with the density ratio $\varLambda = {\rho _{\textit{ice}}}/{\rho _\infty }$ and $\gamma = {2g}/{r^2}$ .

4. Experimental results

4.1. Morphological evolution

First, we will compare the morphological dynamics for the case of saline water to that of fresh water. We observed that ice cylinders melting in fresh water repeatedly capsized, regardless of their initial size, while the ones melting in saline water did not. Because of these different rotational dynamics, we performed an experiment constraining a cylinder in fresh water such that rotations would be prevented, but vertical motion would be allowed. Figure 5 shows the evolution for both cases, and we readily see major differences in the morphological dynamics. Note here that in the freshwater case, the meltwater sinks, while in the saline water case, the meltwater floats. The constrained freshwater case shows an enhanced melt rate in the subsurface region and in the south pole; the saline water cylinder shows reduced melt rate in the subsurface region, and the appearance of a marked minimum of melt rate (showing as a protrusion) that migrates towards the south with time. These morphological features can be explained with arguments based on the boundary layer characteristics in the different cases, see figure 2.

Figure 5. Shape evolution of (a) a cylinder in fresh water whose rotations where prevented (D $_{\textit{initial}}$ = 8.1 cm, Ra $_{\textit{initial}}$ $\approx$ $4.8 \times 10^{7}$ , T $_\infty$ = $20.3\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, S = 0.0 g l⁻¹) and (b) a freely floating cylinder (D $_{\textit{initial}}$ = 8.1 cm, Ra $_{\textit{initial}}$ $\approx$ $1.7 \times 10^{8}$ ) melting in saline water (T $_\infty$ = $18.5\,{}^\circ \textrm{C}\,\pm$ 0.2 K, S = 10.0 g l⁻¹). Contours are spaced by 90 s. Note that the two cylinders are at a different height as a result of the different $\rho _\infty$ .

The feature that appears at the south pole of the freshwater cylinder is an effect of the recirculation zone: for high enough Rayleigh numbers, the sinking flow detaches from the cylinder, creating a turbulent recirculation zone below the object that enhances mixing and heat transfer; in later times (with lower Rayleigh numbers), the flow laminarises, but the footprint of the recirculation zone remains visible. This does not appear in saline conditions because the effective Rayleigh number is decreased as a result of the bidirectional flow, as depicted in figure 2(b ii).

For the case of saline water, the minimum of melt rate is the morphological feature that arises because of the flow inversion identified in figure 2(b ii) and also noted by Yamada et al. (1997). Lastly, the reduced subsurface melt rate is a consequence of the accumulation of cold meltwater at the free surface.

We now come to the more general, unconstrained ice melting in freshwater. In this case, the cylinder capsizes. The shape evolution is shown in figure 6. The different conductive heat transfers (air and water) lead to a marked asymmetry in melt rate, with the top of the ice melting significantly slower than the bottom. This sculpting process can create a horizontal asymmetry in the distribution of mass above and below the water surface, that can result in an imbalance between buoyancy and gravity, thus making the ice unstable against rotations around its axis. The continued erosion makes the ice repeatedly capsize. The interaction between the flow and the ice surface sharpens the initial cylinder to a roughly polygonal prismatic volume.

Figure 6. Shape evolution of an ice cylinder (D $_{\textit{initial}}$ = 8.1 cm, Ra $_{\textit{initial}} \approx 5.1 \times 10^{7}$ ) melting in fresh water (T $_\infty$ = $20.9\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, S = 0.0 g l⁻¹). Each panel is labelled with the time from the beginning of the experiment. The variation of the shape of the cross-section along the cylinder is negligible.

4.2. Convective heat transfers

The scaling of the Nusselt number against the Rayleigh number is provided in figure 7. The data follow a convective Nu $\propto$ Ra $^{1/3}$ scaling that is reported in the literature for laminar-type thermal boundary layers (Churchill & Chu 1975; Yamada et al. 1997; Grossmann & Lohse 2000; Wells & Worster 2008; Hosseini & Rahaeifard 2009). The data follow the 1/3 scaling relation regardless of shape or initial size. A feature appearing in all the panels is the initial increase of Nusselt number. We believe this is due to the temperature profile in the ice becoming less steep as temperature diffuses in the solid.

In the later stages of the melting process, the saline water measurements show a drastic decrease of Nusselt number. As described previously, under saline conditions, the meltwater rises and creates an insulating layer which hinders melting. In terms of the Nusselt versus Rayleigh dependence, this results in a lower Nusselt compared with the predicted 1/3 scaling. The intensity of this deviation depends on the initial size of the cylinder, and not on the Rayleigh number: at a given Rayleigh number, cylinders that have started larger will have melted a larger part of their volume, and consequently the thickness of the subsurface meltwater layer will be larger. Moreover, the mixing between the meltwater and the saline (surrounding) water depends on the salinity itself. This is visible in the early-time deviation from the scaling relation for the 35 g l⁻¹ experiments compared with the 10 g l⁻¹ ones. The deviation shown in few curves in the 0 g l⁻¹ panel, at around $10^{5}\lessapprox {Ra} \lessapprox 10^{6}$ can be ascribed to the transition from a convection dominated regime (Nu $\gg 1$ ) to a diffusion dominated regime (Nu $={O} (1)$ ). The heavy water experiments (D $_2$ O) show a similar Nu versus Ra scaling relation as the others. We consider the top (convex) curve less reliable because, in order to reduce the consumption of the expensive heavy water, we choose to image the melting together with the PIV data. Hence our imaging set-up was not designed for contour recognition.

Figure 7. Scaling of the Nusselt number as a function of the Rayleigh number. The four panels refer to three different salinity cases plus the density matched (heavy water) case. Each experiment is represented by a line with markers. Squares indicate an initial diameter of 5.0 cm, triangles of 8.1 cm and circles of 12.0 cm. As time progresses, the cylinder shrinks, so time direction is from higher Ra to lower Ra. The dashed lines indicate a scaling of Nu $\propto {Ra}^{1/3}$ . The initial increase of Nusselt is explainable by the initial heat diffusion inside the ice. The decrease of the Nusselt number in the last stages of the melting in saline conditions is due to the accumulation of cold fresh water under the surface of the ice.

4.3. Plume measurements

From the PIV measurements, we identified the downward plume by thresholding the vertical velocity. The threshold was chosen of the same magnitude as the converging horizontal flows. Different threshold values do not produce significantly different results. We located the bottom of the ice from the contour information derived from the imaging and considered the velocity field as starting below the ice. Figure 8(a) summarises the characteristics of the plume’s velocity field for a cylinder (initial diameter 50 mm, salinity 35 g l⁻¹). The behaviours described hereafter are consistently seen in all experiments (with different initial size and different salinities).

Figure 8. Characteristics of the plume generated by a melting cylinder (D $_{\textit{equivalent}}$ = 2.04 cm, D $_{\textit{initial}}$ = 5.0 cm, T $_{w\kern-0.3pt\textit{ater}}$ = $17.9\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, S = 35 g l⁻¹, Ra $_{\textit{initial}}$ $\approx 2 \times 10^8$ ). Panel (a) shows the thresholded downward velocity field as obtained from PIV data, with superimposed: streamlines; contour of the ice; water level (horizontal dashed line); the depth at which the profile in (d) has been extracted (red horizontal strip). Panel (b) shows the typical velocity of the plume as a function of depth. The vertical red lines are the velocity that is expected when considering a constant density anomaly of 1 %, 2 % and 3 % above the reference density. Panel (c) shows the width of the plume with a linear fit to the data. In (a), (b) and (c) the dark shading refers to the region outside our region of interest. Panel (d) shows the velocity profile in the plume at the cross-section identified in (a). The horizontal red line is the measured plume width.

As seen in figure 8(a), the actual source of the plume is spread both horizontally and vertically. We determined the velocity profile of the plume by taking the horizontal mean of the velocity field. At a certain depth, the width of the plume was defined as the one that would conserve the downward momentum at that depth, while assuming a top-hat profile. Figure 8(d) shows a typical velocity profile, measured at the transect indicated in red in figure 8(a). Figure 8(c) shows that the plume width increases linearly with depth, in agreement with our line-source plume theory. By performing a linear fit between plume width and depth ( $b=\alpha y$ ) we find an entrainment coefficient of $\alpha = 0.11$ .

A precise measure of the plume density difference as a function of depth is not feasible. A reasonable estimate is obtained by assuming that at the highest point of the plume (just below the ice) the density of the plume was of saline water cooled to zero degrees Celsius. In this case, the density anomaly compared with ambient-temperature saline water is around 3 %. In figure 8(b) we plot the expected plume velocities for density anomalies of 1 %, 2 % and 3 %, calculated with (3.10). We remark that the theory that we presented in § 4.3 assumes a top-hat profile for the plume’s velocity, which is to be compared with the example profile of figure 8(d). This adds on top of the relatively short distance from the source at which our measurements are taken, which relates to the assumption of considering our plumes to be fully developed and turbulent. With these considerations, our measurements show a marked acceleration of the plume with greater depths, and the velocity does not differ significantly from the expected value.

Figure 9(a) plots the Reynolds number as a function of the Rayleigh number of all of our experiments. Straight-line fits of our data are assisted with a Re $\propto$ Ra $^{1/2}$ scaling for the lower region of Rayleigh number (up to Ra $\approx 10^7$ ) and a Re $\propto$ Ra $^{1/3}$ scaling for higher Ra. This is supported by figure 9(a,b). Given that Re depends on $w$ and $R$ , and Ra has its strongest dependency on $R^3$ , a 1/3 scaling implies a constant downward velocity. This is confirmed in figure 9(d), which shows a constant downward velocity for the range of Ra where the 1/3 scaling is valid. Indeed, a roughly constant downward velocity is predicted by (3.8): given that the downward velocity $w$ depends only on the density difference, and the density difference varies only of a few per cent throughout our experiments as a result of different salinities. The equations thus predict a velocity that varies only up to few per cent, much below our experimental accuracy. The scaling relations for Re are in accordance with the Grossmann–Lohse theory (Grossmann & Lohse 2000, 2001), which indeed features a transition from a regime $I_u$ (dominated by boundary layer dissipation, for velocity boundary layer larger than thermal boundary layer) to a regime $IV_u$ (dominated by bulk dissipation, for velocity boundary layer larger than thermal boundary layer).

Figure 9. Reynolds number against Rayleigh number for all the experiments from our PIV data. Each point refers to a PIV dataset (acquisition time 20 s). The colour of the marker (blue, orange or green) refers to the salinity of the surrounding water in the experiment, the shape of the marker (square, triangle, disc) to the initial radius of the cylinder. Panel (a) reports the uncompensated plot, with the dash–dotted line indicating a Re $\propto {Ra}^{1/2}$ scaling and the dashed line a Re $\propto {Ra}^{1/3}$ scaling. Panels (b) and (c) present the compensated plots for these two scaling relations. Panel (d) reports the mean downward velocity (see figure 8 a) as a function of Ra.

Figure 10. Temporal evolution of the stability index for three cylinders (for all, D $_{\textit{initial}}$ = 8.1 cm). They melted in water with salinity S = 0 (circles, T $_{w\kern-0.3pt\textit{ater}}$ = $20.9\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, Ra $_{\textit{initial}}$ $\approx 5.1\times 10^7$ ), S = 10 g l⁻¹ (triangles, T $_{w\kern-0.3pt\textit{ater}}$ = $19.0\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, Ra $_{\textit{initial}}$ $\approx 1.7\times 10^8$ ) and S = 35 g l⁻¹ (squares, T $_{w\kern-0.3pt\textit{ater}}$ = $19.3\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, Ra $_{\textit{initial}}$ $\approx$ $6.9\times 10^8$ ). Vertical dashed lines indicate the times of rotation of the freshwater cylinder (the others do not rotate). For the freshwater cylinder, the stability index increases during a rotation. For the saline water cylinders, the stability index increases monotonically over time.

4.4. Time evolution of the stability of cylinders

Figure 10 reports the temporal evolution of the stability index, as defined in § 3.4, for three cases: fresh water S = 0 g l⁻¹, saline water with S = 10 g l⁻¹ and saline water with S = 35 g l⁻¹. The times of rotations of the freshwater cylinder are marked as vertical dashed lines (the saline cases do not rotate). First, we can support the statement that it is energetically favourable for the cylinder to rotate – i.e. that a rotation lowers the gravitational potential energy of the cylinder – by noting that in the fresh water case the stability index shows a well-defined increase across rotations. (In some rotation events, the increase in stability index is less prominent than in others. This is due to two factors. In some cases the angular rotation is very small, and the vertical displacement of the COM is minimal. In some other case, the cylinder ‘overshoots’, surpassing the closest stable equilibrium point and reaching a further one, which may have a similar stability index.) Second, the stability index of the saline water cases shows a clear increase over time. This is consistent with our experimental observations of the absence of rotations in the saline water case, and supports the statement that the morphology evolution progressively stabilises the cylinder. Most notably, the evolution of the stability of the cylinder is not connected to the flow dynamics, but rather only to the balance between gravity and buoyancy.

4.5. Rotation dynamics

A sequence of images of a typical rotation of a cylinder melting in freshwater is shown in figure 11. Figure 12 shows the angle of orientation with respect to the vertical of an ice cylinder during a capsizing event. We fit the data with (3.15), with $\gamma$ , $C_D$ and $\dot {\theta }(t=0)$ being the fitting parameters. Any experimental error in the measurement of $d(\theta )$ or inaccuracy in the density of the ice or water are accounted for in the errors of $\gamma$ and $C_D$ . A first estimate for $C_D$ is in the range 0.5–1, though a proper estimation is impossible due to the unique shape that the ice has. The equivalent radius $r$ of the cylinder can be evaluated from the contour ((29.5 $\pm$ 3) mm in the case of figure 12), leading to a $\gamma$ equal to (2.1 $\pm$ 0.02) $\times$ 10 $^4$ $\textrm{m}^{-1}\,\textrm{s}^{-2}$ . The fit returns a $C_D$ of 0.41 and a $\gamma$ of (1.9 $\pm$ 0.01) $\times$ 10 $^4$ $\textrm{m}^{-1}\,\textrm{s}^{-2}$ . We believe that the mismatch in the $\gamma$ value is mainly due to the assumption of a constant cross-section of the cylinder along its axis, which is directly influencing the accuracy of the radius and ${d}(\theta )$ .

Figure 11. Images of an ice cylinder (d $_{\textit{equivalent}}$ = 5.8 cm, d $_{\textit{initial}}$ = 8.1 cm, T $_{\textit{fluid}}$ = $19.6\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, S = 0.0 g l⁻¹, Ra $\approx$ $1.6 \times 10^{7}$ ) performing an anticlockwise rotation, shown as an example of the mechanism. The white arrow is added to the image to help the reader track the motion of the cylinder. The images are taken on a later stage of the melting, and the ice is already not circular anymore. The vertical lines visible in the image are the nylon wires to keep the ice in place. The ice is rotating and then oscillating around the new stable equilibrium.

Figure 12. Best fit of the solution of (3.15) to the data of a rotational oscillation of a cylinder (D $_{\textit{equivalent}}$ = 5.9 cm, D $_{\textit{initial}}$ = 8.1 cm, T $_\infty$ = $19.3\,{}^\circ \textrm{C}$ $\pm$ 0.2 K, S = 0.0 g l⁻¹, Ra $_{\textit{initial}}$ $\approx$ $4\times 10^7$ ). The blue dots refer to the angle of rotation around the COM, where the zero is set to the initial position of the ice. The red line is the result of fitting the parameters of (3.15) ( $\gamma$ , $C_D$ and $\dot {\theta }_{t=0}$ ) to match the experimental data.

5. Extension to general two-layer stratifications

In our saline water experiments, the melting cylinder created a two-layer density stratification in the water, which was initially at constant density. In nature, however, there are situations where floating ice melts in already-stratified water, like in ice melanges (see Burton et al. 2018 and references therein), or pack ice (Rothrock 1975), or in fjords (Jackson et al. 2014; FitzMaurice et al. 2016). A common feature of these examples is that in an enclosed basin multiple icebergs melt at the same time. The case of multiple ice blocks melting at the same time goes beyond the scope of this work, but we can easily access the case of an initial two-layer density stratification in the water.

We have conducted two sets of experiments. All of them with the 81 mm diameter cylinders at an initial temperature of ${-}16\,{}^\circ \textrm{C}$ . The first set was such that the fresh water layer left by a molten cylinder was used as the initial condition for the subsequent one. In our tank, the volume of a cylinder corresponds to a water height of approximately 5 mm. Thus, the first cylinder melted in saline water (S = 35 g l⁻¹), the second one with a fresh water layer on top of 5 mm, the third one with 10 mm, the fourth with 15 mm and so on. The temperature of the top layer was not controlled, but it was determined by the influx of meltwater at $0\,{}^\circ \textrm{C}$ , and the heat transfers with the air above and the saline water below. The second set of experiments was done with the same layer thicknesses, but the fresh water layer was created every time with water at $24\,{}^\circ \textrm{C}$ . The sharp stratification was created by carefully pouring water on a floating sponge, that would break the vertical momentum of the poured water. In this configuration, the meltwater accumulates in between the bottom saline layer and the top warm freshwater layer (see figure 2 c).

As visible in figure 13(a), all the first set of experiments resulted in rotationally stable cylinders. In fact, the melting was primarily happening in the bottom saline layer, and the resulting shape was stable. In addition to this, the melting time increased substantially, given that – after the melting of the lower part of the cylinder – the remaining part was immersed in cold freshwater. On the contrary, all the second set of experiments resulted in unstable cylinders (an example visible in figure 13 b). Note here that it was experimentally impossible to artificially create a sharp layer of warm freshwater thinner than 5 mm. Although the morphological changes of these experiments differ from the ones seen in the full fresh water, the reasons for the rotational instability still lie in the horizontal displacement of buoyancy and gravity.

Figure 13. (a) Snapshots of subsequently melting ice cylinders (for all, D $_{\textit{initial}}$ = 8.1 cm). The red dashed line indicates the water surface. The first cylinder was melted in saline water (S = 35 g l⁻¹, T = $19.8\,{}^\circ \textrm{C}$ $\pm$ 0.2 K). The time indicated at the bottom of each panel is the time elapsed from the beginning of each experiment. The number following the hash symbol is the progressive number of the experiments. The meltwater accumulating on top of the denser saline water created a two-layer stratification, with the top layer being fresh and cold and the bottom layer being the original ambient water. The white dashed line indicates the extent of the meltwater layer. Each following experiment was melted in the stratification as left by the previous experiment. The blue shade on the picture is due to food colouring used to chromatically distinguish the bottom layer from the top one. (b) Two shots of an ice cylinder (D $_{\textit{initial}}$ = 8.1 cm) melting in two-layered water, with the top layer being fresh and the bottom layer being at salinity S = 35 g l⁻¹. The temperatures of the two layers were T $_{\textit{top}}$ = $20.6\,{}^\circ \textrm{C}$ , T $_{\textit{bottom}}$ = $24.2\,{}^\circ \textrm{C}$ . The two shots correspond to two instants before and after a rotation of the cylinder. The blue shade on the picture is due to food colouring used to chromatically distinguish the bottom layer from the top one. The red dashed line is the water surface, the white dashed line is the pycnocline. The thickness of the top layer is the lowest for which we can ensure a visually sharp distinction between the two layers.

6. Conclusions and outlook

We have conducted melting experiments in a 115 l tank filled with quiescent water, where ice cylinders were put to float with their axes perpendicular to gravity. We have changed both the density of the ambient water and that of the melting object. The salinity of the water was changed from fresh to 10 g l⁻¹ to oceanic salinity (S = 35 g l⁻¹). We performed a set of experiments with deuterium oxide (heavy water) to explore differently buoyant cylinders. We monitored the morphological evolution and vertical orientation of the cylinders and imaged the convective flow around the cylinder. Our work expands the fundamental research on floating, melting objects. This is of particular significance, given that, in several natural and industrial circumstances, floating ice particles are not fixed in place.

The interaction of the water with the cold cylinder gives rise to density anomalies which drive convective flows in the tank. According to its own density anomaly, the meltwater can either mix or shear with the outer flow. We identified the key physical mechanisms that determine the morphological evolution of the cylinders. Salinity adds a third (solutal) boundary layer to the existing two (thermal and momentum), with the three sharing no common length scale. For sufficiently high ambient water salinities, the flow in the vicinity of the ice can transition from being bidirectional closer to the bottom of the cylinder, to unidirectional closer to the air–water interface. The local melting of the ice by the flow is, consequently, drastically different in saline waters as compared with fresh waters.

We found that the Nusselt number (non-dimensional heat transfer) against the Rayleigh number scales as Nu $\propto$ Ra $^{1/3}$ . This scaling is independent of size and salinity, and aligns with the previous literature findings.

We studied the sinking plume that develops below the cylinders and compared its dynamics with that of line-source plume. Despite being far from an ideal plume, the velocity and width profiles are consistent with the predictions with the theory for a line source plume. The Reynolds number of the velocity of the line source plume was plotted against the cylinder’s Rayleigh number. The data collapses on a line independently of size and salinity, and aligns with known theory for plumes (linearly increasing plume width with depth and constant plume velocity) and convective heat transfers (Re $\propto$ Ra $^{1/2}$ for Ra $\lt \mathcal{O}(10^7)$ and Re $\propto$ Ra $^{1/3}$ for Ra $\gt \mathcal{O}(10^7)$ ).

Cylinders melting in fresh water experience repeated capsizing events, which result from an asymmetric melt rate between the part exposed to air and the part exposed to water. On the contrary, cylinders melting in saline water have a cross-section that is highly stable against rotations, hence do not capsize. Based on Newton’s second law for rotating objects, we developed a measure of the stability of the cylinders, which qualitatively describes the observed phenomena. We modelled the oscillations that occur after a capsize via a damped nonlinear oscillator model, which fits our data satisfactorily.

Finally, we conducted experiments in two-layer stratified water, and showed that the ice is stable when the melt rate is higher in the bottom layer, and unstable when the melt rate is higher in the upper layer. We argue that our experimental findings in the laboratory tank are qualitatively connected to several natural circumstances where icebergs melt in stratified waters, with the substantial difference that real icebergs overturn in saline waters, while our ice cylinders do not. In fact, a few caveats apply to this analogy.

First, the density anomaly that we observed in the laboratory is of the same order as the density difference between fresh and salty water (25 to 30 g l⁻¹), while in observations they report much lower anomalies, of the order of 0.5 g l⁻¹. This is due to the different melting time scales between our experiments and real icebergs. In fact, the ocean temperature at the poles can be as low as zero degrees , and the melting process in such conditions is greatly slowed down, hence the mixing time scales are close to the melting time scale. On the contrary, our ice melts much quicker than the mixing time scales.

Second, the typical thickness of the naturally occurring layer is both hard to determine and impossible to repeat in a laboratory. Yankovsky & Yashayaev (2014) report significant density stratifications up to $-40\,\textrm {m}$ of depth, but they report lack of data in the upper $14\,\textrm {m}$ layer. Direct measurements of icebergs’ keel depth are not usually carried out during oceanographic cruises, but Dowdeswell & Bamber (2007) proposed an estimate of the keel’s depth based on satellite-altimetry data. Their frequency distribution presents three peaks: one at 150–200 m, one around 250–300 m and one at 550–600 m. Detailed values depend on the region of survey. Hence, with the information of a typical plume depth of the order of tens of metres, and a typical iceberg keel depth of the order of hundreds of metres, we conclude that our centimetre thick layer with a 10 cm-thick ice is proportionally not too far off from typical natural occurrence.

Third, a number of mechanisms occur only in nature and not in our laboratory, such as wave erosion and solar radiation. These two typically act on the top layer of an iceberg, resulting in an increased ablation rate of the ice close to the surface of the water. This ablation rate is significantly different from the one of our ice cylinders that melt in well-mixed saline water, but is similar to the one of cylinders melting in stratified water (with the top part melting quicker than the bottom part). For these reasons, the capsizing behaviour of these ice cylinders can be connected with the one of icebergs. Another feature that only appears in the natural icebergs is the presence of (unevenly distributed) sediments and large pockets of air and air bubbles. On top of altering the ablation dynamics, these also generate density inhomogeneities that could determine a significant difference in the stability of real icebergs compared with our cylinders.

In our future experiments, we will vary the temperature of the surrounding water, to lower the corresponding contribution to the melt rate.