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. Reference Dorbolo, Adami, Dubois, Caps, Vandewalle and Darbois-Texier2016 and Schellenberg, Newton & Hunt Reference Schellenberg, Newton and Hunt2023), 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. (Reference Wengrove, Pettit, Nash, Jackson and Skyllingstad2023). 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.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 (Reference Millero and Huang2009), 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 (Reference Straneo and Cenedese2015), 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 Reference Cenedese and Linden2014)) 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.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. (Reference Yamada, Fukusako, Kawanami and Watanabe1997). 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 Reference Churchill and Chu1975; Yamada et al. Reference Yamada, Fukusako, Kawanami and Watanabe1997; Grossmann & Lohse Reference Grossmann and Lohse2000; Wells & Worster Reference Wells and Worster2008; Hosseini & Rahaeifard Reference Hosseini and Rahaeifard2009). 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 Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001), 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.