2. Model

Building on the model presented in Part 1, we study a binary alloy exposed at dimensionless time $t=0$ to a heat sink at dimensionless depth $z=0$, with heat loss via a linearised exchange defined below. We neglect salt diffusion and assume the mushy layer is ideal with equal material properties in each phase. We further assume that the fluid is initially of uniform temperature and salinity, and that the initial velocity is uniformly zero. However, for this part, we will focus our attention on the growing mushy layer specifically and relax the assumption that the fluid always stays at rest.

For brevity, we skip the dimensional form of the model described in Part 1 and non-dimensionalise initially using a fixed, time-independent length scale, $\hat {d}$, but convert to a scale defined by the thermal diffusion length scale later in the analysis. The initial length scale has associated time and velocity scales, ${\hat {d}^2}/{\kappa }$ and ${\kappa }/{\hat {d}}$, respectively. A diagram of the non-dimensional system is shown in figure 1.

Figure 1. Non-dimensional model set-up (cf. dimensional form in figure 1 in Part 1). The temperature profile $\theta$ runs between $0$ at the surface heat sink and $\theta _{\infty }$ at depth, while the liquid salinity $S=1$ and its liquidus temperature is ${\theta }_L=1$. The thermal diffusion length scale is proportional to $\sqrt {t}$ and the constant of proportionality for linearised heat transfer is ${\mathcal {B}_i}$. A full description of the non-dimensional parameters used can be found in the main text.

The non-dimensional temperature and salinity are defined as

(2.1a,b)\begin{equation} \theta = \frac{T - T_c}{T_{L\infty} - T_c},\quad S = \frac{\hat{S}-\hat{S}_s}{\hat{S}_\infty-\hat{S}_s}, \end{equation}

for dimensional temperature $T$ and salinity $\hat {S}$. The temperature scale used $\Delta T = T_{L\infty } - T_c$ is the temperature difference between the salinity-dependent freezing temperature of the liquid $T_{L\infty }$ and the surface heat sink at $T_c$. The salinity scale used is the initial difference between the liquid salinity $\hat {S}_\infty$ and pure solid salinity $\hat {S}_s$.

Under this non-dimensionalisation, conservation of energy and salt yields

(2.2a,b)\begin{equation} \frac{\partial{\theta}}{\partial{t}} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\theta + {\mathcal{S}_t}\frac{\partial{\chi}}{\partial{t}} = \nabla^2\theta,\quad \frac{\partial{\chi S}}{\partial{t}} =- \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}S, \end{equation}

where we have added advective terms to the equations presented in Part 1 with Darcy velocity $\boldsymbol {u}$ describing the volume flux of liquid per unit area. These equations introduce the volumetric liquid fraction $\chi$, and the Stefan number ${\mathcal {S}_t} = {L}/{c_p(T_{L\infty }-T_c)}$ which represents the ratio of latent heat $L$ to specific heat in the system (with specific heat capacity $c_p$). The temperature and salinity are constrained to each other via the liquidus relationship

(2.3)\begin{equation} \theta = 1 - \mathcal{C}(S-1),\quad \mbox{where}\ \mathcal{C}= \frac{\varGamma(\hat{S}_\infty-\hat{S}_s)}{T_{L\infty} - T_c} \equiv \frac{\hat{S}_\infty-\hat{S}_s}{\hat{S}_c-\hat{S}_\infty}, \end{equation}

where $\hat {S}_c$ is the salinity of fluid with liquidus temperature $T_c$ and $\varGamma$ is the dimensional liquidus slope. The dimensionless concentration ratio $\mathcal {C}$ measures the salinity effects due to freezing-point depression, with a larger value indicating a larger impact on the solidification dynamics.

To complete the set of dynamic equations, we use Darcy's law to describe momentum transfer within the porous mush, and note that the fluid is incompressible when solid and liquid densities are assumed equal, yielding

(2.4a–c)\begin{equation} \boldsymbol{u} =-\mathcal{R}_{m}\chi^3[\boldsymbol{\nabla}p + \theta\boldsymbol{{e_z}}],\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0,\quad \mathcal{R}_{m} = \frac{g (\alpha + {\beta}/{\varGamma}) \Delta T\hat{\varPi}_0\hat{d}}{\kappa\nu}. \end{equation}

In the first of these equations, $p$ is the dimensionless pressure, $\boldsymbol {{e_z}}$ is a vertical unit vector pointing downwards and we assume the permeability is $\hat {\varPi }_0\chi ^3$ for constant $\hat {\varPi }_0$, which approximates a measured correlation of sea-ice permeability ${\propto }\chi ^{3.1}$ (Freitag Reference Freitag1999; Notz & Worster Reference Notz and Worster2009). The density anomalies driving convection $\Delta \rho = \rho _0[ \beta (\hat {S}-\hat {S}_c) -\alpha (T -T_c)]$ are assumed to satisfy a linear equation of state with constant values of the thermal expansion coefficient $\alpha$, solutal contraction coefficient $\beta$ and reference density $\rho _0$. The liquidus relationship has been used to express saline buoyancy contributions in terms of the dimensionless temperature. The mushy Rayleigh number $\mathcal {R}_{m}$ describes the relative strength of buoyancy forcing from both thermal and saline expansion/contraction vs the rate of viscous and thermal dissipation.

Following Emms & Fowler (Reference Emms and Fowler1994) we treat the mush–liquid interface as an isothermal, constant-pressure boundary at the freezing temperature which yields

(2.5a,b)\begin{equation} \theta = 1,\quad \boldsymbol{{e_x}}\boldsymbol{\cdot}\boldsymbol{u} = 0,\quad \text{when}\ {z = h}, \end{equation}

where $\boldsymbol {{e_x}}$ is a horizontal unit vector and we have used (2.4a) to link $p$ and $\boldsymbol {u}$. Here, the isothermal condition arises from a marginal equilibrium condition, with the temperature satisfying the liquidus constraint with an assumed uniform salinity in the liquid. Note that the liquid fraction $\chi \rightarrow 1$ at the mush–liquid interface $z\rightarrow h$ in the background state. Emms & Fowler (Reference Emms and Fowler1994) justify use of a constant-pressure boundary condition when the resistance to flow in the liquid region is much smaller than in the porous region, which here occurs when the Darcy number $Da = \hat {\varPi }_0/\hat {d}^2$ is small.

