Enhanced bubble growth near an advancing solidification front

Abstract

Frozen water might appear opaque since gas bubbles can get trapped in the ice during the freezing process. They nucleate and then grow near the advancing solidification front, due to the formation of a gas supersaturation region in its vicinity. A delicate interplay between the rate of mass transfer and the rate of freezing dictates the final shapes and sizes of the entrapped gas bubbles. In this work, we experimentally and numerically investigate the initial growth of such gas bubbles that nucleate and grow near the advancing ice front. We show that the initial growth of these bubbles is governed by diffusion and is enhanced due to a combination of the presence of the background gas concentration gradient and the motion of the approaching front. Additionally, we recast the problem into that of mass transfer to a moving spherical object in a homogeneous concentration field, finding good agreement between our experimental data and the existing scaling relations for that latter problem. Lastly, we address how fluid flow around the bubble might further affect this growth and qualitatively explore this through numerical simulations.

1. Introduction

While slowly freezing an aqueous suspension, the dispersed objects might either be engulfed into the ice or rejected by the advancing solidification front (Shangguan, Ahuja & Stefanescu 1992; Rempel & Worster 2001; Dedovets, Monteux & Deville 2018; Tyagi et al. 2020; Meijer, Bertin & Lohse 2023a). The conditions that govern this are of significant importance in material science as the arrangement of these dispersed particles profoundly shapes the microstructure and, consequently, the functional attributes of the solidified material (Deville et al. 2006; Deville, Saiz & Tomsia 2007; Deville 2008, 2010). Especially during repulsion of objects by the moving front, typically occurring at low freezing velocities, the objects tend to accumulate (Tyagi, Monteux & Deville 2021, 2022), forming a concentration profile of the objects that can evolve in space and time. Similarly, since gases are soluble in a large variety of liquids, e.g. water, silica (Yokokawa 1986), metals (Shapovalov & Boyko 2004) and sapphire (Bunoiu, Duffar & Nicoara 2010; Ghezal et al. 2012), but not in their solid states, during the solidification process, accumulation of gases at the front occurs.

For example, when freezing water, gas bubbles nucleate at the advancing solidification front (Carte 1961; Maeno 1967; Bari & Hallett 1974; Lipp et al. 1987), as the dissolved gases are rejected by the growing ice crystal, and accumulate at the front, leading to a favourable environment for bubbles to grow. While immiscible, ‘soft’ particles, such as drops, are subjected to stresses (Gerber et al. 2022) during the encapsulation into the ice, leading to potential deformation (Tyagi et al. 2021, 2022; Meijer et al. 2023b; Meijer, Kant & Lohse 2024), the complexity during bubble entrapment is amplified due to additional mass transfer. The shapes and sizes of the entrapped gas bubbles in ice, ranging from very small and barely deformed to elongated vertical cylinders, are set by a delicate interplay between the rate of freezing and the rate of mass transfer, and have been studied extensively (Carte 1961; Maeno 1967; Bari & Hallett 1974; Alley & Fitzpatrick 1999; Wei et al. 2000, 2004; Yoshimura, Inada & Koyama 2008; Chu et al. 2019; Shao et al. 2023; Thiévenaz et al. 2024). The initial growth of the bubbles near an advancing solidification front has, however, received only little attention (Bari & Hallett 1974; Lipp et al. 1987) and forms the main focus of this work. Given the presence of gradients in both gas concentration (due to accumulated gases at the front) and temperature (governing the overall freezing process), in combination with the fact that the bubble is approached by a solid interface, a non-trivial growth of these bubbles is ensured. In addition, complex fluid flow structures might emerge around the growing bubble, further affecting its growth. In this work, we aim to delineate the importance of the different processes involved through a combination of experiments and numerical simulations. We will show that this seemingly simple configuration of a bubble near an ice front underlies truly rich physics, making it a prime example of a physicochemical hydrodynamical system out of equilibrium (Lohse & Zhang 2020).

After introduction of the experimental and numerical procedures in § 2, we continue by addressing the front propagation in § 3. We then turn to the main observation of this work which is the bubble growth near the advancing solidification front (§ 4). More specifically, we address the experimentally observed growth at early stages and compare the results with our numerical simulations. This process occurs simultaneously with the overall freezing process but at much smaller time scales. In § 4.3 we then address the enhanced mass transfer in more detail. We end with conclusions and an outlook in § 5. In the Appendix we also discuss the potential effect of fluid flow around the growing bubble and find it to be small.

