2. Theoretical development

We begin by developing a general model for melting or dissolving bodies under stable natural convection. Melting and dissolving bodies can be modelled using a system of equations that are fundamentally the same, but which differ in underlying physical processes. For dissolution (e.g. candy in water), the boundary recedes in response to mass loss through molecular diffusion into the solvent. For melting (e.g. ice in water), the material instead undergoes a phase change induced by thermal conduction. In either case, the recession speed of the boundary is proportional to the net rate of molecular or thermal diffusion transported from or to the interface, respectively (the Stefan condition in the thermal context). In the case of a dissolving solid, diffusion is typically negligible within the interior of the solid itself, and hence the recession speed of the boundary is proportional to the flux of solute into the ambient fluid. In the context of a melting body, thermal conduction can, in principle, be important within the body. In this case, a melting body can either first equilibrate to the melt temperature, or proceed to melt with a thin conductive boundary layer under the interface; either of these cases produces an effective proportionality between exterior heat flux and recession rate (the conditions producing these distinct situations will be clarified as part of the forthcoming paper; Pegler & Davies Wykes Reference Pegler and Davies Wykes2021). With these caveats noted, we henceforth adopt the notation and terminology applicable for dissolving bodies.

We consider first the case of two-dimensional bodies. Let $x$ and $z$ denote the horizontal and vertical coordinates in the laboratory frame, respectively. Suppose that the solid has density $\rho _s$ with surface $z = h(x,t)$ and is immersed in a quiescent, incompressible fluid of density $\rho _{\infty }$ in its pure form, as illustrated in figure 1. The dissolution of the solid forms a solutal boundary layer along the exterior of the solid driven by buoyancy (natural convection). Let $\rho (x,y,t)$ denote the density of the fluid, and let $C(x,z,t)$ denote the normalised concentration of solute in the fluid, defined by

(2.1)\begin{equation} C(x,z,t) \equiv \frac{\rho(x,y,t) - \rho_{\infty}}{\rho_s - \rho_{\infty}}. \end{equation}

The value $C=1$ represents pure solid ($\rho = \rho _s$) and $C=0$ represents pure ambient fluid ($\rho = \rho _{\infty }$). At the interface between the solid and the fluid, we impose

(2.2)\begin{equation} C(x,h,t) \equiv C_{i} = C_{sat}, \end{equation}

where $C_{i}$ is the interfacial concentration, which is assumed to equal the concentration of the solution at saturation, denoted $C_{sat}$. For an aqueous sucrose solution at $25\, ^\circ \textrm {C}$, the saturation concentration is $C_{sat} \approx 0.68$ (Mohos Reference Mohos2010). The value of $C_{sat}$ is less than unity because solvation requires a minimum concentration of water. The saturation concentration $C_{sat}$ depends on temperature, and thus may change if the latent heat generated by dissolution is sufficiently large (thereby requiring simultaneous solutions for both the temperature and concentration fields, rather than only the concentration). By deriving a theoretical prediction for the temperature generated by this effect (appendix A), we estimate the temperature increase to be at most ${\approx }0.8\ \textrm {K}$ for candy dissolving in water, resulting in a negligible change in $C_{sat}$ of ${<}0.2\,\%$. A discussion of the effect of interface temperature on interface concentration is provided in Wells & Worster (Reference Wells and Worster2011). Such two-component effects can in fact be incorporated directly into our model for shape evolution developed below. Nonetheless, we will neglect temperature effects for now for the purpose of simplifying our exposition.

Figure 1. $(a)$ Schematic illustrating the general configuration and notation for a body dissolving under gravitationally stable natural convection. Two coordinate systems are illustrated: ($x,z$) is fixed, while ($s,y$) represents a curvilinear coordinate system in which $s$ is the arc length of the solid surface from the tip and $y$ is the normal distance from the boundary. Panels $(b)$ and $(c)$ illustrate two different limits of the general model (2.16a,b): the shallow limit, where the slope $h_x \ll 1$, and the steep limit, where $h_x \gg 1$, respectively.

The shape evolution is controlled by spatial and temporal variations in the rate of dissolution within the solutal boundary layer. To describe the boundary-layer flow, we use the curvilinear coordinate system $(s, y)$, where $s$ is the arc length along the surface $z=h(x,t)$ measured from the tip $h_0(t) \equiv h(0,t)$, and $y$ is the perpendicular outward distance from the surface, as shown in figure 1$(a)$. The flux of concentration (per unit width) perpendicular to the solid boundary is, in accordance with Fick's law,

(2.3)\begin{equation} q(s,t) ={-}\kappa_i \left. \frac{\partial C}{\partial y}\right|_{y=0}, \end{equation}

where $\kappa _i$ is the diffusivity of the fully saturated solution at the boundary. Let $V(s,t)$ denote the unknown speed at which the surface recedes in the direction of its normal. The speed of recession can be related to the flux (2.3) via the solute-conservation equation

(2.4)\begin{equation} V(s,t) = R q(s,t) ={-}R \kappa_i\left. \frac{\partial C}{\partial y}\right|_{y=0}, \end{equation}

where $R \equiv \rho _i / [\rho _s (1 - C_i)]$ is a dimensionless conversion factor between concentration flux and boundary migration speed, and we have assumed here a linear relationship between mass and volume fraction. Here, we have also neglected flux in concentration in the solid interior, which could be more important in the context of melting bodies. Because the surface moves inwards, we write the vectorial form of the surface velocity as

(2.5a)\begin{equation} \boldsymbol{V} ={-}V\boldsymbol{n}, \end{equation}

where

(2.5b)\begin{equation} \boldsymbol{n} = \frac{\pmb{e}_z - h_x \pmb{e}_x}{ \sqrt{1 + h_x^2}} \end{equation}

is the outward normal to the surface, $\pmb{e}_x$ and $\pmb{e}_z$ are the unit vectors in the horizontal and vertical directions and the $x$ subscript in $h_x$ denotes partial differentiation, $\partial h / \partial x$. Since the solid migrates with the velocity $\boldsymbol {V}$, it follows that