Our treatment of the flow boundary conditions assumes a sharp transition between a pure liquid in $z>h$ and a porous medium in $z< h$. Le Bars & Worster (Reference Le Bars and Worster2006) analysed the transition between purely liquid and porous regions using the Darcy–Brinkmann equation. They found that a split-domain approach with Navier–Stokes and Darcy flow in separate regions is an accurate approximation if the flow boundary conditions at a mush–liquid interface are characterised by instead imposing continuity of pressure and velocity at a displaced dimensional distance $\hat {\delta }=\sqrt {\hat {\varPi }/\chi }$ inside the porous media, where $\hat {\varPi }$ is the dimensional permeability. Applying the permeability used in (2.4a–c) and noting $\chi \sim 1$ near the mush–liquid interface, this corresponds to a dimensionless distance of $\delta =O (Da^{1/ 2})$ from the interface. This displacement from the mush–liquid interface can be neglected whenever it is much smaller than the thickness of the convecting cells (requiring $Da^{1/ 2} \ll h$ for full-depth convection, or $Da^{1/ 2} \ll h\mathcal {C}$ if convection is confined to a fraction of order $\mathcal {C}$ of the full depth.

At the cooling boundary, we apply a linear thermal exchange and an impermeable condition

(2.6a,b)\begin{equation} \frac{\partial{\theta}}{\partial{z}} = {\mathcal{B}_i}\theta,\quad \boldsymbol{{e_z}}\boldsymbol{\cdot}\boldsymbol{u} = 0,\quad \text{when}\ {z = 0}, \end{equation}

where the Biot number ${\mathcal {B}_i} = \mathfrak{h}\hat {d}/k$ characterises the ratio of the effective thermal conductivity $\mathfrak {h}\hat {d}$ of the boundary vs the thermal conductivity $k$ of the mushy layer. The limit ${\mathcal {B}_i} \rightarrow \infty$ represents a perfectly conducting, isothermal, boundary (see Part 1 for further discussion).

We will focus on the development of perturbations in the mushy layer $0< z< h$ using (2.1a,b)–(2.6a,b). However, the background state also depends on properties of the liquid region which has no flow $(\boldsymbol {u}=0)$, uniform salinity $S=1$ and temperature evolving according to (2.2a) with $\chi =1$ and $\theta \rightarrow \theta _{\infty }$ as $z\rightarrow \infty$. Here, the dimensionless far-field temperature of the liquid $\theta _{\infty }={(T_\infty -T_c)}/{(T_{L\infty }-T_c)}$, with $\theta _{\infty }=1$ indicating fluid is at the liquidus temperature for salinity $\hat {S}_\infty$, and $\theta _{\infty }$ increases with the far-field liquid temperature. This parameter and its effects on the dynamics are discussed further in Part 1.

3. Stability method

We use a frozen time stability method (also known as a quasi-static stability analysis) to find the linear stability criteria of the mushy layer in various parameter configurations. This neglects nonlinear products of the perturbations to the base state, and assumes that the base state evolves at a much slower rate than the perturbation time scale, such that the base state can be considered static, or ‘frozen’.

Assuming a two-dimensional convection pattern, incompressibility (2.4b) is satisfied by writing $\boldsymbol {u} = \boldsymbol {\nabla }\times (\psi \boldsymbol {{e_y}})$ using a streamfunction $\psi$ and unit vector $\boldsymbol {{e_y}}$ perpendicular to the plane. This results in a system with three fields, which are expanded as base states and perturbations

(3.1a–c)\begin{equation} \theta = {\theta}_B(z,t) + \epsilon {\theta}_p(x,z,t),\quad \chi= {\chi}_B(z,t) + \epsilon {\chi}_p(x,z,t),\quad \psi = \epsilon \psi_p(x,z,t), \end{equation}

where $\epsilon$ is the small perturbation amplitude and the base-state streamfunction is uniformly zero. The base-state equations of $O(1)$ yielded by this substitution are those solved in Part 1 for diffusive growth, while the linear perturbation equations found by gathering terms of $O(\epsilon )$ are

(3.2a)\begin{gather} \frac{\partial{{\theta}_p}}{\partial{t}} + \frac{\partial{{\theta}_B}}{\partial{z}}\frac{\partial{\psi_p}}{\partial{x}} + {\mathcal{S}_t}\frac{\partial{{\chi}_p}}{\partial{t}} = \nabla^2{\theta}_p, \end{gather}

(3.2b)\begin{gather}-\mathcal{R}_{m}\frac{\partial{{\theta}_p}}{\partial{x}} = \boldsymbol{\nabla}\boldsymbol{\cdot}\left[\frac{1}{{{\chi}_B}^3}\boldsymbol{\nabla}\psi_p\right], \end{gather}

(3.2c)\begin{gather}(\mathcal{C} + 1-{\theta}_B)\frac{\partial{{\chi}_p}}{\partial{t}} - {\theta}_p\frac{\partial{{\chi}_B}}{\partial{t}} - {\chi}_p\frac{\partial{{\theta}_B}}{\partial{t}} - {\chi}_B\frac{\partial{{\theta}_p}}{\partial{t}} = \frac{\partial{{\theta}_B}}{\partial{z}}\frac{\partial{\psi_p}}{\partial{x}}, \end{gather}

where (3.2b) arises from rearranging Darcy's law (2.4a) and taking the curl. We decompose into Fourier modes with wavenumber $\mathrm {k}$ and exponential growth rate $\sigma$

(3.3)\begin{equation} \left[\begin{array}{c} {\theta}_p(x,z,t)\\ \psi_p(x,z,t)\\ {\chi}_p(x,z,t) \end{array}\right]= \left[\begin{array}{c} {\theta}_f(\mathrm{k},z)\\ {\rm i}\mathrm{k}\psi_f(\mathrm{k},z)\\ {\chi}_f(\mathrm{k},z) \end{array}\right] {\rm e}^{{\rm i}\mathrm {k}x}\,{\rm e}^{\sigma t}, \end{equation}

with (3.2) yielding a matrix eigenvalue problem for either $\sigma$ or $\mathcal {R}_{m}$

(3.4a)\begin{gather} {\textit{A}}\boldsymbol{y} - \mathcal{R}_{m} {\textit{B}}\boldsymbol{y} = \sigma{\textit{C}}\boldsymbol{y}, \end{gather}

(3.4b,c) \begin{gather}\boldsymbol{y} = \left[\begin{array}{@{}c@{}} {\theta}_f \\ \psi_f \\ {\chi}_f \end{array}\right],\quad {\textit{A}} =\left[\begin{array}{@{}ccc@{}} \partial^2/\partial z^2 - \mathrm{k}^2 & \mathrm{k}^2 {\partial {{\theta}_B}}/{\partial {z}} & 0\\ 0 & \dfrac{\partial{}}{\partial{z}}\left(\dfrac{1}{{{\chi}_B}^3}\dfrac{\partial{}}{\partial{z}}\right) - \dfrac{\mathrm{k}^2}{{{\chi}_B}^3} & 0\\ -{\partial {{\chi}_B}}/{\partial {t}} & \mathrm{k}^2{\partial {{\theta}_B}}/{\partial {z}} & -{\partial {{\theta}_B}}/{\partial {t}} \end{array}\right], \end{gather}

(3.4d,e) \begin{gather}{\textit{B}} = \left[\begin{array}{@{}ccc@{}} 0 & 0 & 0\\ -1 & 0 & 0\\ 0 & 0 & 0 \end{array}\right],\quad {\textit{C}} =\left[\begin{array}{@{}ccc@{}} 1 & 0 & {\mathcal{S}_t}\\ 0 & 0 & 0\\ {\chi}_B & 0 & -(\mathcal{C} + 1-{\theta}_B) \end{array}\right]. \end{gather}

For computational efficiency, the results below solve the problem directly for the critical Rayleigh number by assuming $\sigma = 0$ at the exchange of stability. Setting $\sigma = 0$ has the effect of making the stability independent of perturbations to the liquid fraction, because (3.4) decouples into a $2\times 2$ matrix eigenvalue problem for ${\theta }_f$ and $\psi _f$ with ${\chi }_f$ decoupled and determined from ${\theta }_f$ and $\psi _f$. However, this also explicitly neglects oscillatory modes from the calculation, with $\Im \{\sigma \} \neq 0$. Oscillatory modes contained within the mushy region have previously been found for models with asymptotically thin mushy layers in the near-eutectic limit (see Anderson & Worster Reference Anderson and Worster1996; Anderson & Guba Reference Anderson and Guba2020), but usually have a relatively minor impact on the critical Rayleigh number for instability. For example, the results illustrated by Anderson & Worster (Reference Anderson and Worster1996) show changes of ${\lesssim }20\,\%$ to the critical Rayleigh number for convective onset when oscillatory modes are the most unstable.

Applying the same decomposition to the boundary conditions leads to

(3.5a,b)\begin{gather} \frac{\partial{{\theta}_f}}{\partial{z}} = {\mathcal{B}_i}{\theta}_f,\quad \psi_f =0,\quad \text{at}\ {z=0}, \end{gather}

(3.5c,d)\begin{gather}{\theta}_f =0,\quad \frac{\partial{\psi_f}}{\partial{z}} =0,\quad \text{at}\ {z=h}, \end{gather}

where for simplicity we have neglected deformation of the mush–liquid interface and assumed no perturbations to $h$. In an analysis of asymptotically thin mushy layers in the near-eutectic limit, Roper, Davis & Voorhees (Reference Roper, Davis and Voorhees2008) showed this assumption does not qualitatively impact the stability results, and only leads to a modest quantitative difference in the critical $\mathcal {R}_{m}$ by ${\lesssim }10\,\%$.

Time now only appears in (3.4) and (3.5) through the base-state temperature and liquid fraction, ${\theta }_B$ and ${\chi }_B$, but we showed in Part 1 that these profiles can be calculated in self-similar coordinates using the variables $\tilde z$ and ${\tilde {\mathcal {B}}_i}$. We therefore repeat the same self-similarity transformation

(3.6a–c)\begin{gather} \tilde z = \frac{z}{\sqrt{t}},\quad \tilde{\mathrm{k}} = \mathrm{k} \sqrt{t},\quad \lambda = \frac{h}{\sqrt{t}}, \end{gather}

(3.6d–g)\begin{gather}\tilde{\mathcal{R}}_{m} =\mathcal{R}_{m} \sqrt{t},\quad {\tilde{\mathcal{B}}_i} = {\mathcal{B}_i} \sqrt{t},\quad \tilde{\theta} =\theta,\quad \tilde\psi_f = \frac{\psi_f}{\sqrt{t}}, \end{gather}

which we noted is equivalent to changing the length scale of the problem from the fixed length $\hat {d}$ to the time-evolving thermal diffusion length scale $\sqrt {\kappa \hat {t}}$. This transformation leaves the decoupled, $2\times 2$ form of (3.4) and (3.5) unchanged.

Finally, we note that varying parameters such as the far-field liquid temperature and Stefan number will change the depth of the mushy layer. This changing vertical extent will have a well-understood impact on the convective instability, with the critical Rayleigh number varying in inverse proportion to the layer depth. To allow us to see past this effect, it is useful to define the scaled Rayleigh number and the corresponding wavenumber

(3.7a,b)\begin{equation} \tilde{\mathcal{R}}_{m}^* = \lambda\tilde{\mathcal{R}}_{m} = \frac{g(\alpha+ {\beta}/{\varGamma})\Delta T\hat{\varPi}_0\hat{h}(\hat{t})}{\kappa\nu} = \frac{\hat{h}(\hat{t})}{\hat{d}_m},\quad \tilde{\mathrm{k}}^* = \tilde{\mathrm{k}} \lambda = \hat{\mathrm{k}} \hat{h}(\hat{t}), \end{equation}

where $\lambda =h(t)/\sqrt {t}$ is the self-similar interface position and the convective length scale $\hat {d}_m=\kappa \nu /[g(\alpha +{\beta }/{\varGamma })\Delta T \hat {\varPi }_0]$. In an experiment or application, the scaled Rayleigh number in (3.7a) will typically increase over time as the mushy layer grows and $\hat {h}$ increases. Our goal is to find the critical threshold value of $\tilde {\mathcal {R}}_{m}^*$ for the onset of instability as a function of the growth conditions that control the structure of the mushy layer, and also as a function of ${\tilde {\mathcal {B}}_i}$ which increases with time when there is imperfect boundary cooling. One can then infer a critical mush thickness or critical onset time for convection by inverting (3.7a) for $\hat {h}$ or $\hat {t}$.

The resulting eigenvalue system is discretised on a Chebyshev grid (see Trefethen Reference Trefethen2000) using between 60 and 200 nodes until the results are well converged. Following the method in Hitchen & Wells (Reference Hitchen and Wells2016) we rewrite the discretised boundary conditions as a constraint on the internal points. The critical $\mathcal {R}_{m}$ is found at a range of wavenumbers using the Matlab ‘eig’ routine, and minimised to find the global stability point for that set of parameters.

5. Effect of imperfect boundary conduction

With an imperfectly conducting boundary, the solution depends on a rescaled Biot number ${\tilde {\mathcal {B}}_i} = {\mathcal {B}_i} \sqrt {t}$ from (3.6) which also acts as a proxy for time. For deep porous layers of uniform porosity that are transiently cooled by an imperfectly conducting boundary, Hitchen & Wells (Reference Hitchen and Wells2016) found three ways that the Biot number could affect the convective stability: by changing the temperature range in the fluid, by changing the structure of the thermal forcing and by changing the boundary resistance applied to the thermal perturbation. For a mushy layer, the evolving Biot number will also affect the porosity of the mushy layer, with an initially high porosity throughout the mush depth, before porosities reduce as the Biot number ${\tilde {\mathcal {B}}_i} = {\mathcal {B}_i}\sqrt {t}$ grows (see Part 1). In this section, we investigate how these four effects impact the convective stability of a mushy layer, focussing on examples in both of the high-liquid-fraction and Stefan regimes.

For the high-liquid-fraction regime, we consider a system with $\mathcal {C} = 12.5$, $\theta _{\infty } = 1.25$ and ${\mathcal {S}_t} = 125$. Such a system initially freezes with a Biot number of ${\tilde {\mathcal {B}}_i} = 0.21 \equiv {\tilde {\mathcal {B}}_f}$ and adjusts to reach a self-similar state for ${\tilde {\mathcal {B}}_i}\gtrsim 100$ at late times (see figure 6 in Part 1). There are significant changes to the growth rate, surface temperature and structure of the mushy layer as the effects of the cooling propagate through the system. Figure 9 shows that the scaled ‘global’ Rayleigh number $\tilde {\mathcal {R}}_{m}^*$, which is defined using the temperature range between the mush–liquid interface and the cooling heat sink and is given by (3.7a), initially decreases sharply with increasing ${\tilde {\mathcal {B}}_i}$ (e.g. increasing time), before reaching a plateau as the mush reaches the self-similar state at large time. The critical wavenumber increases with the Biot number (and with time). With imperfect thermal conduction from the boundary, the surface temperature evolves over time and the temperature difference across the mushy layer grows as ${\tilde {\mathcal {B}}_i}$ increases (figure 9c). To see the effects of this changing temperature difference, and to isolate it from structural and perturbation boundary condition changes, figure 9(a) also shows the evolution of a ‘local’ Rayleigh number

(5.1)\begin{equation} \tilde{\mathcal{R}}_{{m,l}}^*= \frac{g(\alpha+{\beta}/{\varGamma}) [T_{L\infty}-T(\hat{z}=0,\hat{t})]\hat{\varPi}_0\hat{h}}{\kappa\nu}, \end{equation}

defined using the temperature difference $T_{L\infty }-T(\hat {z}=0,\hat {t})$ across the mushy layer only, rather than the temperature difference $\Delta T = T_{L\infty }-T_c$ between the mush–liquid interface and cooling heat sink. Here, hatted variables are dimensional. In contrast to the global Rayleigh number $\tilde {\mathcal {R}}_{m}^*$ (blue curve), the local Rayleigh number $\tilde {\mathcal {R}}_{{m,l}}^*$ (green curve) shows a modest increase with ${\tilde {\mathcal {B}}_i}$ and remains of similar order with $18.4\leq \tilde {\mathcal {R}}_{{m,l}}^* \leq 33.8$ for ${\tilde {\mathcal {B}}_f}\leq {\tilde {\mathcal {B}}_i} < \infty$. This indicates that the changing thermal range is the most significant effect of varying Biot number which increases over time. The global and local Rayleigh numbers converge for large ${\tilde {\mathcal {B}}_i}$ or late times, as the boundary approaches the temperature of the heat sink.

Figure 9. (a) Critical Rayleigh numbers, (b) critical wavenumbers and (c) temperature difference $1 - \theta (\tilde z =0,{\tilde {\mathcal {B}}_i})$ across the mushy layer against Biot number for $\mathcal {C} = 12.5$, $\theta _{\infty } = 1.25$ and ${\mathcal {S}_t} = 125$. These conditions place the system well into the high-liquid-fraction regime, with $\chi >0.92$ for all ${\tilde {\mathcal {B}}_i}$. Two different definitions for the Rayleigh number are used. In blue is the ‘global’ Rayleigh number, $\tilde {\mathcal {R}}_{m}^*$ defined in (3.7a) in terms of the temperature difference from the mush–liquid interface to the temperature of the cooling heat sink and using the mushy-layer depth as the length scale. The ‘local’ Rayleigh number, $\tilde {\mathcal {R}}_{{m,l}}^*$, in green is defined in (5.1) using the temperature difference $1 - \theta (\tilde z =0,{\tilde {\mathcal {B}}_i})$ across the mushy layer only. Shown in dashed grey are the first solidification at ${\tilde {\mathcal {B}}_i} = 0.21$ and the approximate transition from the adjustment period to the self-similar period at ${\tilde {\mathcal {B}}_i} \sim 100$ identified in Part 1.

It is interesting to consider how the other three factors affect the local Rayleigh number. The liquid fraction is uniformly one at first freezing and then decreases over time as the surface temperature adjusts and the surface solidifies (see figure 6 in Part 1). Decreased porosity at the surface decreases the permeability and increases the resistance to flow, providing a tendency for the critical Rayleigh number to increase. But these effects remain small in the high-liquid-fraction regime (here, the minimum liquid fraction $\chi =0.92$). Similarly, Hitchen & Wells (Reference Hitchen and Wells2016) showed that increasing the Biot number applied to a thermal perturbation leads to an increase in the critical Rayleigh number for thermal convection in a porous medium. This is because a larger Biot number increases the resistance from the boundary to changes in temperature, as seen by a number of previous studies such as Kubitschek & Weidman (Reference Kubitschek and Weidman2003) and Barletta, Tyvand & Nygård (Reference Barletta, Tyvand and Nygård2015). Finally, the thermal gradients are more uniform when the ice initially forms but at later times the temperature profile becomes more nonlinear with a sharper gradient near the surface and a shallower gradient at depth. Figure 7 demonstrates that a more nonlinear temperature profile tends towards an increased Rayleigh number, so changes to the thermal structure will also result in an increased Rayleigh number as the Biot number increases. Thus changes to the liquid-fraction structure, perturbation boundary condition and thermal structure all act to the increase the local Rayleigh number as the Biot number increases, as indicated by the green curve in figure 9(a). However, these changes are not significant enough to overcome the reduction in the global Rayleigh number caused by increased thermal forcing, as shown by the blue curve in figure 9(a).

We next consider the Stefan regime, which has $\mathcal {C}\ll 1$ and relatively low liquid fractions throughout much of the mush. Figure 10 shows the effects of imperfect conduction on convective stability in this regime. Initially the porosity is high throughout the layer (figure 10c) but a large low-porosity region begins to form for ${\tilde {\mathcal {B}}_i}\gtrsim 2$ and eventually occupies roughly $90\,\%$ of the mushy layer (see figure 9 of Part 1 for the full time evolution). In this instance, the porosity has a significant impact on stability due to convective confinement. While in figure 9(a) the local Rayleigh number $\tilde {\mathcal {R}}_{{m,l}}^*$ roughly doubled between first growth at small ${\tilde {\mathcal {B}}_i}$ and the self-similar limit for large ${\tilde {\mathcal {B}}_i}$, in figure 10(a) it increases over three orders of magnitude. To provide further insight we also calculate a global critical Rayleigh number using the full temperature range but applying a uniform porosity $\chi =1$ in the stability calculation, shown in orange in figure 10(a), which is initially close to the scaled Rayleigh number $\tilde {\mathcal {R}}_{m}^*$, but rapidly drops two and a half orders of magnitude. The critical wavenumber increases by an order of magnitude in figure 10(b) for calculations with variable porosity (blue curve), indicating a dramatic narrowing of the convective cells relative to the mushy-layer depth, but hardly changes for calculations with uniform porosity. By noting the preference for square aspect ratios in convective patterns, this is further evidence of the convective confinement becoming more significant as the Biot number increases.

Figure 10. (a) Critical Rayleigh numbers and (b) critical wavenumbers against Biot number for $\mathcal {C} = 0.017$, $\theta _{\infty } = 1.7$ and ${\mathcal {S}_t} = 500$ in the Stefan regime, with a large region of low porosity forming at late times. Two different definitions for the Rayleigh number are used. In blue is the ‘global’ Rayleigh number defined in (3.7a) in terms of the temperature difference from the mush–liquid interface to the cooling heat sink while the ‘local’ Rayleigh number in green is defined in (5.1) in terms of the temperature difference across the mushy layer only. The global critical Rayleigh number and wavenumber are also calculated with a uniform porosity, and shown in orange for comparison. Shown in dashed grey are the first solidification at ${\tilde {\mathcal {B}}_i} = 0.53$ and approximate transition from the adjustment period to the self-similar period at ${\tilde {\mathcal {B}}_i} \sim 1000$. Due to low-level noise in the interface position being amplified by the stability calculations, $\lambda$ was smoothed with a running average in ${\tilde {\mathcal {B}}_i}$ before being used in the stability analysis. The red-dashed curve shows an empirical fit $\tilde {\mathcal {R}}_{m}^* = 4000 + 1800[\ln ({{\tilde {\mathcal {B}}_i}}/{2.9})]^2$ around the minimum global $\tilde {\mathcal {R}}_{m}^*$, with numerical prefactors determined using a least-squares method. (c) The depth variation of porosity $\chi$ evolves over time as $\tilde{\mathcal{B}}_i$ changes, where depth $\tilde{z}/\lambda$ has been scaled by the mush thickness.

The above calculations suggest that the critical Rayleigh number $\tilde {\mathcal {R}}_{m}^*$ (blue curve in figure 10a) is being set by two competing effects. To illustrate this we highlight the convective profiles for three different Biot numbers in figure 11, which show the patterns before, during and after the minimum of the global Rayleigh number. As the Biot number increases over time, the increasing temperature difference across the mushy layer increases the thermal forcing and therefore increases the amount of potential energy available to drive convection. This acts to decrease the global critical Rayleigh number for onset of instability, and leads to the sharp initial decline seen for low values of ${\tilde {\mathcal {B}}_i}$ in figure 10(a) for the Stefan regime and also in figure 9(a) for the high-liquid-fraction regime. This regime is illustrated in figure 11(a) for ${\tilde {\mathcal {B}}_i}=1.0$ where there is very little energy available as $\theta (0)>0.99$, but there is also little resistance to flow with a minimum permeability of $\varPi =0.52$ and the convective cells occupy the whole domain to maximise the use of the little energy available. However, as the Biot number increases over time, the impact of the increasing temperature difference is later overcome as increased solidification reduces the permeability and increases the resistance to flow. Hence the growth of the low-porosity region significantly affects the vertical extent of the convection. We see this in figure 11(c) with ${\tilde {\mathcal {B}}_i}=10.0$ where the surface temperature has now dropped to $\theta (0)=0.84$, but roughly 70 % of the domain has a permeability $\varPi <0.1$. The reduced permeability prohibits utilisation of the potential energy in that region to drive convection, and hence convection is localised into the high-porosity interface region. A balance forms between these two effects which results in a minimum in the scaled Rayleigh number $\tilde {\mathcal {R}}_{m}^*$ at ${\tilde {\mathcal {B}}_i}\approx 2.81$, with corresponding surface values $\chi (0) = 0.42$ and $\theta (0) = 0.97$, with the background structure and perturbations illustrated in figure 11(b). The temperature difference across the mushy layer only represents $3\,\%$ of the potential total temperature difference, $T_{L\infty }-T_c$, at this point but the surface permeability has now dropped to $\varPi (0)=0.07$. This is starting to provide appreciable resistance to flow at the surface, and comparing with (a) we can see that the convection has shifted slightly downwards. However, the convection is still reaching the majority of the domain and so able to use almost all of the potential energy available, thus creating the optimum conditions for convection and resulting in the lowest critical Rayleigh number.

Figure 11. Background profiles of liquid fraction, permeability and temperature (left column) and convective perturbations (right column) illustrating three different regimes from figure 10 for (a) ${\tilde {\mathcal {B}}_i} = 1.0$ (descending branch), before the minimum Rayleigh number, (b) ${\tilde {\mathcal {B}}_i} = 2.9$, approximately aligned with the minimum Rayleigh number, (c) ${\tilde {\mathcal {B}}_i} = 10.0$ (ascending branch), after the minimum Rayleigh number but before the long-time behaviour. The other parameter values $\mathcal {C}= 0.017$, $\theta _\infty = 1.7$ and ${\mathcal {S}_t} = 500$ are the same as figure 10. In the right column, the colour scale shows the thermal perturbation ${\theta }_p$ and streamfunction contours are shown at $0.02$, $0.05$, $0.1$ and $0.2$, with black contours representing clockwise circulation and grey contours representing anti-clockwise circulation.

6. Application to sea-ice growth

We apply the results to a case study of sea-ice growth in Arctic conditions and consider conditions relevant to the field experiment in Notz & Worster (Reference Notz and Worster2008), for initial ice growth in a hole cut into land-fast sea ice in Svalbard. Sea water with salinity $\hat {S}_\infty =35\ {\rm g}\ {\rm kg}^{-1}$ was cooled from $T_\infty = -1\,^{\circ }$C with atmospheric temperature $T_c = -30\,^{\circ }$C in the boundary layer, which gives the dimensionless parameters $\mathcal {C} = 0.11$, $\theta _{\infty } = 1.1$ and ${\mathcal {S}_t} = 6.2$, using $L/c_p=168\,^{\circ }$C appropriate to ice properties, $\hat {S}_s=0$ and $\varGamma =0.085\ \mathrm {K}\ (\mathrm {g}\ {\rm kg}^{-1})^{-1}$. We also use a permeability reference $\hat {\varPi }_0 = 2\times 10^{-8}\ {\rm m}^2$ (Freitag Reference Freitag1999), $g=9.81\ \mathrm {m}\ \mathrm {s}^{-2}$, $\alpha =4\times 10^{-5}\ \mathrm {K}^{-1}$, $\beta =8 \times 10^{-4}\,(\mathrm {g}\ \mathrm {kg}^{-1})^{-1}$, $\kappa =1.2 \times 10^{-6}\ \mathrm {m}^2\ \mathrm {s}^{-1}$, $\nu =1.8 \times 10^{-6}\ \mathrm {m}^2\ \mathrm {s}^{-1}$ and $k=2.2\ \mathrm {W}\ \mathrm {m}^{-1}\ \mathrm {K}^{-1}$ (see chapter 1 of Hitchen (Reference Hitchen2017) for discussion of these estimates).

Figure 12(a) plots the critical depth for onset of convection in this set-up against Biot number, found by rearranging (3.7a) to give $\hat {h}_c$. While this set-up is in the transitional regime rather than the pure Stefan regime, we still observe the non-monotonic behaviour in the scaled Rayleigh number and critical depth. In this instance the minimum scaled Rayleigh number is $550$, corresponding to a critical depth of 2.4 cm. This occurs when ${\tilde {\mathcal {B}}_i} = 0.81$ and just after half of the surface layer has solidified, $\chi (0,t) = 0.43$, with a dimensionless surface temperature of $0.86$, which corresponds to $-6.9\,^{\circ }$C. The scaled Rayleigh number then rises to a long-time limit of $1165$, which is equivalent to a depth of 5.1 cm for onset of convection at large ${\tilde {\mathcal {B}}_i}$ or long times.

Figure 12. Convective stability characteristics for mush growth in sea water of salinity $35\ \mathrm {g}\ \mathrm {kg}^{-1}$ cooled from $-1\,^{\circ }$C to $-30\,^{\circ }$C, with all other properties described in the text. (a) The critical depth $\hat {h}_c$ (black curve) and time-evolving ice depth $\hat {h}$ (coloured curves) against self-similar Biot number ${\tilde {\mathcal {B}}_i}={\mathfrak {h}\sqrt {\kappa \hat {t}}}/{k}$ for four different surface wind speeds, using $\mathfrak {h}=\rho_ac _{pa}{C_s}u_w$ for cooling by parameterised turbulent atmospheric heat fluxes. Shown on the alternate vertical scale are the scaled Rayleigh numbers corresponding to these ice depths from (3.7a). The ice depth monotonically increases with Biot number, and hence time, but the critical depth exhibits the non-monotonic variation first observed in figure 10 for the Stefan regime which arises from the formation of a low-porosity region and subsequent convective confinement. (b) The onset time and (c) onset depth (blue) and half-wavelength (red) against wind speed. Faster wind speeds result in convective onset occurring sooner, at shorter wavelengths and with thinner ice for $u_w < 8.6\ {\rm m}\ {\rm s}^{-1}$, but thereafter ice thickness increases slightly with wind speed.

The time-evolving ice-thickness trajectories depend on ${\tilde {\mathcal {B}}_i}$ and hence the efficiency of cooling. We consider a case where the dominant heat losses are due to turbulent sensible heat fluxes parameterised in terms of a 10-metre wind speed $u_w$ by $k \partial T/\partial \hat {z} = \rho_{a}c_{pa}C_s u _w(T - T_c)$ at $\hat {z} = 0$, where the air density $\rho _a=1.3\ \mathrm {kg}\ \mathrm {m}^{-3}$, air heat capacity $c_{pa}=10^3\ \mathrm {J}\ \mathrm {kg}^{-1}\ \mathrm {K}^{-1}$ and heat transfer coefficient $C_s=2.5 \times 10^{-3}$. Increasing the wind speed enhances the efficiency of heat transfer and increases ${\tilde {\mathcal {B}}_i}$. Figure 12(a) shows the resulting time-evolving trajectories of ice depth for four different surface wind speeds of similar magnitude to those observed by Notz & Worster (Reference Notz and Worster2008). For $u_w \leq 8.6\ \mathrm {m}\ \mathrm {s}^{-1}$, increasing the wind speed decreases the depth at which the mushy layer becomes convectively unstable. This can be seen by the critical depth and the physical growth trajectories intercepting at a smaller thickness. A wind speed of $8.6\ \mathrm {m}\ \mathrm {s}^{-1}$ represents the ‘optimal’ rate of cooling, where the mushy layer becomes linearly unstable with the least possible ice growth. Increasing the wind speed further delays the convective onset until the ice has reached a greater depth.

Calculating the Biot number at which $\hat {h}$ intercepts the stability curve $\hat {h}_c$ allows one to infer the onset time, $\hat {t}_i$, and the wavelength of convection ${2{\rm \pi} \hat {h}_c}/{\mathrm {k}_{c}^*}$. These are shown in figure 12(b,c) for a range of realistic wind speeds. Whilst an increase in the wind speed increases the value of ${\tilde {\mathcal {B}}_i}$ needed for convective onset, the actual onset time decreases as the wind speed increases. This occurs because higher wind speeds allow faster heat transfer from the surface, which enhances the cooling of the mushy layer. As a result, the mushy layer grows faster and reaches the critical depth sooner. The critical depth is also lowered by the increasing wind speed for $u_w< 8.6\ {\rm m}\ {\rm s}^{-1}$. Increasing the wind speed also decreases the wavelength of convective patterns which form, with similar order of magnitude to the decreasing critical depth.

The observations of Notz & Worster (Reference Notz and Worster2008) show a noticeable decrease in the bulk salinity at the first measurement point, indicating that convection has started sometime during the first six hours of ice growth for an ice depth $\hat {h}\leq 4$ cm. Using figure 12 and an average observed wind speed of $2\ {\rm m}\ {\rm s}^{-1}$ yields convective onset 170 min after the start of the experiment with a 3.1 cm thick mush which lie within the observational bounds. The corresponding linearly unstable wavelength is 9.5 cm.

Finally, we recall that in § 2 we assumed that the ratio $Da^{1/ 2}/h$ is small vs the minimum of $1$ or $\mathcal {C}$, in order to approximately apply the flow boundary conditions at the mush–liquid interface rather than displaced appreciably into the high-porosity mush near the interface. Using $\hat {\varPi }_0 = 2\times 10^{-8}\ {\rm m}^2$ and the critical ice thickness for convection of $0.03\ \mathrm {m}$ in this case study yields $Da^{1/ 2}/h \approx 5 \times 10^{-3}$ which is smaller than $\mathcal {C}\approx 0.1$ and justifies the assumption made previously.

7. Summary and conclusions

This study has used a linear stability analysis to consider the impact of temperature and liquid fraction profiles on the onset of convection in a transiently growing mushy layer, from either a perfectly cooled isothermal boundary or an imperfectly cooled boundary. We focussed on convective modes driven within the mushy layer, with the neighbouring liquid treated as maintaining an isothermal, constant-pressure boundary at the mush–liquid interface. We determined the scaled Rayleigh number $\tilde {\mathcal {R}}_{m}^* = {g(\alpha +{\beta }/{\varGamma })(T_{L\infty }-T_c)\hat {\varPi }_0\hat {h}(\hat {t})}/{\kappa \nu }$ for linear instability, corresponding critical depth $\hat {h}_c$ and scaled wavenumber $\tilde {\mathrm {k}}^* = \hat {\mathrm {k}}\hat {h}(\hat {t})$, as a function of the salinity and cooling conditions described by dimensionless parameters $\mathcal {C}$, $\theta _{\infty }$, ${\mathcal {S}_t}$ and ${\tilde {\mathcal {B}}_i}$.

In Part 1 we identified that the internal mush structure shows different limiting behaviour depending on the concentration ratio $\mathcal {C} ={(\hat {S}_\infty -\hat {S}_s)}/{(\hat {S}_c-\hat {S}_\infty )} =\varGamma (\hat {S}_\infty -\hat {S}_s)/(T_{L\infty } - T_c)$, here determined by the freezing-point depression vs the maximal temperature difference across the mush. We here find a corresponding impact on convective stability for cooling by an isothermal perfectly conducting boundary. The scaled critical Rayleigh number approximately satisfies $\tilde {\mathcal {R}}_{m}^* \approx 39 + {24}/{\mathcal {C}^2}$ and the scaled critical wavenumber $\tilde {\mathrm {k}}^* \approx 2.4+{0.34}/{\mathcal {C}}$, with comparatively weaker dependence on the Stefan number ${\mathcal {S}_t}$ which measures the ratio of latent to specific heat. A qualitatively similar pattern of increase had been observed previously by Hwang (Reference Hwang2013), but the scaling had not been characterised quantitatively. The behaviour can be reconciled with a constant local Rayleigh number evaluated across the depth of the high-porosity region, and wavelength proportional to this depth. When $\mathcal {C}\gg 1$, the freezing-point depression is much larger than the range of temperatures in the system. The liquid fraction is close to one throughout and convection can penetrate the full depth of the mushy layer, with a local Rayleigh number based on the temperature difference across the mush yielding $\tilde {\mathcal {R}}_{m}^* \approx \mathrm {constant}$. However, the system is comparatively fresh and a large low-porosity region forms for the Stefan regime with $\mathcal {C}\ll 1$, which severely increases the resistance to flow. Convection is instead confined to the high-porosity interfacial region, which has a fractional thickness proportional to $\mathcal {C}$ and where the temperature changes by a fraction $\mathcal {C}$ of the total temperature difference between the boundary and mush–liquid interface. A constant local Rayleigh number across this regions yields a scaling $\mathcal {C}^2 \tilde {\mathcal {R}}_{m}^* \approx \mathrm {constant}$ for $\mathcal {C}\ll 1$. This confinement also explains the increase in wavenumber $\tilde {\mathrm {k}}^*$ for small $\mathcal {C}$, since convection has a preference for square aspect ratios cells so that the wavelength $2{\rm \pi} /\tilde {\mathrm {k}}^* \propto \lambda \mathcal {C}$ for $\mathcal {C}\ll 1$. The vertical extent of the flow is being reduced, so the convective cells narrow.

The effect of the Stefan number on convective instability depends on the concentration ratio. For $\mathcal {C}\gg 1$ the mushy layer has a high liquid fraction throughout its depth and there is a modest increase in the critical Rayleigh number as ${\mathcal {S}_t}$ increases. This is consistent with an interaction of the convective perturbation with resistance to flow from its boundary condition; increasing the Stefan number focusses the thermal forcing nearer to the restrictive no-normal-flow surface condition and away from the liberal constant-pressure interface condition. In contrast, there is a decrease in the critical Rayleigh number with increasing ${\mathcal {S}_t}$ in the Stefan regime when $\mathcal {C}\ll 1$. The largest effect comes from the impact of ${\mathcal {S}_t}$ on the liquid-fraction profiles, because increasing the Stefan number increases the size of the high-porosity region close to the interface. This promotes convection by reducing the resistance to flow. This effect is much larger than that caused by the changing temperature profiles.

For growth with an imperfectly conducting boundary, the dominant impact on the convective instability is due to the impact of imperfect cooling on the temperature difference across the mushy layer and the porosity structure. Changes to the internal temperature structure or arising from the imperfectly conducting boundary modifying the perturbation of surface temperature are both less significant. Increasing the temperature difference across the mushy layer directly increases the amount of potential energy available to drive convection and leads to a corresponding decrease in the critical $\tilde {\mathcal {R}}_{m}^*$. In the high-liquid-fraction regime, the change in liquid fraction as ${\tilde {\mathcal {B}}_i}$ is varied is comparatively small, and the temperature difference between the surface layer and the interface is the main driver of change to $\tilde {\mathcal {R}}_{m}^*$. However, in the Stefan regime the imperfectly conducting boundary leads to significant evolution of the liquid fraction. There is initially high porosity throughout the mushy layer, but at later times low-porosity regions develop which greatly increase the resistance to flow. This acts in competition with the total amount of available potential energy increasing over time, due to both increasing mushy-layer thickness and increasing temperature difference across the mushy layer as the surface temperature adjusts. As a result of this competition, $\tilde {\mathcal {R}}_{m}^*$ shows non-monotonic behaviour as ${\tilde {\mathcal {B}}_i}$ increases with $\mathcal {C}\ll 1$. Initially $\tilde {\mathcal {R}}_{m}^*$ decays, indicating reduced stability as potential energy accumulates in the initial high-porosity layer. But eventually $\tilde {\mathcal {R}}_{m}^*$ increases, with enhanced stability as a low-porosity zone develops and increases the resistance to flow. This is in contrast to the monotonic variation of $\tilde {\mathcal {R}}_{m}^*$ with ${\tilde {\mathcal {B}}_i}$ when $\mathcal {C}\gg 1$ and internal solidification and the corresponding permeability changes are relatively weak. The non-monotonic behaviour of $\tilde {\mathcal {R}}_{m}^*$ (and hence $\hat {h}_c$) combined with the steady growth of the physical ice depth, $\hat {h}(\hat {t})$ raises the possibility that certain conditions could exhibit a convective shutdown and restart. This could occur if the critical depth increases faster than the physical ice depth after an initial onset, but further investigation is needed to determine if this is possible.

For sea-ice growth driven by cooling from parameterised turbulent atmospheric heat transfer we find that the convective onset time, ice thickness and wavelength all initially decrease as the wind speed and efficiency of turbulent cooling increases. The ice thickness at the onset of convection then increases above a critical wind speed. Onset times are several hours, with ice depths of 2–8 cm for winter sea-ice growth in Arctic conditions with an atmospheric temperature of $-30\,^{\circ }$C. However, there will be quantitative uncertainties propagating through these estimates because there are challenges in accurately determining the sea-ice permeability for highly porous regions of ice (see discussion in Wells et al. Reference Wells, Hitchen and Parkinson2019). Such uncertainty in the permeability may also impact whether there may be non-trivial coupling between flow in pure liquid and porous regions in a transition zone near the mush–liquid interface (Le Bars & Worster Reference Le Bars and Worster2006).

The impact of porosity localisation on convective onset may have implications in several fields. Convection is a primary mechanism for the desalination of sea ice, controlling salt fluxes to the ocean, the material properties of the evolving ice and nutrient transport for biogeochemical processes (Notz & Worster Reference Notz and Worster2008; Wells, Wettlaufer & Orszag Reference Wells, Wettlaufer and Orszag2011; Vancoppenolle et al. Reference Vancoppenolle2013). Several parameterisations activate convection in regions where a local depth-dependent Rayleigh number $g \Delta \rho (z) (\hat {h}-\hat {z})\bar {\varPi }(z)/\kappa \mu$ exceeds a critical value, for an appropriately averaged permeability $\bar {\varPi }(z)$ (Worster & Rees Jones Reference Worster and Rees Jones2015). Our linear stability bounds for $\mathcal {C}\ll 1$ lend support to this approach, as they suggest convective localisation and a stability bound based off the properties across high-porosity regions where the permeability $\bar {\varPi }(z)$ is relatively large. Meanwhile in metallurgy, the onset of convection can lead to the formation of freckle defects from the generation of chimney features via the nonlinear development of convection and impact on porosity (Anderson & Guba Reference Anderson and Guba2020). Our results suggest that material grown with $\mathcal {C}\ll 1$ may restrict the onset of convection to a comparatively narrow depth range near the mush–liquid interface, although a fully nonlinear analysis is required to infer the quantitative consequences for development of chimneys and freckles.

The present analysis neglected certain effects which will be important in some settings. We focussed on mushy-layer-mode convection, neglecting flow generated by buoyancy forcing in the liquid regions. We have thus neglected the short wavelength boundary-layer mode of convection which can be significant in certain parameter regimes (Worster Reference Worster1992). It would also be interesting to examine the interaction of shear-driven flows and turbulent ocean mixing with the evolving convective confinement considered here. Perhaps most significantly, we have focussed on a linear analysis to provide physical insight into the localisation of convective flow, but have neglected nonlinear effects. These include nonlinear interactions with deformation of the mush–liquid interface, nonlinear advective feedbacks on convective flow and nonlinear evolution of the porosity leading to the formation of chimneys. Many studies of the nonlinear dynamics of chimney convection have focussed on the parameter regime $\mathcal {C}\gg 1$ (see Anderson & Guba Reference Anderson and Guba2020), but it would also be interesting to examine the impact of the convective localisation discussed here on nonlinear convection in mushy layers in the Stefan regime with $\mathcal {C}\ll 1$.