2. Experimental procedure and numerical methods

2.1. Experimental procedure

The aim of the experimental set-up is to freeze a sessile water drop on a cold substrate. During the freezing process, bubbles will naturally nucleate and grow at and near the advancing solidification front. To avoid lensing effects we are interested in freezing only a thin slice of purified water (Milli-Q), which we will refer to as ‘drop’ in the remainder of the text. In order to achieve such a quasi-two-dimensional drop, an aluminium mount is placed on top of a freezing stage (BFS-40 MPA, Physitemp) that allows two acrylic plates to be pressed against a thin metal strip (see figure 1). The gap between the plates is ${1}\ {\rm mm}$ and the temperature of the substrate close to the base of the drop, $T_b$, is measured by a thermocouple that is placed inside a groove at the side of the metal strip. We make sure that the desired substrate temperature has been reached well within ${\pm }{0.1}\ {\rm K}$ for several minutes before starting the experiment (see Appendix A). A needle (Nordson) and a syringe pump (PHD 2000 Infusion, Havard Apparatus) are used to deposit drops of equal volume ($V_d = {25}\ {\mathrm {\mu }}{\rm l}$) between the plates, resulting in drop heights of roughly ${3}\ {\rm mm}$. To guarantee freezing as soon as the (room temperature) water touches the substrate and to avoid supercooling, we only deposit the drop once ice crystals have formed on top of the thin metal strip. A gentle flow of nitrogen along the outsides of the plates prevents fog and frost formation that otherwise would obscure the view. The drop is illuminated with a diffused cold LED to avoid local heating. The freezing process is recorded in side view using a camera (Nikon D850) connected to a long working distance lens (Thorlabs, MVL12X12Z). Once the freezing process is complete, the plates are dismounted, the frozen drop is removed and the experimental set-up is cleaned. Experiments are then repeated for different substrate temperatures.

Figure 1. The experimental Hele-Shaw set-up (a) and schematic side views of a freezing sessile drop sandwiched between the two acrylic plates (b). A thermocouple is used to measure $T_b(t)$ and the substrate temperature is fixed at $T_s$. For the analysis we distinguish between front propagation driven by thermal effects (§ 3), and the more rapid (diffusive) bubble growth near the front (§ 4).

Whereas the overall freezing process, driven by thermal effects, and the front propagation are extensively discussed in § 3, the growth of gas bubbles near the front, occurring at much smaller time scales, forms the focus of § 4. Before reporting on all the experimental observations, we briefly go over the technical details of our numerical simulations, for the interested reader, that are used in § 4.2 and Appendix C, and that mimic the bubble growth near the moving front.

2.2. Numerical simulation

For a more detailed analysis of the growth of the bubbles near an advancing solidification front (discussed in § 4.2 and Appendix C), we rely on axisymmetric numerical simulations that account for all the relevant physical mechanisms and allow for a comparison with the experimental observations. To this end, we use a sharp-interface arbitrary Lagrangian–Eulerian finite element method from the pyoomph package (Diddens & Rocha 2024), based on oomph-lib (Heil & Hazel 2006) and GiNaC (Bauer, Frink & Kreckel 2002). The domains here are represented by a triangular mesh.

In our simulations, we consider a rectangular domain with a hole representing the bubble. The bubble has an initial radius of $R_0 = {10}\ {\mathrm {\mu }}{\rm m}$ and is centred at the axis of symmetry, located ${125}\ {\mathrm {\mu }}{\rm m}$ from the bottom boundary. The sidewalls of the domain are far enough from the bubble (25$R_0$) to neglect any boundary effects. We assume the dissolved gases in the water to be a single species (see § 4.2). The gas concentration, i.e. the mass of dissolved gas per volume, denoted by $C$, evolves according to an advection–diffusion equation, see (2.1). If thermal effects are considered, the temperature field $T$ is governed by the heat equation, see (2.2). We neglect thermophoresis effects (Piazza 2008) as our tests indicate that their inclusion is of no significance. The fluid motion in the liquid surrounding the bubble is described by the incompressible Navier–Stokes equations (see (2.3a,b)), leading to a system of equations that reads

(2.1)$$\begin{gather} \partial_t C + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} C = D \nabla^2 C, \end{gather}$$

(2.2)$$\begin{gather}\partial_t T + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \kappa \nabla^2 T, \end{gather}$$