(2.6)\begin{equation} \left( \frac{\partial }{\partial t} + \boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \left[ z - h(x,t) \right] = 0. \end{equation}

On simplification and use of (2.5a,b), we obtain the interface evolution equation

(2.7)\begin{equation} \frac{\partial h}{\partial t} = \boldsymbol{V} \boldsymbol{\cdot} \left( \pmb{e}_z - \pmb{e}_x \frac{\partial h}{\partial x} \right) ={-}\frac{\partial s}{\partial x}V, \end{equation}

where $s_x = \sqrt {1 + h_x^2}$. Substituting (2.4) into (2.7), we obtain the evolution equation

(2.8)\begin{equation} \frac{\partial h}{\partial t} ={-}Rq \frac{\partial s}{\partial x} = R \kappa_i \frac{\partial s}{\partial x} \left.\frac{\partial C}{\partial y}\right|_{y=0}, \end{equation}

providing an equation for the rate of vertical descent of the solid boundary in terms of the dissolution pattern. Thus, to develop a closed model, it remains to relate the normal gradient $\partial C / \partial y$ along the surface to the shape profile $h(x,t)$.

Let $u$ and $v$ define the components of fluid velocity in the $s$ (tangential) and $y$ (normal) directions, respectively. In modelling the boundary-layer flow, we assume that the curvature of the interface is sufficiently small and that the flow is quasi-steady on the timescale of evolution of the free boundary (the criteria for these assumptions will be detailed in § 2.3). The equations describing steady free-convective laminar boundary layers read (Ostrach Reference Ostrach1953)

(2.9)\begin{align} \frac{\partial u}{\partial s} + \frac{\partial v}{\partial y} &= 0, \end{align}

(2.10)\begin{align} \rho \left( u \frac{\partial u}{\partial s} + v \frac{\partial u}{\partial y} \right) &= \frac{\partial }{\partial y}\left( \mu \frac{\partial u}{\partial y} \right) + {\rm \Delta} \rho_s g \sin[\alpha(s,t)] C, \end{align}

(2.11)\begin{align} u \frac{\partial C}{\partial s} + v \frac{\partial C}{\partial y} &= \frac{\partial }{\partial y} \left( \kappa \frac{\partial C}{\partial y} \right), \end{align}

where $\alpha (s,t)$ is the angle of the slope to the horizontal (the angle $\alpha$ is shown red in figure 1), and ${\rm \Delta} \rho _s \equiv \rho _s - \rho _{\infty }$. Equations (2.9)–(2.11) represent the conservation of total fluid mass, momentum and solute mass, respectively.

We allow here for general dependences of viscosity $\mu (C)$, diffusivity $\kappa (C)$ and density $\rho (C)$ on the concentration $C$. For many fluids, including rock candy dissolving in water, the viscosity of the solution $\mu (C)$ can vary by many orders of magnitude (for a solution of rock candy, $\mu$ will vary by at least two orders of magnitude; Quintas et al. Reference Quintas, Brandão, Silva and Cunha2006). Since viscosity and molecular diffusion are related approximately inversely in accordance with the Stokes–Einstein relation, $\kappa \propto \mu ^{-1}$ (e.g. Mohos Reference Mohos2010; Price, Mattsson & Murray Reference Price, Mattsson and Murray2016), order-of-magnitude variations in molecular diffusivity $\kappa (C)$ will also apply. Fortunately, these complications can all be encapsulated concisely within a single parameter.

The boundary-layer system (2.9)–(2.11) is fifth order and hence requires five boundary conditions. One is provided by the Dirichlet condition on the interface concentration (2.2). Further conditions are provided by the no-slip condition and normal velocity condition at the migrating solid boundary given by

(2.12a,b)\begin{equation} u = 0, \quad v = (\rho_s/\rho_i)V, \quad \text{on}\ y = 0, \end{equation}

respectively (e.g. Wells & Worster Reference Wells and Worster2011). The normal velocity condition (2.12b) represents the sum of two contributions: $([(\rho _s - \rho _i) / \rho _i]V + V) = (\rho _s/\rho _i)V$. The first is the speed accounting for the net sourcing of fluid caused by the density difference between the solid and the adjacent solution. The second contribution accounts for the fact that $y=0$ is moving inwards at velocity $V$, and thus $V$ must be added to account for this frame change. Finally, we impose

(2.13a,b)\begin{equation} C \to 0, \quad u \to 0 \quad \text{as } y \to \infty, \end{equation}

completing the set of conditions needed to close the boundary-layer system.

For any given surface profile $h(x,t)$, the system (2.9)–(2.13) forms a parabolic boundary-value problem for $C(s,y,t)$. Once $C$ is known, the normal gradient $\partial C / \partial y$ on $y=0$ can be evaluated and the surface descent rate $\partial h / \partial t$ determined using (2.8). The surface can then be evolved in time and the process repeated. The coupled model (2.8)–(2.13) generalises models of natural convection along dissolving or melting boundaries (Schlichting & Gersten Reference Schlichting and Gersten2017) to incorporate the two-way coupling between boundary-layer dynamics and the evolution of the boundary. Conversely, the model forms a variation of classical (purely diffusive) Stefan problems in which the migration of a melting or solidifying boundary is modelled using conductive transport processes (Stefan Reference Stefan1891; Jaeger & Carslaw Reference Jaeger and Carslaw1959) to assume instead that the shape evolution is driven by the buoyancy-driven convection along the exterior of the body.

2.1. Analytical reduction for large Schmidt number

For general situations in which all terms in the momentum equation (2.10) contribute, one could consider the boundary-layer system (2.9)–(2.13) alongside (2.8). However, we note that if the Schmidt number $Sc \equiv \mu /\rho \kappa$ is sufficiently large, a major simplification of the boundary-layer system (2.9)–(2.13) is available that allows a single integro-differential hyperbolic equation describing the evolution of $h(x,t)$ to be developed. We consider this situation first as an inroad for clearly demonstrating many of the simplified scalings, properties and self-similar forms underlying dissolving shapes. The effects of reinstating inertia are discussed in § 6.

Convective boundary layers described by (2.9)–(2.13) will generally involve both a viscous sublayer, in which viscous stresses become important, and a diffusive sublayer in which the transfer of solute is localised. The relative size of the viscous sublayer compared to the solutal layer is characterised by the Schmidt number $Sc \equiv \nu /\kappa$. If $Sc \gg 1$, the viscous sublayer is considerably larger than the solutal layer, such that the solutal layer forms an inner sublayer inside the viscous sublayer. Exterior to the viscous sublayer lies an inertial outer region in which the transition to the far-field stagnancy condition (2.13) occurs. More specifically, for the purpose of determining the solutal flux along the boundary, it is sufficient to consider just the viscous sublayer subject to a matching condition of zero stress on its far field, a system we review in appendix C. As shown by Kuiken (Reference Kuiken1968), the neglect of inertia in the solutal region introduces an error of ${<}20\,\%$ in the predicted flux for all $Sc > 1$, implying that even order-unity values of $Sc$ are well approximated by the $Sc \to \infty$ limit. Since the Schmidt number is comprised purely of material parameters, and is independent of the shape of the body, for example, it is straightforward to evaluate for a given configuration. For example, $Sc \approx 2\times 10^3$ (dilute) to $8 \times 10^6$ (saturated) using the kinematic viscosity and diffusivity for dilute and saturated solutions of sucrose in water, respectively, and hence $Sc \gg 1$ is firmly satisfied in this context.

The analytical reduction that arises for $Sc \to \infty$ is based on a transformation of the boundary-layer subsystem describing the viscous sublayer that eliminates the interface slope (Acrivos Reference Acrivos1960a), leaving a transformed boundary-layer subsystem that is independent of the surface profile (the transformed system is equivalent to the similarity system describing free convection along a vertical wall, for example). By considering the similarity solution to the transformed system and then recasting the system back in terms of the original variables, one obtains the dissolution profile along an arbitrary two-dimensional shape as

(2.14)\begin{equation} q(s,t) ={-}\kappa_i \left.\frac{\partial C}{\partial y} \right|_{y=0} = \frac{ D \sin[\alpha(s,t)]^{1/3} }{\displaystyle \left( \int_0^s \sin[\alpha(s,t)]^{1/3} \,{\textrm{d}} s \right)^{1/4} }, \end{equation}

where $D \equiv \gamma \kappa _i C_i ({\rm \Delta} \rho _s g C_i/\kappa _{\infty } \mu _{\infty })^{1/4}$ is a constant material parameter, $\mu _{\infty }$ and $\kappa _{\infty }$ are the viscosity and diffusivity of the solvent and $\gamma$ is a dimensionless number. The result above is derived by Acrivos (Reference Acrivos1960a), and a step-by-step review of the derivation is provided in appendix C. As discussed in the appendix, the dimensionless prefactor $\gamma$ depends on the material properties of the fluid and solid under consideration, including any dependences of $\mu$, $\kappa$ and $\rho$ on the concentration $C$, as well as any two-component temperature effects. The presence of $\sin \alpha$ in the numerator of (2.14) implies that a larger along-slope component of buoyancy drives faster convective transport locally and hence more efficient dissolution. The integral in the denominator is a ‘memory term’ representing the total dissolution between the tip and the position $s$. The increase of this integral with $s$ represents the insulating effect of solute accumulation, which reduces the dissolution rate as the boundary layer grows downstream.

Using (2.14), the recession rate of the boundary can be determined as

(2.15)\begin{equation} V(s,t) = Rq = \frac{ A \sin[\alpha(s,t)]^{1/3} }{\displaystyle \left(\int_0^s \sin[\alpha(s,t)]^{1/3} \,{\textrm{d}} s \right)^{1/4} }, \end{equation}

where $A \equiv RD \equiv \gamma R \kappa _i C_i ({\rm \Delta} \rho _s g C_i/\kappa _{\infty } \mu _{\infty })^{1/4}$. The parameter $A$ represents all effects on dissolution associated with the material properties of the fluid and the solid, and is referred to herein as the convective strength. This is distinct from the effects of shape, which are controlled by the terms in expression (2.15) depending on $\alpha (s,t)$. In particular, the number $A$ encapsulates, via $\gamma$, the full functional dependences of viscosity $\mu$, diffusivity $\kappa$ and density $\rho$ on the concentration $C$. The parameter $A$ could, in principle, be determined by measuring the functions of $\mu (C)$, $\kappa (C)$ and $\rho (C)$ for the fluid of interest and solving the similarity system (C 5)–(C 6) for $\gamma$ and, in turn, $A$. Alternatively, one could measure $A$ empirically for a given system of interest using a ‘one off’ observation of an instantaneous recession rate, circumventing the need to measure $\mu (C)$, $\kappa (C)$ and $\rho (C)$ explicitly and solving (C 5)–(C 6). This direct empirical method of inferring $A$ will be discussed further in § 5. For a vertical wall (for which $\alpha = {\rm \pi}/2$ and $s = -z$), (2.15) reduces to $V = A (-z)^{-1/4}$, giving the $(-z)^{-1/4}$ dissolution profile along a wall of constant inclination (e.g. Ostrach Reference Ostrach1953). The convective strength $A$ can thus be interpreted as the prefactor to $(-z)^{-1/4}$ that would apply in an expression for the instantaneous recession speed along a vertical wall for the fluid–solid configuration under consideration. In this case, the recession speed at the top of the vertical wall ($z=0$) is, at least instantaneously, infinite. However, once the evolution of the boundary is also incorporated into the model (see below), the shape at $z=0$ would be instantaneously smoothed in such a way as to produce a finite recession speed for all subsequent times.

Using (2.15), we couple the recession speed $V(s,t)$ fully to the profile $h(x,t)$ via (2.7). Substituting (2.15) into the solute conservation equation (2.7) and noting that $\sin \alpha = h_x/s_x$ and ${\textrm {d}} s = s_x \,{\textrm {d}} x$, we obtain the shape evolution equation

(2.16a,b)\begin{equation} h_t = {\displaystyle \frac{\displaystyle -A \left|h_x s_x^2 \right|^{1/3} }{\displaystyle \left( \int_0^x \left|h_x s_x^2 \right|^{1/3} \,{\textrm{d}} x \right)^{1/4} }}, \quad s_x = \sqrt{1 + h_x^2}. \end{equation}

This result yields a new partial integro-differential equation that independently describes the evolution of a dissolving solid interface $h(x,t)$. The equation shows that the evolution of a dissolving or melting boundary in free convection is governed by a hyperbolic system relating the descent rate $h_t$ to the slope $h_x$. The numerator represents a dependence of the descent rate on the local value of the slope $h_x$, showing that the dissolution rate increases with the local slope. The integral in the denominator represents a global dependence of the descent rate at a position $x$ on all the values of $h_x$ between the tip and the position $x$. This contribution represents the insulating effect of concentration accumulation between the tip and the position $x$.

2.2. Tip shape and its relation to instantaneous descent rate

Before solving (2.16a,b), we first note several implied properties of the local tip shape and its relationship to the instantaneous rate of tip descent. To seek the admissible tip shapes, we try a power-law $h \propto x^{a}$, where $a>1$ is a constant, in (2.16a,b) and take the near-tip limit $x \to 0$. We find that (2.16a,b) then simplifies to

(2.17)\begin{equation} h_t \propto x^{(a-2)/4} \quad (x \to 0). \end{equation}

The rate of tip descent $h_t$ is therefore finite if and only if the tip is parabolic, $a=2$. Otherwise, the tip would move infinitely quickly ($a<2$) or remain stationary ($a>2$). The tip will therefore assume a parabola as it descends at a finite, non-zero speed. This provides a theoretical answer to the question posed by Nakouzi et al. (Reference Nakouzi, Goldstein and Steinbock2015) on whether the shapes of buoyancy-driven descending tips converge to a universal shape, and explains the emergence of a paraboloidal shape indicated by their experimental observations. It should be noted, however, that the parabola only applies close to the tip. The precise asymptotic meaning of how close, and an identification of the different shapes that arise exterior to the parabolic region, will be determined in our later analyses. Indeed, the parabolic region can be so small in the case of steep bodies that an entirely different near-tip shape is more representative. Previous theoretical work on forced convection has found that the local shape of a tip dissolving, melting or eroding in a fluid flow is that which produces a boundary layer with a locally smooth recession rate (Ristroph et al. Reference Ristroph, Moore, Childress, Shelley and Zhang2012; Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013; Huang et al. Reference Huang, Moore and Ristroph2015; Moore Reference Moore2017). We note that the parabola is the unique shape that satisfies the property of a locally uniform recession rate under stable natural convection in the limit $h_x \to 0$ (this property is shown by the expression for $q$ in (C 9a--c) given in appendix C).

The implied parabolic profile near the tip can be written

(2.18)\begin{equation} h \sim h_0(t) - {\frac{1}{2}} S(t) x^2 \quad (x \to 0), \end{equation}

where $S(t) = |h_{xx}(0,t)|$ is the tip curvature or sharpness. Substituting the implied slope $h_x = -S(t)x$ into (2.16a,b), we determine the general law relating tip speed to tip sharpness,

(2.19)\begin{equation} \frac{{\textrm{d}} h_0}{{\textrm{d}} t} ={-}A\left( \frac{4}{3} S(t) \right)^{1/4}. \end{equation}

Sharper tips thus descend faster, at a rate directly proportional to the convective strength $A$ and the $1/4$-power of the tip curvature $S(t)$.

2.3. Assumptions underlying the boundary-layer model

A foundational assumption of the boundary-layer model of (2.9)–(2.11) is that the curvature of the interface is sufficiently small and the boundary-layer sufficiently thin that $\delta h_{ss} \ll 1$, where $\delta$ is the boundary-layer thickness. In this limit, the derivatives with respect to $y$, i.e. the normal shear stress and normal diffusion, yield the dominant contributions to the diffusive and viscous stress terms, as assumed in (2.10) and (2.11). This approximation is strongly satisfied in our laboratory experiments, for example, where the boundary-layer thickness in the laminar regime is small enough to not be visible (figure 14) and is estimated to be ${<}1\ \textrm {mm}$. A further underlying assumption is that the boundary layer remains laminar and attached to the surface. This aspect will be discussed further in § 6.

The system (2.9)–(2.11) is underlain by two further assumptions relating to assumed quasi-steadiness of the flow and the form of the buoyancy force. The first of these arises from the neglect of the time derivatives in (2.10) and (2.11) on the basis that the speed of recession of the surface is much slower than the adjustment time scale of the boundary layer. As an example, the typical speed of the interface in our experiments is ${\sim }1\ \textrm {cm}\,\textrm {h}^{-1}$, while the boundary-layer flow rate is ${\sim }1\ \textrm {cm}\,\textrm {s}^{-1}$, indicating that quasi-steadiness is strongly satisfied. To provide an explicit non-dimensional index to assess the satisfaction of quasi-steadiness and the relative importance of ambient convection during the dissolution of the body, we define a dimensionless number

(2.20)\begin{equation} \Gamma \equiv V/u, \end{equation}

representing the ratio of the recession speed to the rate of convective transport. If $\Gamma \ll 1$ across the majority of the shape, the boundary layer will adjust effectively instantaneously to the solution of the steady boundary-layer equations (2.9)–(2.11). If instead $\Gamma = O(1)$, the boundary-layer flow will evolve on the same time scale as the interface. In this case, a generalisation of the boundary-layer equations (2.9)–(2.11) to include time dependence could be considered. If $\Gamma \gg 1$, the interface will recede considerably faster than the boundary-layer flow. In this limit, convection is negligible and the melt rate is controlled by time-dependent diffusion in the manner of a classical (purely diffusive) Stefan problem. Values of $\Gamma \ll 1$ apply for sufficiently large bodies. To show this, we consider the value of $\Gamma$ for two illustrative examples. For these estimates, we use the scalings for the boundary-layer velocity $u$ given in § C.1. For the case of a parabolic shape, (2.19) yields the scaling $V \sim AS^{1/4}$ and (C 9a--c) yields $u \sim \kappa _{\infty } (SG_*)^{1/2} x$, where $G_* = {\rm \Delta} \rho _s g C_i / \kappa _{\infty } \mu _{\infty }$. Substituting these expressions into (2.20), we obtain $\Gamma \sim l_1/x$, where $l_1 \equiv (R^4 / S G_*)^{1/4}$ is an intrinsic length scale. Thus, there is a location $l_1$ such that the recession speed is faster than the flow for $x \lesssim l_1$ and the recession speed is slower than the flow for $x \gtrsim l_1$. Therefore, if $l/l_1 \gg 1$ where $l$ is the horizontal length scale of the body, the shape evolution is predominantly controlled by the quasi-steady equations. Our laboratory experiments presented in § 5 form a steep body with approximately linear surfaces. To evaluate $\Gamma$ for this case, we note that (2.15) yields the scaling $V \sim A|z|^{-1/4}$ and (C 10a--c) yields $u \sim \kappa _{\infty } (G_* |z|)^{1/2}$, such that $\Gamma \sim (l_2/|z|)^{3/4}$, where $l_2 \equiv (R^4 / G_* )^{1/3} = (A^2 / \kappa _i^2 G_* )^{2/3}$ is a length scale. This scale is ${<}1\ \mathrm {\mu }\textrm {m}$ for our experiments, indicating that quasi-steadiness is strongly satisfied throughout the majority of the flow.

The form of the buoyancy term on the right-hand side of (2.10) is standard in the modelling of free convection along vertical and inclined walls (Schlichting & Gersten Reference Schlichting and Gersten2017). However, we wish to allow for shallow slopes, $h_x \ll 1$, for which the buoyancy force can also be affected by horizontal gradients in hydrostatic pressure caused by the stratified density field (Stewartson Reference Stewartson1958). By considering a scaling analysis of a generalised system of boundary-layer equations detailed in appendix B, we discuss the transition of the flow from a gradient-driven region near the tip, in which gradients in hydrostatic pressure contribute to driving the flow, to a slope-driven region wherein these gradients are negligible and the flow is driven dominantly by the along-slope component of buoyancy, as assumed in (2.10). Assuming a surface boundary of the general power-law form $h \sim -c_nx^n$, where $c_n$ and $n$ are constants, the transition between the gradient and slope driven regimes is determined to occur over the length scale

(2.21a)\begin{equation} \mathcal{L}\equiv \left( \mathcal{L}_*^3 c_n^{{-}5} \right)^{{1}/({5n-2})}, \end{equation}

where

(2.21b)\begin{equation} \mathcal{L}_* \equiv \left( \frac{\mu \kappa}{{\rm \Delta} \rho_s g} \right)^{1/3}. \end{equation}

The case $n=1$ of the above was derived by Stewartson (Reference Stewartson1958), who also discussed its interpretation as the scale on which a gradient-driven region applies for $x \ll \mathcal {L}$. Assuming a representative linear surface, $c_1 = m$, we estimate that $\mathcal {L} = (2 \times 10^{-5}\ m^{-5/3})$ metres using values representative of our laboratory experiments. The length is ${\approx }1\ \mathrm {\mu }\textrm {m}$ for our shallowest laboratory experiment ($m \approx 5.8$). The gradient-driven region is therefore indicated to lie very close to the tip. In order for the gradient-driven region to be appreciable at the laboratory scale (taking $\mathcal {L} = 1\ \textrm {cm}$ for instance), the slope must be $m \lesssim 2 \times 10^{-3}$. The gradient-driven regime (Stewartson Reference Stewartson1958; Neufeld, Goldstein & Worster Reference Neufeld, Goldstein and Worster2010) therefore applies only briefly, even for a shallow slope. On this basis, we anticipate that the along-slope form of buoyancy represented by the model of (2.10) is likely to represent the most widely applicable regime describing dissolving and melting bodies for appreciable shape changes at the laboratory scale and larger.

3. Mathematical analysis

Equation (2.16a,b) describes the evolution of dissolving two-dimensional bodies, subject to specification of any given initial shape $h(x,0)$ with a single maximum. Across the scope of this paper and the companion, we will allow for general initial shapes of the form $h(x,0) \propto -|x|^{n}$, where $n>0$. This form encompasses initially rectangular cross-sectioned surfaces ($n=\infty$), smooth parabolas ($n=2$), wedges ($n=1$) and sharp concave surfaces (e.g. $n=1/2$). General values of $n$ will be found to exhibit a variety of asymptotic transitions between different regimes, owing to the existence of an explicit time scale in the problem $t \sim AL^{5/4}$ representing the time scale for the tip to descend a distance comparable to the initial size of the body $L$ (Pegler & Davies Wykes Reference Pegler and Davies Wykes2021). In the unique case $n=1$, no such time-independent length scale exists (to be shown in § 3.1 below), resulting in similarity solutions that apply for all time. As a stepping stone, we begin here by analysing these similarity solutions, reserving a consideration of general $n$ for the forthcoming paper. Thus, we specify the initial condition

(3.1)\begin{equation} h(x,0) ={-}m|x|, \end{equation}

where $m$ is the initial magnitude of the slopes of the two adjoining planes, forming the unique dimensionless parameter in the problem.

As a preliminary illustration of the model solutions, we present a selection of direct numerical solutions of (2.16a,b) subject to (3.1). For this, we applied the second-order upwind finite-difference scheme detailed in appendix D. A typical progression is shown in figure 2 for the case $m=1$ for increasing times of $(At) = 1, 2, \ldots , 5$, where the initial shape is shown red. The solution illustrates the progressive rounding of the tip with time superposed with an approximately uniform global recession rate. The evolution of the tip $h_0(t) \equiv h(0,t)$ is shown for $m=0.2, 1$ and $5$ in figure 2$(b)$. In each case, the tip descends in proportion to $t^{4/5}$ for all time with a prefactor that increases with steepness $m$. The steeper case dissolves faster because the larger along-slope component of buoyancy drives a faster boundary-layer flow. The persistent $|h_0| \propto t^{4/5}$ power law indicates the existence of self-similar evolutions arising for every value of the initial slope $m$.

Figure 2. Numerical solutions to the full model (2.16a,b) subject to the initial condition (3.1). Panel $(a)$ shows the evolution for the case in which the shape is initialised as two intersecting planes of slope $m=1$. The progression is shown at (scaled) times of $At=0, 1, \ldots , 5$, illustrating the rounding of the tip. The initial shape is shown in red. Panel $(b)$ shows tip displacement with time for $m=0.2, 1$ and 5, demonstrating the power law $|h_0| \sim (At)^{4/5}$ implied by the similarity scaling of (3.2).

3.1. Similarity solutions

We begin by attempting to construct an intrinsic length scale from the two equations (2.16a,b) and (3.1) using scaling analysis. Depending on whether this length scale must incorporate time $t$ explicitly, one can assess whether the system is likely to support self-similar solutions. The initial condition (3.1) yields the direct scaling between horizontal and vertical length scales, $h \sim x$. The governing equation (2.16a,b) yields the further scaling relationship

(3.2)\begin{equation} x \sim h \sim (At)^{4/5}, \end{equation}

indicating that the evolution of the body remains invariant on scaling both the vertical and horizontal dimensions by a factor of $(At)^{4/5}$. With all the scaling relationships in the system now applied, we deduce that no intrinsic time-independent length scale can be constructed. Thus, the system is anticipated to support self-similar evolutions. Since the system depends on the dimensionless number $m$ in (3.1), different similarity solutions can be anticipated to arise for each value of $m$. To explore these solutions, let us recast the system in terms of the similarity variables motivated by (3.2), namely,

(3.3a,b)\begin{equation} x = (At)^{{4}/{5}} \xi, \quad h(x,t) = (At)^{4/{5}} f(\xi), \end{equation}

where $\xi$ is the similarity coordinate and $f(\xi )$ is the shape function. The governing equation (2.16a,b) and initial condition (3.1) transform to

(3.4a,b)\begin{equation} \frac{4}{5} \left(\,f - \xi f'\right) = \frac{- |\, f' s'^2 |^{1/3} }{\displaystyle \left( \int_0^{\xi} |\, f' s'^2 |^{1/3} \,{\textrm{d}} \xi \right)^{1/4}}, \quad f(\infty) \sim{-}m \xi, \end{equation}

respectively, where $s' = \sqrt{1 + (\,f')^2}$. The large-$\xi$ asymptotic condition (3.4) arises from the initial condition (3.1) because the limit $t \to 0$ corresponds to the limit $\xi = (At)^{-4/5} x \to \infty$ for all $x>0$. Physically, this condition applies because the dissolution rate tends to zero as the boundary-layer thickness grows with $x$. Hence, the shape asymptotes to its initial shape in the limit $x \to \infty$. The similarity system above forms an ordinary integro-differential system for $f(\xi )$, giving a different solution for each value of $m$.

We solved (3.4a,b) numerically using a shooting method. For a given value of $m$, we begin by choosing a trial value of the unknown tip height, $f_0 \equiv f(0)$. We then integrate (3.4a) forwards from this trial value to a large value of $\xi$, denoted $\xi _{\infty }$. The integration was conducted using a Runge–Kutta integrator. Since (3.4a) cannot be expressed explicitly in terms of $f'$, we employed the implicit Matlab integrator ode15i. The numerically predicted value of the gradient $f'(\xi )$ was then compared with the far-field value $-m$ required by (3.4b). Depending on whether the value was greater than or less $m$, a new trial value of $f_0$ was formulated. The bisection method was applied until the error $|\,f'(\xi _{\infty }) + m|<10^{-6}$, for typical values of the integration length taken as $\xi _{\infty } > 10^3\ m^{-1}$.

Our resulting numerical solutions are shown in figure 3. Figures 3(a) and 3(c) illustrate the similarity solution arising for the shallow slope $m=0.2$ and the steep slope $m=5$, respectively. The initial shapes are shown as red curves. The self-similar shapes exhibit rounded tips, with the steeper case retaining a relatively sharper tip that descends faster. In accordance with the similarity scaling of (3.3a,b), the tip descends according to

(3.5)\begin{equation} h_0(t) = f_0 (At)^{4/5}, \end{equation}