(2.3a,b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0 \quad \text{and} \quad \rho_{l} \left( \partial_t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} \right) =- \boldsymbol{\nabla} p_l + \mu \nabla^2 \boldsymbol{u}, \end{gather}$$

where $\boldsymbol {u}$ is the fluid flow velocity, $D$ and $\kappa$ the gas and thermal diffusivity, respectively, $\rho _l$ is the liquid density, $\mu$ its dynamic viscosity and $p_l$ the pressure in the liquid.

At the interface of the bubble, we impose a constant gas concentration, specifically the saturation concentration $C_{{sat}}$. Additionally, we take the fluid velocity at the interface to be continuous and consider a dynamic boundary condition, including Marangoni stresses

(2.4)$$\begin{gather} p + \mu \boldsymbol{n} \boldsymbol{\cdot}( \boldsymbol{\nabla}\boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^{\rm T}) \boldsymbol{\cdot}\boldsymbol{n}= \sigma \mathcal{K}, \end{gather}$$

(2.5)$$\begin{gather}\mu \boldsymbol{n} \boldsymbol{\cdot} ( \boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^{\rm T})\boldsymbol{\cdot}\boldsymbol{t} = \boldsymbol{\nabla}_{\varGamma} \sigma \boldsymbol{\cdot}\boldsymbol{t}. \end{gather}$$

Here, $p = -p_l + p_g$, with $p_l$ the pressure in the liquid and $p_g$ the gas pressure in the bubble; $\sigma (T)$ is the (temperature-dependent) surface tension, $\mathcal {K}$ is the curvature of the interface, $\boldsymbol {n}$ and $\boldsymbol {t}$ are the outwards-pointing normal and tangent to the bubble–liquid interface, respectively, and $\boldsymbol {\nabla }_{\varGamma }$ is the surface gradient operator.

Despite not modelling the flow inside the bubble, it is crucial to consider the pressure exerted by the bubble on the surrounding fluid. The mass of gas within the bubble, denoted by $m_g=\rho _g V_g$, where $\rho _g$ is the gas density and $V_{g}$ the bubble volume, increases due to the gas mass transfer from the surrounding supersaturated liquid, i.e. locally $C > C_{{sat}}$, into the bubble through diffusion. More specifically,

(2.6)\begin{equation} \partial_t m_g = \int_{\partial \varOmega} -j\,{\rm d}A, \quad \text{with}\ j=\rho_g D\boldsymbol{\nabla} C\boldsymbol{\cdot}\boldsymbol{n}, \end{equation}

where $\partial \varOmega$ is the bubble–liquid interface domain, $j$ is the mass transfer rate and $A$ is the interfacial area of the bubble. The kinematic boundary condition satisfies

(2.7)\begin{equation} \rho_g(\boldsymbol{u} - \partial_t \boldsymbol{R})\boldsymbol{\cdot}\boldsymbol{n}=-j, \end{equation}

with $\partial _t \boldsymbol {R}$ the Lagrangian interface velocity. It is important to note that $\rho _g$ may not remain constant during the bubble growth. We therefore assume it to evolve according to the ideal gas law, $p_g V_g = m_g \bar {R} T_g$, where $\bar {R}$ is the gas constant, and $T_g$ the temperature within the bubble, considered constant.

Lastly, to account for the planar solidification front that approaches the bubble during its growth, we let the bottom boundary advance with a prescribed velocity $V_{{front}}$ (see § 3.2). Its temperature is kept constant and is set to the melting temperature of water $T_m = {0}\,^{\circ }{\rm C}$. We assume no significant action of volume-change convection at the wall (Davis 2001), i.e. $\boldsymbol {n} \boldsymbol{\cdot} \boldsymbol {u}= V_{{front}} (1 - \rho _i/\rho _l ) = 0$, thereby disregarding any density difference between water ($\rho _l$) and ice ($\rho _i$). The top boundary, which is sufficiently distanced from the bubble, moves with the bottom boundary to mimic an infinite domain at each time step. Its temperature and gas concentration are kept constant at $T_{{top}}$ and $C_0$, respectively (see Appendices A and B).

The evolution of the gas pressure, volume and mass is implemented using three Lagrange multipliers, each representing one of these quantities. These are coupled as described in the system of equations and solved monolithically. While the liquid–gas interface is allowed to grow, its centre position is fixed via an additional Lagrange multiplier. We do so, in order to account for the experimentally observed pinning of the bubble (see § 4.2), which likely occurs on the acrylic plate. This imposes a force in the liquid bulk to offset the hydrostatic pressure, ensuring the bubble remains neutrally buoyant. Before solving the system of equations transiently, we first satisfy the stationary diffusive conditions for the gas concentration and temperature fields around the bubble. This is achieved by solving (2.1) and (2.2) with the appropriate mentioned boundary conditions while disregarding the velocity and time-dependent terms, i.e. $\nabla ^2 C = 0$ and $\nabla ^2 T = 0$, respectively. We track the resulting bubble radius at each time step by calculating it from the value of the Lagrange multiplier representing its volume and compare with the experimental results in the subsequent sections.

3. Front propagation

3.1. General experimental observations

When the room temperature drop is deposited onto the cooled, frosted substrate, it immediately starts to freeze from the bottom up. The freezing processes for three different substrate temperatures are shown as sequences of images in figure 2. In all cases, we observe the immediate emergence of a shaded region around the drop that slowly fades away as the drop solidifies, which is known as a ‘frost halo’ (Jung, Tiwari & Poulikakos 2012). Due to the rapid phase change at very early times (shock freezing), an excess amount of released latent heat causes the drop to evaporate. The water vapour consequently condensates at the inside of the cold acrylic plates before slowly evaporating once more. Apart from the significant differences in time it takes the drops to solidify (${203}\ {\rm s}, {95}\ {\rm s}, {35}\ {\rm s}$ for $T_s = {-5}\,^{\circ }{\rm C}, {-10}\,^{\circ }{\rm C}, {-20}\,^{\circ }{\rm C}$, respectively), discussed in more detail in § 3.2, we observe that the base width of the deposited drop, $w_{{b}}$ (see figure 1), becomes smaller when the substrate temperature is reduced. As the drop is deposited from a certain height, which might slightly vary from case to case, the spreading dynamics is governed by an interplay between contact line motion and phase transition phenomena (Koldeweij et al. 2021; Grivet et al. 2022), leading to an earlier arrest of the contact line on colder substrates. During the solidification process the roughly planer ice front advances at a velocity $V_{{front}} = \mathrm {d}y_{{front}}/\mathrm {d}t$ (see figure 1 and § 3.2). As the water solidifies into a crystalline structure, i.e. ice, the gases dissolved in water are expelled and accumulate at the moving front (Tiller et al. 1953). This lowers the threshold for bubbles to nucleate at the front, after which they start to grow. The initial growth of such bubbles is the main focus of the second part of this paper, see § 4. The number of nucleated bubbles and their final shapes and sizes when engulfed in the ice are governed by the rate of freezing, and hence the substrate temperature (see figure 2(a)–(c)), as well as the locally available gas content (Carte 1961; Maeno 1967; Bari & Hallett 1974). Whereas many small bubbles are trapped in the ice when $T_s = {-20}\,^{\circ }{\rm C}$, fewer and more elongated ones are formed when increasing the substrate temperature to $T_s = {-10}\,^{\circ }{\rm C}$ and eventually $T_s = {-5}\,^{\circ }{\rm C}$ (see most right panels in figure 2(a)–(c)). At even lower freezing rates bubbles do not nucleate anymore and clear ice is formed (Bari & Hallett 1974). Finally, at the latest stage of the freezing process, we recover the iconic pointy tip of a frozen water drop (Marin et al. 2014), albeit less pronounced compared with the three-dimensional case due to confinement.