where $f_0 \equiv f(0)$ is the tip-descent prefactor. The plot of the tip-descent prefactor $|\,f_0|$ is shown as a function of $m$ in figure 3. The prefactor $|\,f_0|$ increases with $m$, consistent with steeper shapes receding faster. The power-exponent of $f_0$ switches from $2/5$ for small $m$ to $4/5$ for large $m$, indicating the existence of distinct asymptotic regimes arising for shallow and steep dissolving bodies.

Figure 3. Similarity solutions describing the shape evolution of a body initially comprised of two intersecting planes, determined by numerically solving the ordinary integro-differential system (3.4). Panels $(a,c)$ show the solutions for the shallow case $m=0.2$ and the steep case $m=5$, respectively. Panel $(b)$ shows the universal relationship between m and the tip-descent prefactor $|\,f_0|$ which controls the rate of descent of the tip in accordance with (3.5). The asymptotes in the limits of $m \to 0$ and $m \to \infty$ are shown as a black dotted line and a black dashed line, respectively, as predicted by the reduced asymptotic models derived in § 4. The locations of the examples in $(a,c)$ are indicated by red circles.

Before considering these limiting regimes, we consider the tip sharpness predicted by the self-similar evolutions. In terms of (3.3a,b), the tip sharpness is

(3.6)\begin{equation} S(t) = |h_{xx}| = \Sigma (At)^{{-}4/5}, \end{equation}

where $\Sigma \equiv |\,f''(0)|$, implying that the surfaces blunt progressively as $t^{-4/5}$. This is consistent with the blunting exhibited by our numerical solution in figure 2. Surprisingly, dissolving shapes do not always blunt with time, a result that will be demonstrated in the companion paper (Pegler & Davies Wykes Reference Pegler and Davies Wykes2021). It will be shown in particular that the tips of more general shapes (e.g. initially rectangular or cylindrical objects) instead sharpen progressively.

Recasting (2.19) in similarity form, we obtain

(3.7)\begin{equation} \Sigma = \frac{3}{4} \left( \frac{4}{5} f_0 \right)^4, \end{equation}

which allows us to evaluate $\Sigma$ automatically in terms of the already determined $f_0$. As shown in figure 4, the sharpness prefactor increases with $m$, implying that initially steeper shapes retain sharper tips. It should be noted that slope and sharpness need not always increase together, as will be demonstrated in the forthcoming paper.

Figure 4. The tip-sharpness prefactor $\Sigma$ defined by (3.6) as a function of the initial slope $m$, determined using (3.7). The shallow and steep asymptotic predictions of (3.13) and (3.20b) are shown by the dotted and dashed black lines, respectively.

3.2. The shallow limit

If the slope is small, $h_x \ll 1$, then the along-slope component of gravity is given to leading order directly by the slope of the surface, $\sin \alpha \approx -h_x$. Under this approximation, the governing equation (2.16a,b) simplifies to

(3.8)\begin{equation} \frac{\partial h}{\partial t} \approx \frac{-A|h_x|^{1/3}}{\displaystyle\left( \int_0^x |h_x|^{1/3} \,{\textrm{d}} x \right)^{1/4}}, \end{equation}

producing a simplified shallow version of the general theory. In this limit, the similarity system (3.4a,b) reduces to

(3.9a,b)\begin{equation} \frac{4}{5} \left(\,f - \xi f'\right) = \frac{-|\,f' |^{1/3}}{\displaystyle \left( \int_0^{\xi} |\,f' |^{1/3} \,{\textrm{d}} \xi \right)^{1/4}}, \quad f(\infty) \sim{-}m \xi. \end{equation}

With the approximation of shallowness applied, it is now possible to eliminate the initial slope $m$ from the system. To do this, we recast the reduced system above in terms of the scaled variables

(3.10a,b)\begin{equation} \eta = m^{3/5} \xi, \quad \phi = m^{{-}2/5} f, \end{equation}

(3.9a,b) become

(3.11a,b)\begin{equation} \frac{4}{5} \left( \phi - \eta \phi'\right) = \frac{-|\phi ' |^{1/3}}{\displaystyle \left( \int_0^{\eta} |\phi' |^{1/3} \,{\textrm{d}} \eta \right)^{1/4}}, \quad \phi(\infty) \sim{-}\eta, \end{equation}

which has rendered the system free of the only parameter $m$. The solution to this reduced system thus yields a universal similarity solution describing the shape of dissolving shallow bodies. The solution, determined using the same numerical scheme as that used for the full model, is shown as a dashed black curve in figure 5 alongside a selection of solutions to the full model (3.4) with decreasing $m$, illustrating the approach of the full solutions towards the shallow regime as $m \to 0$. The (scaled) tip prefactor determined from the numerical solution is $\phi _0 \approx -1.49$. Using (3.10a,b), we determine the tip-descent prefactor as

(3.12)\begin{equation} |\,f_0| \sim 1.49\, m^{2/5}. \end{equation}

This result is shown as a dotted black line in figure 3, where it is confirmed to describe the tip-descent prefactor as $m \to 0$.

Figure 5. The universal limiting similarity solution arising for shallow slope ($m \to 0$), obtained by solving (3.11a,b) numerically (dashed black curve) plotted in terms of the rescaled coordinates (3.10a,b). Numerical solutions to the full similarity system (3.4a,b) (blue curves) are shown for $m=1.2, 0.8$ and $0.2$, illustrating convergence to the shallow theory as $m \to 0$.

Substituting (3.12) into (3.6), we obtain the prediction for the sharpness prefactor

(3.13)\begin{equation} \Sigma = 1.50 \, m^{8/5}. \end{equation}

The sharpness of tips formed by dissolving shallow bodies therefore scales approximately as the square of the slope scale, $m^{1.6}$.

3.3. The steep limit

Let us now consider the converse limit of steep surfaces, $h_x \gg 1$, given by $m \to \infty$. In this limit, the along-slope component of gravity acts vertically to leading order with negligible dependence on the slope, with $\sin \alpha \approx 1$ and $s_x \approx |h_x|$. In this limit, the governing equation (2.16a,b) simplifies to

(3.14)\begin{equation} \frac{\partial h}{\partial t} = \frac{-A|h_x|}{\displaystyle\left( \int_0^x |h_x| \,{\textrm{d}} x \right)^{1/4}}, \end{equation}

forming a steep version of the theory. In this case, the dissolution pattern is given by $V \approx A [h_0(t) - z]^{-1/4}$, where $h_0(t)$ is the tip position, as derived from (2.15) on setting $\alpha \approx {\rm \pi}/2$. Thus, the dissolution pattern does not depend on the slope $h_x$ at leading order. However, a non-trivial time dependence of the problem is retained in the coupling between dissolution and the tip position $h_0(t)$. To see this explicitly, it is helpful to recast the evolution equation (3.14) in terms of the horizontal position of the boundary, $x=X(z,t)$, which yields

(3.15a)\begin{equation} \frac{\partial X}{\partial t} ={-}V = \frac{-A}{(h_0(t) - z)^{1/4}}, \end{equation}

where

(3.15b)\begin{equation} X(h_0, t) = 0. \end{equation}

Thus, the dissolution pattern evolves in response to descent of the tip $h_0(t)$ which is, in turn, coupled to the surface profile via (3.15b).

The similarity form of (3.14), along with (3.4b), yields the reduced system

(3.16a,b)\begin{equation} \frac{4}{5} \left(\,f - \xi f'\right) = \frac{-|\, f'|}{\displaystyle \left( \int_0^{\xi} |\, f' | \,{\textrm{d}} \xi \right)^{1/4}}, \quad f(\infty) \sim{-}m \xi. \end{equation}

Recasting these equations in terms of the new set of variables

(3.17a,b)\begin{equation} \eta = m^{1/5} \xi, \quad \phi = m^{{-}4/5} f, \end{equation}

we obtain

(3.18a,b)\begin{equation} \frac{4}{5} \left( \phi - \eta \phi'\right) =\frac{ - | \phi'|}{\displaystyle \left( \int_0^{\eta} | \phi' | \,{\textrm{d}} \eta \right)^{1/4}}, \quad \phi(\infty) \sim{-} \eta. \end{equation}

The new system above is independent of $m$, implying that there exists a universal limiting steep solution subject to scaling. This solution can be determined analytically (appendix E) as

(3.19)\begin{equation} \frac{4}{5} (-\phi_0)^{1/4} \eta = \frac{\phi}{\phi_0} \int_{1}^{\phi/\phi_0} \frac{\chi^{{-}2}}{(\chi - 1)^{1/4}} \,{\textrm{d}} \chi, \end{equation}

where $\phi _0 = -(5 {\rm \pi}/8 \sqrt {2} )^{4/5} \approx -1.300$ is the prefactor to the descent position. The solution (3.19) is plotted in terms of the original similarity coordinates $(\xi ,f)$ as a dashed blue curve in figure 6$(a)$ where, at the scale shown, it effectively overlays the full numerically determined solution shown as a solid black curve. The implied asymptote for the descent position and sharpness prefactor (3.6) are therefore

(3.20a,b)\begin{equation} f_0 \sim 1.300 \, m^{4/5}, \quad \Sigma = 0.878 \, m^{16/5}, \end{equation}

respectively. The descent position (3.20a) is shown as a dashed black line in figure 3$(b)$, and matches the numerically determined prefactor for $m \to \infty$. The result of (3.20b) implies a rapid near cubic increase of the sharpness with the slope scale $m^{3.2}$, implying a particularly rapid increase in tip sharpness with the initial slope.

Figure 6. The asymptotic structure of the steep similarity solutions, as detailed in § 3.3, shown for the illustrative steep case $m=5$. The plots show, at different scales, the full numerical solution to (3.4) (black solid curve). In panel $(a)$, the inner solution (3.19) (dashed blue curve) is compared alongside the full numerical solution at a large scale, illustrating the transition over a region of $O(m^{-1/5})$ from an outer region wherein the shape asymptotes towards the initial profile to a region in which the shape assumes a $4/3$-power. In panel $(b)$, the full solution is compared alongside the solution to the inner–inner (3.24) (red dotted curve). Through this inner–inner sublayer, the profile transitions from the $4/3$-power shape to a parabola on a scale of $O(m^{-16/5})$. The variable ${\rm \Delta} f(\xi ) = f(\xi ) - f_0$ represents the shape of the body zeroed at the tip, where $f_0 = f(0)$. In plotting the inner–inner solution, we have chosen to fix the reference tip position $f_0$ using the numerically determined solution, in order to show a clear comparison between the near-tip shape of the inner–inner region and the numerical solution in panel $(b)$.

According to (3.19), the shape near the tip is given by

(3.21)\begin{equation} \phi \sim \phi_0 - \left(-{\frac{3}{5}} \phi_0 \eta \right)^{4/3} \quad (\eta \to 0), \end{equation}

indicating a $4/3$-power shape. The result of (3.21) implies infinite tip sharpness ($\phi ''(0) = \infty )$ and appears to conflict with the general property derived in § 2.2 that the tip is parabolic. The discrepancy arises because the assumption of a steep slope ($h_x \gg 1$) underlying (3.19) and (3.21) breaks down as $\eta \to 0$ as the profile transitions to being locally horizontal. Indeed (3.21) predicts that $\phi ' \to 0$ as $\eta \to 0$, self-contradicting the assumption that $h_x \gg 1$ that underlies it. The steep limit $m \to \infty$ is therefore singular.

3.3.1. Matching to the parabolic tip through an inner non-steep region

To resolve the conflict between (2.18) and (3.21), we consider an asymptotic sublayer near the tip through which we allow the surface to be non-steep, $|h_x| \sim 1$. Let $\eta \sim \varepsilon (m)$ denote the as-yet-undetermined horizontal extent of this non-steep sublayer. Recasting the full similarity equation (3.4) in terms of the steep similarity variables (3.17a,b), but now retaining all the terms, we obtain

(3.22)\begin{equation} \frac{4}{5} \left( \phi - \eta \phi' \right) = {\displaystyle \frac{\displaystyle -|\phi' (m^{{-}2} + \phi'^2)|^{1/3}}{\displaystyle \left( \int_0^{\eta} |\phi' (m^{{-}2} + \phi'^2 )|^{1/3} \,{\textrm{d}} \eta \right)^{1/4}}}. \end{equation}

Here, the $m^{-2}$ terms can be interpreted as the non-steep correction to the steep solution (3.19), with these terms intervening wherever $\phi ' \lesssim m^{-1}$. Neglect of the $m^{-2}$ terms in (3.22) indeed recovers (3.18a,b). To determine the size of the non-steep sublayer, we note that the steep outer (3.21) predicts its own inconsistency once $\eta \sim m^{-3} \equiv \varepsilon$. The relevant inner–inner coordinate $\zeta$ and variable $\psi$ can now be written down as

(3.23a,b)\begin{equation} \eta = m^{{-}3} \zeta, \quad \phi = \phi_0 + m^{{-}4} \psi(\zeta), \end{equation}

where the latter follows on considering where $\phi ' \sim m^{-1}$. Substituting (3.23a,b) into (3.22) and neglecting higher-order terms, we obtain the inner–inner equation

(3.24)\begin{equation} \frac{4}{5} \phi_0= {\displaystyle \frac{\displaystyle -|\psi' (1 + \psi'^2)|^{1/3}}{\displaystyle \left( \int_0^{\zeta} |\psi' (1 + \psi'^2 )|^{1/3} \,{\textrm{d}} \zeta \right)^{1/4}}}. \end{equation}

The single neglected term in this inner equation compared to the full equation (3.22) is the $\eta \phi '$ term. Thus, the solution undergoes an asymptotic transition from an outer steep region that is simplified by neglect of the non-steep correction $m^{-2}$, to an inner non-steep region that is simplified by neglect of the $\eta \phi '$ term.

The solution to the inner (3.24) provides the universal form of a tip matching to a steep body upstream. The solution to (3.24) for $\psi (\zeta )$ was determined numerically using a modification of our numerical scheme used for (3.4). To illustrate the general asymptotic structure, we show in figure 6, at two levels of magnification, the full numerical solution to (3.4) for $m=5$ (black solid curve), plotted in terms of the original similarity variables $(\xi , f)$. In figure 6$(a)$, the full solution is compared alongside the inner solution (3.16a,b) (blue dashed curve). In figure 6$(b)$, the full solution is compared with the inner solution to (3.24) (red dotted curve). To provide an effective comparison of the near-tip shape in the inner–inner region in figure 6$(b)$, we set the height of the tip of the inner–inner solution to coincide with that of our numerically determined value of $f_0$. The plot shows the switch from the intermediate $4/3$-shape in the inner limit of the inner solution (3.19) to a parabola near the tip. In terms of the original similarity variables, the size of the non-steep zone at the tip, derived using (3.23a,b) and (3.17a,b), has equivalent horizontal and vertical dimensions given by $\xi \sim m^{-16/5}$. The $4/3$-power intermediate shape (3.21), applicable over the region $m^{-16/5} \ll \xi \ll m^{-1/5}$, is a factor of $m^{3}$ larger than the parabola in $\xi \ll m^{-16/5}$.

4. Axisymmetric bodies

The dissolution of axisymmetric bodies is governed by different equations and produces different similarity solutions. We define the surface of an axisymmetric body by $z = h(r,t)$, where $r$ is the horizontal distance from the vertical axis of the body. The boundary-layer equations describing axisymmetric natural convection read (Merk & Prins Reference Merk and Prins1953)

(4.1)\begin{align} \frac{\partial (ru)}{\partial s} + \frac{\partial (rv)}{\partial y} &= 0, \end{align}

(4.2)\begin{align} \rho \left( u \frac{\partial u}{\partial s} + v \frac{\partial u}{\partial y} \right) &= \frac{\partial }{\partial y}\left( \mu \frac{\partial u}{\partial y} \right) + {\rm \Delta} \rho g \sin[\alpha(s,t)] C, \end{align}

(4.3)\begin{align} u \frac{\partial C}{\partial s} + v \frac{\partial C}{\partial y} & = \frac{\partial }{\partial y} \left( \kappa \frac{\partial C}{\partial y} \right), \end{align}

where the coordinate system $(s,y)$ is defined similarly to that used for the two-dimensional case, with $s$ representing the arc length moving radially along the outside of the body, and $y$ the normal distance from the surface. The only difference between these equations and those of (2.9)–(2.11) is the presence of $r$ in (4.1). One interesting implication is to introduce a new dependence of the boundary-layer system on the surface profile within the governing boundary-layer system. Even for steep bodies, an axisymmetric dissolution pattern depends on the shape profile and not only the tip position, as applied in the two-dimensional case (3.15). The presence of $r$ in (4.1) also precludes the transformation of Acrivos (Reference Acrivos1960a) from eliminating the surface profile, as was required to derive the dissolution pattern (2.15) (appendix C). However, it is possible to transform the axisymmetric equations first onto a new system (Mangler Reference Mangler1948) that is equivalent in mathematical form to the two-dimensional equations. On applying these transformations, and then recasting the surface mass diffusion rate back in terms of the original variables (appendix F), the axisymmetric analogue of (2.14) is determined as

(4.4)\begin{equation} q = \frac{ D (r\sin \alpha )^{1/3} }{\displaystyle \left(\int_0^s (r^4 \sin \alpha)^{1/3} \,{\textrm{d}} s \right)^{1/4} }, \end{equation}

derived in Acrivos (Reference Acrivos1960b). Compared to its two-dimensional analogue (2.14), this equation includes new dependences on the radial coordinate $r$ representing the effects of radial spreading.

Substituting (4.4) into (2.7), we obtain the new asymmetric model for shape evolution

(4.5a,b)\begin{equation} h_t = \frac{ -A |rh_r s_r^2|^{1/3} }{\displaystyle \left( \int_0^r |r^4 h_r s_r^2|^{1/3} \,{\textrm{d}} r \right)^{1/4} }, \quad s_r = \sqrt{1 + h_r^2}, \end{equation}

forming the axisymmetric version of (2.16a,b). The structure of (4.5a,b) is similar to the two-dimensional equation (2.16a,b), with the exception that it includes new factors of $r$ representing the effects of radial spreading. These new factors cancel in a scaling sense, implying that all the fundamental similarity scalings discussed in § 3.1 are preserved. Indeed, this could have been anticipated since the original introduction of $r$ in (4.1) did not introduce any new scalings into the problem. Although the presence of $r$ does not alter the fundamental temporal scales of evolution, it does alter the cross-section shapes of the evolving body.

Repeating the analysis of the admissible tip shape and tip-descent rate conducted in § 2.2 for the axisymmetric equation (4.5a,b), we obtain the near-tip ($r \to 0$) asymptotic shape

(4.6)\begin{equation} h \sim h_0(t) - {\frac{1}{2}} S(t) r^2. \end{equation}

Thus, the shape of an axisymmetric tip dissolving under buoyancy is paraboloidal. Moreover, we obtain the axisymmetric tip-descent law

(4.7)\begin{equation} \frac{\partial h_0}{\partial t} ={-}A \left( \frac{8}{3} S(t) \right)^{1/4}. \end{equation}

The prefactor to the tip-descent law (4.7) is a factor of $2^{1/4} \approx 1.2$ larger than that of the corresponding two-dimensional result of (2.17), implying that axisymmetric tips dissolve approximately 20 % faster for the same cross-sectional tip curvature $S(t)$.

To determine the shape functions of the axisymmetric similarity solutions, we recast (4.5a,b) in terms of the similarity variables, $r = (At)^{{4}/{5}} \xi$ and $h(r,t) = (At)^{4/{5}} f(\xi )$, giving

(4.8a,b)\begin{equation} \frac{4}{5} \left(\,f - \xi f'\right) = \frac{- | \xi f' s'^2 |^{1/3} }{\displaystyle \left( \int_0^{\xi} | \xi ^4 f' s'^2 |^{1/3} \,{\textrm{d}} \xi \right)^{1/4}}, \quad f(\infty) \sim{-}mr, \end{equation}

forming the axisymmetric analogue of (3.4). We solve this equation using a modification of the numerical scheme used for (3.4) to incorporate the new $\xi$ terms. The resulting shape functions for cones of initial slope $m=0.2$ and 5 are shown in figure 7(a,c). The shapes are qualitatively similar to those found for the two-dimensional solutions shown in figure 3. The most appreciable difference is the alteration of the descent prefactor $|\,f_0|$ shown in figure 7$(b)$. The descent prefactor is, across the complete range of $m$, a near-universal factor of ${\approx }1.16$ larger than those of the corresponding two-dimensional values of $|\,f_0|$ shown in figure 3$(b)$. Axisymmetric bodies thus dissolve approximately 16 % faster for the same slope scale $m$. The similarity form of the tip-descent law (4.7) is given by

(4.9)\begin{equation} \Sigma = \frac{3}{8} \left( \frac{4}{5} f_0 \right)^4, \end{equation}

which is one half the value of the two-dimensional analogue (3.7). Axisymmetric tips are therefore half as sharp as two-dimensional tips with the same descent prefactor $f_0$.

Figure 7. Similarity solutions describing the shape evolution of an initially axisymmetric conic body, determined by numerically solving the ordinary integro-differential system (4.8). Panels $(a,c)$ show the solutions for the shallow case $m=0.2$ and the steep case $m=5$, respectively. Panel $(b)$ shows the general relationship between tip-descent prefactor $|\,f_0|$ and the slope $m$ for an initially conic body. The asymptotes in the limits of $m \to 0$ and $m \to \infty$ are shown as a black dotted line and a black dashed line, respectively, as predicted by the reduced asymptotic models derived in § 4. The locations of the examples in $(a,c)$ are indicated by red circles. The solutions are determined by numerically solving the ordinary integro-differential system (4.8), forming the axisymmetric analogue of figure 3.

We develop shallow and steep asymptotic limits as follows. The limit of a shallow axisymmetric body ($h_r \ll 1$) yields

(4.10)\begin{equation} \frac{4}{5} \left(\,f - \xi f'\right) = \frac{- | \xi f' |^{1/3} }{\displaystyle \left( \int_0^{\xi} | \xi ^4 f' |^{1/3} \,{\textrm{d}} \xi \right)^{1/4}}. \end{equation}

Following the same rescaling of (3.11), we obtain $\phi _0 \approx -1.726$. The associated asymptote $f = -1.726\, m^{2/5}$ is shown as a dotted line in figure 7$(b)$, and accurately represents the small-$m$ limit for $m < 0.5$.

The corresponding similarity equation in the steep limit ($h_r \gg 1$) is given by

(4.11)\begin{equation} \frac{4}{5} \left(\,f - \xi f'\right) = \frac{- \xi|\, f' | }{\displaystyle \left( \int_0^{\xi} \xi ^4|\, f' | \,{\textrm{d}} \xi \right)^{1/4}}. \end{equation}

Unlike the two-dimensional case, the steep axisymmetric model of (4.11) cannot be integrated analytically. This is a reflection of the fact that the axisymmetric dissolution pattern depends on the surface profile as a consequence of radial spreading from the tip. Recasting (4.11) in terms of the scaled variables (3.17a,b) and solving the resulting system numerically, we obtain the universal steep solution applicable to an axisymmetric body. The resulting prefactor is $\phi _0 \approx 1.505$ and hence $f_0 \sim 1.505 \, m^{4/5}$, which is shown as a dashed black line in figure 7. Trying a power-law solution of the form $\phi = \phi _0 + c \eta ^a$ in (4.11) and taking the limit $\eta \to 0$, we obtain

(4.12)\begin{equation} \phi \sim \phi_0 - \left(-{\frac{3}{5 \sqrt{2} }} \phi_0 \eta \right)^{4/3} \quad (\eta \to 0), \end{equation}

forming the modified axisymmetric analogue of (3.21). Thus, the shape of a steep axisymmetric body also forms a $4/3$-power in the lead in to the tip with a slightly modified prefactor. To resolve the transition to a parabolic tip with a finite sharpness prefactor (3.20), one could construct the same kind of asymptotic framework as was done in the two-dimensional case in § 3.3.1. The same essential scalings for the size of the inner non-steep region arise, and thus we omit the analysis here for brevity.

5. Laboratory study

In order to test the predictions of our theoretical model, we conducted a series of laboratory experiments. For this, we used bodies of solidified sugar (known as rock candy or boiled sweets) immersed in freshwater, as illustrated schematically in figure 8 and as a sequence of photographs in figure 9. The experiment was contained in a square-based tank of length 62.5 cm, width 62.5 cm and height 61 cm. The tank was backlit using LED light tape. The top of the tank was covered to suppress evaporative cooling. The water temperature was measured during experiments using an alcohol thermometer inside the tank. The temperature remained between 23 and $25\, ^\circ \textrm {C}$ during each experiment.

Figure 8. Schematic of our experimental apparatus, as described in § 5. An alcohol thermometer was used to measure the temperature of the water inside the tank. The top of the tank was covered to suppress evaporative cooling during experiments. Two stands were used to raise the candy cone off the floor of the tank: the one illustrated in this figure was 12 cm high and was used for the relatively shorter cones used for experiments A and B, whereas a stand of 2.5 cm was used for the relatively taller cones used for experiments C and D. The stand was used to create a space for the heavy solution to collect at the floor of the tank so that the candy remains immersed in fresh ambient fluid.

Figure 9. Dissolution of a candy cone at 40 min intervals, illustrating the descent of the tip. The images are from experiment C, which had an initial slope of $m = 12.3$. A scale bar is provided on the left, indicating descent rates of several cm per hour. Supplementary movie of this experiment is available at https://doi.org/10.1017/jfm.2020.507.

The candy was made from an 8 : 3 : 2 volume mixture of sugar, corn syrup, and water, and coloured using blue food dye. The mixture was heated to $150\, ^\circ \textrm {C}$, representing the ‘hard crack’ temperature at which the solution has crossed the glass transition and will set as solid candy. The molten candy was then sand cast, using a plastic conic mould and left to cool and harden.

In order to prevent the candy from moving or rocking over the course of the experiment, the base of each cone was glued to a plastic dish (10 cm diameter) using silicone sealant. The dish was rested on a raised stand 12 cm from the base of the tank for the relatively shorter cones used for experiments A and B, and a piece of PVC to a height of 2.5 cm for the relatively taller cones used for experiments C and D. This was done in order to provide a region below the base of the cone where dissolved sugar can collect, thus suppressing the development of a stratification above the base of the cone that would change the ambient conditions. The experiment was recorded using a camera facing the front of the experimental tank with photographs taken at regular intervals of 30 s.

The convective strength $A$ was estimated for each experiment by measuring the spatial pattern of the rate of change of the surface position in the horizontal direction, $V_r \equiv |\boldsymbol {V} \boldsymbol {\cdot } \pmb{e}_x|$. As illustrated for experiment C in figure 10, the profiles of $V_r(z,t)$ collapse together onto a curve that remains approximately steady in time. For a near-conic surface, (4.4) predicts a dissolution pattern of the form

(5.1)\begin{equation} V_r = A_c \left[h_0(t) - z\right]^{{-}1/4}, \end{equation}

where $A_c$ satisfies

(5.2)\begin{equation} A = \left({\frac{3}{7}}\right)^{1/4} \left( \sin \alpha \right)^{{-}3/2} A_c, \end{equation}

and $\sin \alpha = m / \sqrt {1+m^2}$. The experimental value of $A$ is inferred by fitting (5.1) to the measured dissolution profiles. The resulting curve for experiment C is illustrated by the black dashed curve in figure 10. The experimental dissolution patterns sampled at regular time intervals show good agreement with the predicted pattern of (5.1). For each sampled time, we obtained an empirical estimate of the prefactor $A_c$ and hence the convective strength $A$ via (5.2). A single value of $A$ for each experiment was evaluated by averaging the sampled values. A slight variation in the value of $A$ between different experiments shown in table 1 was likely caused by slight differences in temperature, or an effect of slight non-axisymmetry caused by some leaning of the cone.

Figure 10. The horizontal velocity of the surface of the candy, $V_r$, as a function of distance below the tip ($h_0(t) - z$) for a selection of sampled times over the course of experiment C (coloured solid curves). The horizontal velocity was evaluated using the formula $V_r(z) = {\rm \Delta} r/{\rm \Delta} t$, where ${\rm \Delta} r (z)$ is the horizontal difference in surface profiles between the profiles at times $t$ and $t + {\rm \Delta} t$, where ${\rm \Delta} t = 5\ \textrm {min}$. The value of $|V_r|$ is found for both the left and right sides of the shape separately and then averaged. The function (5.1) is fitted to each time with $A_c$ and $h_0$ treated as fitting coefficients. The dashed line represents the curve produced by the mean value of the values obtained from the individual sampled times, yielding the values listed in table 1.

Table 1. The initial slope $m$ and convective strength $A$ for each of our experiments. The bounds on $A$ represent twice the standard deviation of values of the set of values of $A$ inferred over the course of each experiment, as described in § 5 and illustrated in figure 10.

5.1. Experimental results

We performed four experiments spanning initial cone slopes $m$ from $5.8$ to $17.6$ and convective strengths $A$ from $1.4 \times 10^{-4}$ to $1.6 \times 10^{-4}\ \textrm {cm}^{5/4}\,\textrm {s}^{-1}$, as listed in table 1. Due to small errors in alignment, the cones were angled slightly with respect to gravity (by ${<}1^{\circ }$). In order to correct for small variations in each experiment in the angle of the cone with respect to the camera, the images were rotated by a small angle (around a $1^{\circ }$) such that the cone is oriented vertically in the frame of each image. A sequence of images from experiment C are shown in figure 9, illustrating descent rates of a few cm per hour. Each of the four experiments are available to view as supplementary movies.

The profiles of the solid boundaries were extracted digitally from the experimental images using a thresholding method of the red channel of the RGB images. The evolution of the extracted profiles is shown in the left-hand column of figure 11 with each row representing a different experiment listed in table 1. The profiles are observed to produce relatively smooth shapes that descend with time. In experiment D, a localised notch in the side of the candy was found to grow and eventually cause the tip to break off, as shown in figure 11$(d)$. We believe this phenomenon to be associated with a localised instability whereby a small notch in the side of the dissolving shape rapidly grows in response to accelerated dissolution along the underside of the notch. This may be related to the instability mechanism giving rise to pits on the underside of melting or dissolving bodies (e.g. Vanier & Tien Reference Vanier and Tien1970; Kerr Reference Kerr1994; Mcleod et al. Reference Mcleod, Riley and Sparks1996; Sullivan, Liu & Ecke Reference Sullivan, Liu and Ecke1996; Davies Wykes et al. Reference Davies Wykes, Huang, Hajjar and Ristroph2018). Despite the emergence of this feature, we have nonetheless retained this experiment to illustrate the phenomenon, as well as to show that the descent of the tip is, before the notch extends across the width of the candy, largely unaffected by its emergence.

Figure 11. Experimental results of dissolving candy cones. Successive rows represent experiments A–D, as listed in table 1, corresponding to slopes of $(a\text {--}c)$$m = 5.8$, $(d\text {--}f)$$m = 10.8$, $(g\text {--}i)$$m = 12.3$ and $(\,j\text {--}l)$$m = 17.6$. $(a,d,g,j)$ Show the evolution of the shape extracted digitally from our recordings. $(b,e,h,k)$ Show the same profiles in a coordinate system in which the horizontal and vertical dimensions are each rescaled by $t^{4/5}$. $(c,f,i,l)$ Show the experimental tip position against time (circles) and the associated theoretical prediction (5.3) (red solid line), where a time shift has been applied to account for the fact that the initial shapes are not perfect cones. Insets in $(c,f,i,l)$ show the outlines of the initial shapes for each experiment, illustrating the varying steepnesses.

Naturally, the initial experimental shapes were not perfect cones. In particular, experiments B and D began with relatively blunted tips. To draw a fairer comparison between experiments in view of these blunted tips, we chose the reference height $z=0$ and time origin $t=0$ such as to be consistent with the assumed initial condition of the similarity solution. For this, the location of $z = 0$ was set as the height of the intersection of two lines fitted to either side of the initial shape, corresponding to the initial position of the tip that would apply if the initial shapes were perfectly conic. To allow for small variations between the time that the cones were placed in the tank and the time at which the first images were taken, the time origin was selected by picking a single measurement of the tip height $h_0$ for each experiment and requiring that $t = (h_0/f_0)^{5/4}/A$. The time chosen for this was $30$ minutes after each experiment was initialised. Increasing this time made no appreciable difference to the inferred time origin. The inferred time origins were 2 to 6 minutes before the first image was taken.

As an initial test, we examine whether the experimental shapes collapse onto a self-similar shape on suitable stretching of the coordinates in terms of our similarity variables (3.3). According to the theory, the shape evolves as

(5.3)\begin{equation} h(r,t) = (At)^{4/5} f\left(\frac{r}{(At)^{4/5}} \right), \end{equation}

where $f(\xi )$ is a shape function determined theoretically in § 4. Therefore, if the vertical and horizontal coordinates are both rescaled by $t^{-4/5}$, the outlines of the body should collapse to lie on top of one another. To test this collapse, we plotted the outlines of the experimental shapes with the horizontal and vertical dimensions both scaled by $t^{4/5}$. As shown in the central column of figure 11, each outline collapses to excellent approximation onto a shape that is fixed with respect to the rescaled coordinates. This confirms the essential self-similarity predicted by the theory.

Next, we present direct comparisons between the theoretical prediction for the tip evolution $h_0(t)$ and the experimental results. The theory of § 4 predicts that the tip descends according to

(5.4)\begin{equation} h_0(t) = f_0 (A t)^{4/5} ={-}1.505 \ m^{4/5} (A t)^{4/5}, \end{equation}

where, in view of the steepness of our experiments ($m > 5.8$), we have in the second equation used the prefactor that applies to the limiting steep theory (3.20a).

The tip positions of each experiment inferred using the thresholding method are plotted as functions of time in the right-hand column of figure 11. In each case, we have overlaid the theoretical prediction of (5.3). Generally excellent agreement is observed. To test the universal collapse of the solutions, we have plotted each of the tip evolutions together in figure 12 in forms scaled by $|\,f_0|A^{4/5}$. According to the theoretical prediction of (5.3), the curves should, under this rescaling, all collapse onto the universal curve $t^{4/5}$. An excellent general collapse of the experiments onto this curve is found. With the exception of the choice of time origin used to account for the slightly non-conic initial shapes, no adjustable parameters have been used in this comparison. The agreement indicates that the general theoretical framework we have derived accurately captures the descent rates of dissolving bodies. The self-similar shapes predicted by the axisymmetric theory from § 4 are plotted alongside the collapsed profiles from each experiment A–D in the respective panels of figure 13. Again, there is an excellent agreement between experiment and theoretically predicted shapes.

Figure 12. Collapse of the tip position from all four of our experiments onto a universal curve. The experimental tip positions $|h_0 (t)|$ (shown as symbols) are normalised by the quantity $|\,f_0|A^{4/5}$, where $|\,f_0| = 1.505 \ m^{4/5}$, which, according to the theoretical prediction (5.3), should collapse onto the universal curve $t^{4/5}$ (shown as a black solid curve). A time shift has been applied to each experiment to account for the fact that the initial shapes were not perfect cones. An error estimate was evaluated as twice the error associated with the measurement of $A$ (shown in table 1). The error range was smaller than the sizes of the symbols, and has thus been omitted.

Figure 13. Experimental profiles rescaled into similarity variables, $\xi = (At)^{-4/5}r$ and $f(\xi ) = (At)^{-4/5} h(r,t)$, for each of our experiments A to D (panels a to d, respectively). The theoretical prediction (5.3) is shown as a dashed black curve in each case, with $f(\xi )$ is determined by solving (4.8) using the value of $m$ relevant to each experiment. There is a generally good collapse of the experiments onto the theoretically predicted self-similar shape. With the exception of the choice of time origin, no adjustable parameters were used to formulate this comparison. The initial theoretical shape is shown by the solid blue lines in each case.

6. Summary and discussion

We have developed a model describing the dissolution and melting of solid bodies under free convection that addresses the open question of the shapes formed by dissolving or melting bodies in the case of gravitationally stable free convection, and validated the predictions in the laboratory. Analysis of the model produces explicit predictions for the rates at which melting or dissolving shapes recede. The model development and demonstration of analytical approaches provides a basis for various generalisations of the model to understand dissolving shapes in more complicated and naturally occurring settings, a number of which will be discussed below. Our observed experimental shape profiles collapse to excellent approximation onto a consistent theoretically derived descent position and shape, giving confidence in the essential theoretical framework.

The models for both two-dimensional and axisymmetric bodies were derived by coupling the boundary-layer equations describing stable free convection to an evolution equation for the boundary. In the case of large-$Sc$ number, the model was reduced further to a single integro-differential equation. Analysis of this equation shows that bodies comprised of intersecting planes or cones result in similarity solutions that ‘stretch’ in both the horizontal and vertical dimensions as $t^{4/5}$ with respect to the initial tip position. A general relationship between tip sharpness and descent rate was derived, showing that the tip is always parabolic sufficiently close to the tip and the descent rate is always proportional to the tip curvature to the one-quarter power.

Differing reduced asymptotic descriptions of shallow (near-horizontal) and steep (near-vertical) surfaces were derived. For shallow (near-horizontal) dissolving surfaces with initial slope $m$, the limiting parabolic tip occupies the region $x \ll \ m^{-3/5} (At)^{4/5}$. In the limit $m \to 0$, a full asymptotic description of the body surface $h(x,t)$ as a whole was derived (§ 3.2), which is summarised by

(6.1)\begin{align} (At)^{{-}1/5} h(x,t) \sim \begin{cases} f_0 - 0.75 \, m^{8/5} \xi^2 & \xi \ll m^{{-}3/5}, \\ m^{2/5} \phi(m^{3/5} \xi) & \xi = O(m^{{-}3/5}), \\ -m \xi & \xi \gg m^{{-}3/5}, \end{cases} \end{align}

where $\xi = (At)^{-4/5}x$ is a scaled $x$-position, $f_0 = -1.49 \ m^{2/5}$, and $\phi$ is a universal (parameter-free) function. Shallow bodies thus form a parabolic region that matches to the initial shape far downstream over a length scale of order $m^{-3/5} (At)^{4/5}$.

For a steep body ($m \to \infty$), the shape instead forms a more intricate double-decked boundary layer structure containing both a parabolic and a $4/3$-power near-tip region (§ 3.3). The full structure is summarised by

(6.2)\begin{align} (At)^{{-}1/5} h(x,t) \sim \begin{cases} f_0 - 0.44 \, m^{16/5} \xi^2 & \xi \ll m^{{-}16/5}, \\ f_0 + m^{{-}16/5} \psi(m^{16/5} \xi) & \xi = O(m^{{-}16/5}),\\ f_0 - 0.72\, m^{16/75} \xi^{4/3} & m^{{-}16/5} \ll \xi \ll m^{{-}1/5}, \\ m^{4/5} \phi(m^{1/5} \xi) & \xi = O(m^{{-}1/5}), \\ -m \xi & \xi \gg m^{{-}1/5}, \end{cases} \end{align}

where $f_0 = -1.30 \, m^{4/5}$, and $\psi$ and $\phi$ are universal functions describing the inner–inner and inner decks, respectively (the function $\phi$ here differs from that of (6.1)). The ‘very’ near-tip limit of the inner–inner region, defined by $\xi \ll m^{-16/5}$, forms a parabola. However, within the intermediate region between the inner and the inner–inner, the shape forms a $4/3$ power. The $4/3$ region is considerably larger than the parabolic region (it is a factor of order $m^3$ larger). Thus, while there are two near-tip limits in the case of steep bodies, the $4/3$-shape applies over a much larger region, though the ultimate shape of the tip at the finest scale is parabolic.

The model and its predictions apply generally to situations where the viscosity, diffusivity and density are permitted to vary as any given functions of the concentration (or temperature in the case of melting). Even for mixtures with order-of-magnitude variations in fluid properties over the flow domain, the shapes produced are equivalent. For example, the shape evolution for a melting block of ice is the same as that which arises for a dissolving body of candy in water, despite the viscosity of the air remaining approximately constant in the case of an ice cube, but the viscosity of candy solutions varying by at least two orders of magnitude. The generality stems from a favourable property of similarity solutions in free convection that allows the effects of concentration-dependent variables to be encapsulated in the single parameter $A$. For this reason, equivalent shapes would even arise for two-component effects of compositional convection such as a block of ice melting in salty water, and for flows involving the reversal of buoyancy.

While, as noted above, the results do generalise to fluids with strongly concentration-dependent viscosity, one can question whether the generalisation applies for non-Newtonian fluid with a shear-dependent viscosity. The situation may apply, for example, to the case of a polymer dissolving under free convection into a solution with a different concentration. It was determined by Acrivos (Reference Acrivos1960b) that a direct generalisation of the result of (2.15) can be derived for power-law fluids. Therefore, a generalised form of the high Schmidt number reduction of our model of the shape evolution (2.16a,b) could be straightforwardly developed and analysed using the same techniques considered here. It can be anticipated that the self-similar shapes and temporal exponents of the descent rates associated with the similarity solutions would be altered.

The full model of (2.9)–(2.13) allows for the effects of inertia. However, our analysis focused on the predictions of a reduced model applicable under high Schmidt (or Prandtl) number, which affords the consideration of a purely viscous inner sublayer. The approximation of neglecting inertia in fact remains good even for Schmidt numbers of order unity (Kuiken Reference Kuiken1968), with only minor corrections to the thermal or mass-transfer rates (the adjustment from the $Sc = \infty$ case is just $\approx 20\,\%$ for $Sc = 1$, for example). However, for situations where inertia may be required as a leading-order term in the solutal layer (e.g. for gases), the expression (2.15) can be anticipated to break down. In considering this case, we note that retaining the inertial terms introduces one new scaling into the system that is identical in form to the one already produced by scaling the advective and diffusive terms in (2.11), namely, $u/x \sim \nu / y^2$. Therefore, all the similarity scalings we derived from the simplified high-$Sc$ number models of (2.16a,b) or (4.5a,b) will also apply generally with inertia, i.e. general $Sc$. While the similarity scalings for tip-descent rates would be preserved (e.g. $h_0 \sim t^{4/5}$ for initially conic bodies), the introduction of inertia will alter the evolving shapes. The shapes could be determined by transforming the full model of (2.9)–(2.13) in terms of the similarity variables. The shape could then be determined using a one-off solution of the resulting two-variable heat equation in $(\xi , y)$. Thus, similarity solutions still exist for general $Sc$, but require the consideration of a two-variable system in $(\xi , y)$, rather than the purely ordinary system of the form (3.4) that applies in the large-$Sc$ limit.

An overarching assumption in our analysis is that the boundary layer remains laminar. If the boundary-layer flow becomes unstable and transitions to turbulence, the rates and profiles of convective transfer at the interface will differ from what we have assumed. We were surprised to find that this transition occurs even in the relatively viscous context of our laboratory experiments, as shown in the shadowgraph of one of our experiments in figure 14. Figure 14$(c)$ illustrates the development of a regular pattern of convective rolls which flow along the surface of the body, characterised by a much thicker boundary-layer region compared to the laminar region shown in figure 14$(b)$. Experiments of turbulent composition free convection of ice melting into ambient salty water (Kerr & McConnochie Reference Kerr and McConnochie2015; McConnochie & Kerr Reference McConnochie and Kerr2018) indicate that the rate of heat or solutal transfer in this regime becomes independent of the distance moved downstream of the body, but still depends on the local slope. It is anticipated that an evolution equation analogous to the kind we have developed (2.16a,b) could be derived to couple a transfer law applicable to turbulent free convection to the shape evolution and analysed to derive shape evolutions in this case. The result would likely reveal different shapes and descent laws to those we have reported here. Recent research has also considered the different turbulent regime of a plume flowing along the underside of an ice shelf and maintained by a constant source of buoyancy at the grounding line by an input of relatively fresh subglacial discharge (e.g. Dallaston, Hewitt & Wells Reference Dallaston, Hewitt and Wells2015; Slater et al. Reference Slater, Nienow, Goldberg, Cowton and Sole2017; Hewitt Reference Hewitt2020). This regime differs in being driven by the introduced buoyancy flux as opposed to the cooling by the wall, but similarly results in the potential for interplay between buoyancy-driven flow and shape evolution.

Figure 14. Shadowgraph image of dissolving candy cone. $(a)$ A shadowgraph of experiment C ($m = 12.3$), $(b)$ a section near the top of the cone, where the boundary layer cannot be seen as it remains laminar and $(c)$ shows a section lower down the cone, where intensity fluctuations suggest that the boundary layer has transitioned to turbulence. The image forms a still from a supplementary movie, which illustrates the unsteadiness of the turbulent region.

We speculate that an adaptation of the model could be developed to describe the evolution of solidifying bodies under free convection. A situation of this kind has been used to model the growth of icicles (Makkonen Reference Makkonen1988; Short, Baygents & Goldstein Reference Short, Baygents and Goldstein2006; Camporeale & Ridolfi Reference Camporeale and Ridolfi2012). In this situation, a thin film introduced at the top of the icicle flows down along its exterior and freezes at a spatially variable rate controlled by an upwelling convective boundary layer generated by latent heating (Makkonen Reference Makkonen1988). By imposing a constraint that the shape of the growing icicle remains preserved, Short et al. (Reference Short, Baygents and Goldstein2006) derive a solution that asymptotes to a $4/3$-power shape far away from the tip. This form of shape is interpreted as a far-field asymptote in this case, contrasting with our prediction of (4.12) that a $4/3$-power shape applies in a limit close to the tip. Analysis of the formation of dendrites in the (non-fluid-mechanical) context of a two-dimensional solid solidifying into an ambient under thermal diffusion has revealed the development of parabolic solutions (Nash & Glicksman Reference Nash and Glicksman1979; Ivantsov Reference Ivantsov1985). Similarity with the parabolic regimes we identify appears to be incidental because the formation of the parabola in the context of our results is due to an inherent fluid-mechanical property of free-convective boundary layers that the parabola is the unique shape for which the boundary layer predicts a locally smooth recession rate. In principle, our model could be adapted to describe a solidifying body under free convection, and would likely require additional physics such as surface tension in order to suppress morphological instabilities.

The model we have considered has focused on the case of a purely free-convective (buoyancy-driven) boundary layer. We have neglected any effect of an externally imposed flow, i.e. forced convection, which has been investigated by others (Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013; Huang et al. Reference Huang, Moore and Ristroph2015; Moore Reference Moore2017). In many naturally occurring problems involving melting bodies, such as icebergs (e.g. FitzMaurice, Cenedese & Straneo Reference FitzMaurice, Cenedese and Straneo2017), it can be anticipated that a mixture of forced convection and natural convection can both contribute. Our analysis has addressed one limit of this spectrum. Future theoretical work could consider the extension to allow for mixed (forced and free) convection, potentially with the two acting in differing directions on the body (as considered experimentally by Hao & Tao Reference Hao and Tao2001, Reference Hao and Tao2002). For the context of the shaping of rock spires, the buoyancy-driven dissolution may take the form of a thin-layer flow fed at a constant flux (rainfall) over the surface. This situation creates a buoyancy-driven gravity current over the surface fed by source conditions, a case which could be addressed within the analytical framework developed here.

Finally, we note that the analysis of the theoretical model developed in this paper has focused on situations where the initial shape is conic or formed by two intersecting planes. These cases result in special $t^{4/5}$ similarity solutions for all time. For more general bodies, a rich variety of regime transitions in the temporal exponents of tip-descent rates, occur. This includes the case of initially rectangular bodies, as is representative of the initial shape of ice bergs or ice cubes. These broader themes, including the determination of the intrinsic scales of melting or dissolution of such bodies, will form the focus of our companion paper (Pegler & Davies Wykes Reference Pegler and Davies Wykes2021).