Figure 2. Sequences of images (supplementary movies 1–3 available at https://doi.org/10.1017/jfm.2024.777) highlighting the differences when freezing a sessile drop of equal volume ($V_d = {25}\ {\mathrm {\mu }}{\rm l}$) sandwiched between two acrylic plates at three different substrate temperatures, i.e. (a) $T_s = {-5}\,^{\circ }{\rm C}$, (b) $T_s = {-10}\,^{\circ }{\rm C}$, and (c) $T_s = {-20}\,^{\circ }{\rm C}$ from top to bottom, respectively. All scale bars (right panels) are ${1}\ {\rm mm}$.

3.2. Rate of freezing

Before addressing the bubble growth near the ice front, we first briefly turn to the solidification dynamics of the drop, quantified by the position of the ice front $y_{{front}}(t)$. For several drops at three different substrate temperatures we measure the ice thickness as a function of time and observe repeatably that, as expected, faster freezing occurs at lower substrate temperatures (see figure 3). Given the quasi-two-dimensional nature of the experiments, we model the rate of freezing by simply considering the heat balance at a planar solid–liquid interface in the absence of bulk flow, which is given by Fourier's law of heat conduction (Davis 2001)

(3.1)\begin{equation} \rho_i \mathcal{L}V_{{front}} = \left(k_i \boldsymbol{\nabla} T_i - k_{\ell} \boldsymbol{\nabla} T_{\ell} \right) \boldsymbol{\cdot} \boldsymbol{n}. \end{equation}

Here, $\rho _i$ is the mass density of ice, $\mathcal {L}$ the latent heat of solidification, $k$ the thermal conductivity of ice ($i$) and water ($\ell$), respectively, $T$ the temperature in both phases and $\boldsymbol {n}$ the normal vector to the interface. Assuming that all the heat conducted through the ice contributes to its growth and that heat conduction through the liquid can be neglected, (3.1) simplifies to

(3.2)\begin{equation} \rho_i \mathcal{L}\frac{\mathrm{d}y_{{front}}}{\mathrm{d}t} =k_i \frac{T_m - T_b}{y_{{front}}} = \lambda \left( T_b - T_s \right), \end{equation}

where $T_m$ is the melting temperature of water, $T_b$ the bottom temperature and $T_s$ the imposed temperature on the substrate. We model the heat flux from the substrate to the boundary using a heat-transfer coefficient $\lambda = k_s/d$, where $k_s$ is the thermal conductivity of the (aluminium) substrate and $d$ its (effective) thickness. The second equation in (3.2) is solved to determine $T_b(y_{{front}})$ (see Appendix A), and to extract values for $\lambda$. The first equation in (3.2) then yields

(3.3)\begin{equation} \frac{\mathrm{d}y_{{front}}}{\mathrm{d}t} = \frac{k_i}{\rho_i \mathcal{L}} \frac{1}{y_{{front}}} \left[T_m - \frac{k_i T_m + \lambda T_s y_{{front}}}{k_i + \lambda y_{{front}}} \right], \end{equation}

which is solved numerically and compared with the experimental observations (see red lines in figure 3(a)), resulting in a good agreement when $T_s = {-20}\,^{\circ }{\rm C}, {-10}\,^{\circ }{\rm C}$ and an over-prediction when $T_s = {-5}\,^{\circ }{\rm C}$. We argue that this discrepancy stems from the additional heat loss, potentially through the acrylic plates, that becomes more important at later times during slower freezing, and is not taken into account in the model. As $y_{{front}} \rightarrow \infty$, the right-hand term in the brackets of (3.3) approaches $T_s$ and it follows that $y_{{front}}(t) \approx (2\kappa _i \,St\, t )^{1/2}$, where $\kappa _i = k_i/(\rho _i c_p) \approx 1\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$ is the thermal diffusivity in the ice, and where the Stefan number, $St = \mathcal {L}/c_p(T_m - T_s)$ is defined as the ratio of latent to sensible heat, with $c_p$ the heat capacity of ice. At the later stages of the freezing process we converge to the appropriate scaling of $y_{{front}} \sim t^{1/2}$ (see figure 3(b)). The observed deviations between experiment and theory at the end, especially for the case where $T_s = {-20}\,^{\circ }{\rm C}$, are most likely a consequence of the non-planer geometry of the drop.

Figure 3. Position of the ice front, $y_{{front}}$, as a function of time on linear ($a$) and double-logarithmic ($b$) scale, for repeated experiments (different symbols) at three different substrate temperatures (different colours). The solid red lines in ($a$) correspond to (3.3).

4. Bubble growth

In this section we elaborate on the gas accumulation at the moving ice front and the initial growth of bubbles near the advancing front, which forms the main focus of this paper.

4.1. Gas accumulation at the moving ice front

Gases are naturally dissolved in water. As water turns into ice, these gases are expelled and accumulate at the advancing solidification front. Under steady conditions, at times long compared with the diffusive time $D/(KV_{{front}}^2 )$, where $K$ is the partition coefficient, and for constant front propagation $V_{{front}}$, the expression of the steady concentration profile of the accumulated gases at the front reads (Tiller et al. 1953; Pohl 1954; Carte 1961)

(4.1)\begin{equation} C(y) = C_0 \left( 1 + \frac{1-K}{K} \exp\left[-\frac{V_{{front}}}{D} y \right] \right), \end{equation}

where $D \approx 1\times 10^{-9}\ {\rm m}^2\ {\rm s}^{-1}$ is the gas diffusion coefficient (Bari & Hallett 1974), $K = 0.037 \pm 0.007$ is the partition coefficient for the gases in a water–ice system and $C_0 = (9 \pm 3 ) \ {\rm mg}\ {\rm l}^{-1}$ the gas concentration in water far away from the front, both determined experimentally (see Appendix B).

Nucleated bubbles grow in this enriched environment near the moving ice front at a sufficiently high rate, as we will show, that it allows for the quasi-steady approximation on the front propagation in (4.1) to be used. Additionally, any thermal effects on the bubble growth itself are neglected, as the thermal gradient in the liquid ahead of the ice front is sufficiently small (see Appendix A). Given the typical size of the bubble, the water is therefore assumed to be isothermal at a value close to $T_m$, during its growth. Thermal effects only govern the overall freezing process of the droplet, which occur at a much larger time scale.

4.2. Bubble growth near the approaching ice front

Bubbles nucleate at the moving ice front since the gases dissolved in water are expelled by the growing ice crystal, accumulate at the front and consequently locally lower the nucleation threshold. It is important to note that the content of these bubbles can differ significantly compared with the composition of air in the atmosphere (${\approx }0.79/0.21/4\times 10^{-4}$ vol. for ${\rm N}_{2}/{\rm O}_{2}/{\rm CO}_{2}$). The reason for this is the substantial difference in solubility of the individual gases in water in combination with the different partial pressures of the respective components under standard conditions (Lohse & Hilgenfeldt 1997). Using Henry's law we estimate a bubble composition that is enriched in ${\rm O}_2$ and ${\rm CO}_2$ compared with air (${\approx }0.59/0.38/0.03$ vol. for ${\rm N}_{2}/{\rm O}_{2}/{\rm CO}_{2}$), in line with experimental observations (Bruns 1937; Matsuo & Miyake 1966; Tsurikov 1979). For the remainder of the text we will refer to the content of the bubble as ‘gas’ and assume parameter values for the individual (dominant) components ${\rm N}_2$ and ${\rm O}_2$ (see Appendix C).

Although the majority of bubbles nucleate at the solidification front, where the rough ice surface functions as favourable nucleation site, occasionally, but not too sparsely, we observe bubble nucleation near, but not at, the advancing solidification front (see figure 4(a)). At a certain initial distance from the front a nucleation site might be present in the form of a scratch in the acrylic plate, a dust particle or some other type of impurity, hence assuming some kind of nucleation crevice (Atchley & Prosperetti 1989). The size of this initial crevice is below our optical resolution (see supplementary movie 4). As the ice front approaches, and with that the accumulated gases, this nucleation site experiences a rapidly changing and gas-enriched environment, causing it to grow into a spherical bubble, whose size is characterised by $R_{{bub}}$ (see figure 4(b)). Initially, when $R_{{bub}} \ll \delta _C$, where $\delta _C = D/V_{{front}} \approx 100\ {\mathrm {\mu }}{\rm m}$ is the typical length of the diffusion boundary layer at the front (see (4.1)), it might be assumed that the bubble experiences a homogeneous background concentration field (Bari & Hallett 1974; Lipp et al. 1987). During this stage, the bubble growth is dominated by diffusion and is well described by the steady Epstein–Plesset model (Epstein & Plesset 1950)

(4.2)\begin{equation} \frac{\mathrm{d}R_{{bub}}}{\mathrm{d}t} = \frac{D \Delta C}{\rho_g}\frac{1}{R_{{bub}}}, \quad \text{and thus} \ R_{{bub}}^2 = \frac{2D \Delta C}{\rho_g} \left(t -t_0 \right), \end{equation}

with $\rho _g$ the mass density of the gas and $\Delta C = C(H_0) - C_{{sat}}$ the initial concentration difference between that at the bubble interface and the far field. We define $t-t_0<0$ as the time when $R_{{bub}}$ would have been zero and find, rather surprisingly, that for all cases the bubble-front distance $H$ (see figure 4(d)) linearly decreases according to the front velocity $V_{{front}}$. Consequently, the centres of the bubbles considered in this work remain stationary during their growth and show no migration towards or away from the front, and do not rise due to buoyancy. The reason presumably is some sort of pinning to the sidewalls.

Figure 4. ($a$) Sequence of images (supplementary movie 4) showing the growth of a bubble ahead of the solidification front, as well as ($b$) the bubble radius $R_{{bub}}(t)$ as function of time $t$, for several different bubbles (different colours). ($c$) Bubble growth normalised by the final bubble size $R_{{bub},\infty }$, before it is engulfed into the ice, as compared with (4.2) (dashed black line). We define $t-t_0<0$ as the time when $R_{{bub}}$ would have been zero. The inset shows the initial concentration difference, $\Delta C$, governing the initial growth rate (see (4.2)) as a function of the advancing velocity $V_{{front}}$. ($d$) Schematic of the bubble near the ice moving front in a supersaturated environment. For all cases we find that $H \approx H_0 - V_{{front}}t$, with $H_0$ the initial bubble-front distance, meaning that bubble centres remain stationary during their growth. The typical length of the diffusion boundary layer at the front is denoted by $\delta _C$.

The initial concentration difference, $\Delta C$, depends on (among others) the initial size of the bubble, setting the saturation concentration $C_{{sat}}$ at the interface, the initial bubble-front distance ($H_0 \approx {125}\ {\mathrm {\mu }}{\rm m}$, see Appendix B), and the advancing velocity of the front (see inset figure 4(c)). We obtain $\Delta C$ by fitting the slopes of the experimental data according to (4.2). Rescaling by the size of the bubble as it touches the ice, $R_{{bub,\infty }}$, and the typical diffusion time, $t_{{diff}} = \rho _g R_{{bub},\infty }^2 / (2 D \Delta C)$, shows that the initial bubble growth for all cases indeed follows (4.2) (see figure 4(c)). However, as the bubble grows and the front approaches, its growth deviates from that of pure diffusion under initial conditions and is enhanced.

To pinpoint the origin of the enhanced bubble growth, we perform numerical simulation to study their dynamics under various conditions (see § 2.2 for technical details). We ensure that all physical parameters are known from the literature or obtained experimentally, and that we do no rely on any adjustable fitting parameters. Initially, advection is omitted and we only consider growth through diffusion. The effect of advection around the bubble is addressed in more detail in Appendix C.

For the numerical simulations we consider three distinct cases. Firstly, we take the case of a bubble growing in a homogeneous background concentration, with the front, located at a distance $H_0 = {125}\ {\mathrm {\mu }}{\rm m}$, remaining static. We retrieve a behaviour very similar to (4.2) (see black line figure 5(a)), where the deviation at the end is a consequence of the presence of the static interface, bounding the domain, and slightly altering the concentration field around the bubble (see panel I in figure 5(b)). Secondly, we assume a background concentration gradient according to (4.1), resulting in a significantly faster growth of the bubble (see blue line and panel II in figure 5(b)). Thirdly, we consider the front to be in motion with velocity $V_{{front}}$ and to approach the growing bubble, which further enhances its growth (see green line and panel III in figure 5(b)). The determined bubble growth of the latter agrees very nicely with our experimental observations, where the faster growth originates from a combination of the background gas concentration gradient and the motion of the interface, causing an accumulation of the iso-concentration lines at the bottom of the bubble (see panel III in figure 5(b)). This alters the local concentration gradient close to the interface of the bubble, leading to an enhanced mass transfer and thus a more rapid growth.

Figure 5. (a) Comparison between experiments (data points) and numerical simulations (solid lines) (supplementary movie 5) on the growth of a bubble near a solidification front, namely when considering (b) no advection for (I) a homogeneous background concentration with a static front, and for a background concentration gradient according to (4.1) with a static (II) or moving (III) front. The panels show the concentration profiles around the bubbles 5 s after initialisation. The initial bubble radius is $R_0 = {10}\ {\mathrm {\mu }}{\rm m}$, $C_{{sat}} = {0.029}\ {\rm kg}\ {\rm m}^{-3}$ and $V_{{front}} = {10}\ {\mathrm {\mu }}{\rm m}\ {\rm s}^{-1}$ (for III only).

4.3. Enhanced mass transfer

To further rationalise the observed enhanced mass transfer, we recast our problem of a bubble growing near an advancing solidification front into that of the well-studied case of mass (or equivalently heat) transfer to a moving spherical object in an infinite domain with a constant background concentration (or temperature). The dimensionless diffusional flux to the moving object can then be expressed as (Levich 1962)

(4.3)\begin{equation} Sh \propto \left( Re \, Sc \right)^{1/2}, \end{equation}

where the dimensionless Sherwood, Reynolds and Schmidt numbers are defined as

(4.4a–c)\begin{equation} Sh = \frac{2R_{{bub}} \dot{R}_{{bub}} \rho_g}{D \Delta C}, \quad Re = \frac{\rho U R}{\mu}, \quad \text{and} \quad Sc = \frac{\mu}{\rho D}, \end{equation}

with $\dot {R}_{{bub}} = \mathrm {d}R_{{bub}}/\mathrm {d}t$ the growth rate of the bubble and $U$ the relative fluid velocity at the equator of the spherical object. Whereas the Sherwood number is defined as the ratio of the convective mass transfer rate to diffusive mass transport, the product of Reynolds and Schmidt numbers is equivalent to a Péclet number, i.e. $Re \, Sc = Pe$, which is defined as the ratio of the rate of advection to that of diffusion. In the limit that $(Re \, Sc ) \rightarrow 0$ the solution to this classical case is $Sh = 2$, where the concentration field looks similar to that depicted in panel I of figure 5(b). As the object moves through the (homogeneous) concentration field, the iso-concentration lines tend to accumulate at one end (Steinberger & Treybal 1960; Frossling 1963; Ihme, Schmidt-Traub & Brauer 1972; Oellrich, Schmidt-Traub & Brauer 1973), similar to panel III of figure 5(b), leading to an enhanced mass transfer according to (4.3).

In our case, the bubble is stationary and the background concentration is not constant. To account for the moving background concentration gradient the bubble experiences, we introduce an effective velocity as

(4.5)\begin{equation} U_{{eff}} = \left\vert \frac{\mathrm{d} C}{\mathrm{d}y} \right\vert \frac{R }{C_0} V_{{front}}. \end{equation}

Considering that $U = U_{{eff}}$, we obtain that our experimental observations are in line with (4.3) as we approach the limit $(Re \, Sc ) \rightarrow \infty$ (see dashed black line in figure 6). In the opposite limit, when the bubble and the background concentration gradient are both small, i.e. $(Re \, Sc ) \rightarrow 0$, the system undergoes a transition towards a different scaling, approaching $Sh \approx 2$, where the notable scatter arises from inaccuracies in determining the bubble growth rate $\dot {R}_{{bub}}$ for smaller bubbles.

Figure 6. Value of $Sh$ vs $Re\,Sc$ for different bubbles growing near the advancing solidification front, represented by different colours. The data tend towards the limits $Sh = 2$ and $Sh \sim (Re \, Sc )^{1/2}$ (dashed black lines) (Levich 1962).

5. Conclusions and outlook

To conclude, we have investigated the freezing of a quasi-two-dimensional sessile drop on a cooled substrate. During the freezing process, gases dissolved in water are rejected by the growing ice crystal and accumulate at the moving ice front, creating a supersaturated region which favours the nucleation of gas bubbles. These bubbles then grow, before eventually being engulfed into the ice. Their final shape and size in the ice are governed by a delicate interplay between the rate of freezing and the rate of mass transfer towards the bubble (Bari & Hallett 1974; Lipp et al. 1987).

In this work, we focused on the experimental and numerical investigation of the initial growth of those bubbles that nucleate not at but near the solidification front. We show that at the early stages the growth is dominated by diffusion and is enhanced due to a combination of the presence of the background gas concentration gradient and the motion of the approaching front. We rationalised that, based on the numerical simulations, which agree nicely with the experimental observations, the iso-concentration lines tend to accumulate at the bottom of the bubble, leading to an enhanced mass transfer across the interface of the bubble, and hence a faster growth. We have additionally shown that our problem can be recast into that of mass transfer to a moving spherical object in a homogeneous concentration field, finding good agreement between our experimental data and the existing scaling relations. Finally, through numerical simulations we have qualitatively addressed how fluid flow around the bubble might further affect its growth (see Appendix C).

Our findings shed new light on the diverse processes that might govern the growth of bubbles near a moving interface, subjected to gradients in both concentration and temperature. Besides any solidification process of pure or multi-component liquids, such as metals, silica or sapphires, where these findings might be relevant, they might also bring new insight into the formation of gas bubbles in hailstones (Bari & Hallett 1974) and lake ice (Swinzow 1966; Gow & Langston 1977) under various conditions. As continuation, it might be interesting to investigate how the bubble growth is altered when freezing liquids that are saturated with a more (less) soluble gas, such as ${\rm CO}_2$ (Ar), to further explore the relevance of our findings.