2.1 Mass balance for salt

The strong form of the mass balance equation for salt is

(2.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}u(t,x)=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}_{u},\quad \text{for }t\geqslant 0,x\in \unicode[STIX]{x1D6FA}(t),\end{eqnarray}$$

where $\boldsymbol{J}_{u}$ is the flux vector. Similar to Kutschan, Morawetz & Gemming (Reference Kutschan, Morawetz and Gemming2010), Thoms et al. (Reference Thoms, Kutschan, Morawetz and Cemming2014), we neglect the expansion effects due to phase change and assume further that the fluid is at rest. As with the spatially isothermal assumption, this simplification can be justified strongly in 1D and 2D for pipes and slits, but requires further justifications in 3D. With Fickian diffusion, the salt flux thus takes the form

(2.4) $$\begin{eqnarray}\displaystyle \boldsymbol{J}_{u}(t,x)=-D\unicode[STIX]{x1D735}u(t,x),\quad \text{for }t\geqslant 0,x\in \unicode[STIX]{x1D6FA}(t), & & \displaystyle\end{eqnarray}$$

where $D$ is the diffusion coefficient, here taken as a constant. With this, equation (2.3) becomes

(2.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}u(t,x)=D\unicode[STIX]{x0394}u(t,x),\quad \text{for }t\geqslant 0,x\in \unicode[STIX]{x1D6FA}(t).\end{eqnarray}$$

2.2 Salt expulsion at the boundary

Since the ice is assumed salt free, salt is expelled from the ice formed at the brine–ice interfaces. In order to conserve the mass of salt at the brine–ice interface, the normal component of the salt-expulsion flux $\boldsymbol{J}_{bc}\boldsymbol{\cdot }\boldsymbol{n}$ and normal component of the diffusive flux $\boldsymbol{J}_{u}$ must equal

(2.6) $$\begin{eqnarray}-\boldsymbol{J}_{bc}\boldsymbol{\cdot }\boldsymbol{n}+\boldsymbol{J}_{u}\boldsymbol{\cdot }\boldsymbol{n}=0,\quad \text{for }t\geqslant 0\text{ and }x\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}(t).\end{eqnarray}$$

Here $\boldsymbol{n}$ is the normal vector at the boundary $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ pointing inside the fluid. Letting $\boldsymbol{s}$ denote the position of the ice boundary, the conservation of salt as seen from the ice domain implies that the expulsion flux is given by

(2.7) $$\begin{eqnarray}\boldsymbol{J}_{bc}\equiv u_{crit}\frac{\text{d}\boldsymbol{s}}{\text{d}t}.\end{eqnarray}$$

Together equations (2.6) and (2.7) become a Stefan condition at the brine–ice boundary

(2.8) $$\begin{eqnarray}-u_{crit}\frac{\text{d}\boldsymbol{s}}{\text{d}t}\boldsymbol{\cdot }\boldsymbol{n}+\boldsymbol{J}_{u}\boldsymbol{\cdot }\boldsymbol{n}=0,\quad \text{for }t\geqslant 0\text{ and }x\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}(t).\end{eqnarray}$$

Figure 1 is a 1D illustration of the advancement $\text{d}s$ of the ice–water boundary over an infinitesimal time $\text{d}t$ , due to the increase of the critical salinity at the boundary from $u_{crit}$ to $u_{crit}+\text{d}u_{crit}$ . This results in the salt-expulsion flux $J_{bc}$ at the boundary and the corresponding change in the overall salt profile from solid to dashed line as a result of the diffusion in the water domain.

Figure 1. A sketch of the ice–water boundary movement caused by an increase of the critical salinity during freezing.

2.3 Nucleation

Salinity and phase transition are two processes that are strongly connected. As seen from the discussion above, a decrease in the temperature can lead to freezing associated with existing ice regions through an increase of the critical salinity $u_{crit}$ . A second important mechanism is the nucleation of new ice regions. The exact laws of nucleation are difficult to obtain. We expect supercooling to occur in the system that would trigger the nucleation process when below a (salinity dependent) nucleation temperature. Conversely, it can be stated that nucleation only happens below a (temperature dependent) nucleation salinity. To be consistent with the treatment of freezing, we adopt this latter viewpoint and thus we write $u_{nucl}=u_{nucl}(T(t))=u_{nucl}(t)$ . This allows us to naturally define a dimensionless number associated with nucleation,

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D707}(t)=\frac{u_{nucl}(t)}{u_{crit}(t)},\end{eqnarray}$$

where $0\leqslant \unicode[STIX]{x1D707}<1$ has the interpretation of a nucleation multiplier. Once nucleation is encountered, salinity decreases locally due to ice formation (as discussed above), and a measurable region of ice is formed at nucleation sites.

In general, the nucleation multiplier $\unicode[STIX]{x1D707}$ may vary with salinity. For the mathematical analysis in the next section, it is convenient to consider the lowest-order approximation of a constant multiplier. The nucleation multiplier $\unicode[STIX]{x1D707}$ can also be interpreted as a dimensionless parameter related to the degree of super-saturation of micro-crystals permitted in the fluid for the current salinity, and is thus related to the free energy released for a nucleation event. The interpretation of a constant nucleation multiplier $\unicode[STIX]{x1D707}$ is thus in principle an approximation of a nucleation probability dependent on the deviation in free energy from the equilibrium state (Mullin Reference Mullin2001). One should note that the first-order approximation of (2.9) around the regime of interest might not be extended to arbitrary values of concentration; our model for instance does not capture undercooling of freshwater.

2.4 The complete mathematical model

In the continuation, we will always consider the space, the time and the concentration dimensionless. More precisely, we use the natural non-dimensionalizations $x/x_{\ast }\rightarrow x$ , $u_{crit}(t)/u_{crit}(0)\rightarrow u_{crit}(t)$ and $t/t_{\ast }\rightarrow t$ where $t_{\ast }=(x_{\ast })^{2}/D$ . Here $x_{\ast }$ is the diameter of the domain (length, in 1D) and $t_{\ast }$ a characteristic time. By this, the non-dimensionalized $u_{crit}(0)=1$ and the dimensionless diffusion coefficient becomes 1.

The mathematical model discussed so far has the salinity $u$ and the location of the ice–brine boundary $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ as unknowns. The process is governed by the critical salinity at the ice–brine interface $u_{crit}$ , which is temperature and thus time dependent and defines the boundary conditions at $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}(t)$ .

(2.10) $$\begin{eqnarray}\left.\begin{array}{@{}ll@{}}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}u(t,x)=\unicode[STIX]{x0394}u(t,x), & \text{for }t\geqslant 0\text{ and }x\in \unicode[STIX]{x1D6FA};\\ u(t,x)\geqslant u_{nucl}(t) & \text{for }t\geqslant 0\text{ and }x\in \unicode[STIX]{x1D6FA};\\ u(t,x)=0 & \text{for }t\geqslant 0\text{ and }x\notin \unicode[STIX]{x1D6FA};\\ \displaystyle u_{crit}(t)\frac{\text{d}}{\text{d}t}\boldsymbol{s}(t)\boldsymbol{\cdot }\boldsymbol{n}=(-\unicode[STIX]{x1D735}u(t,x))\boldsymbol{\cdot }\boldsymbol{n}, & \text{for }t\geqslant 0\text{ and }x\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA};\\ u(t,x)=u_{crit}(t) & \text{for }t\geqslant 0\text{ and }x\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA};\\ u_{nucl}(t)=\unicode[STIX]{x1D707}u_{crit}(t), & \text{for }t\geqslant 0.\end{array}\right\}\end{eqnarray}$$

The system is completed by initial conditions.

2.5 Sherwood number

In the absence of spontaneous nucleation, the qualitative behaviour of the ice–water system can be characterized by a relation between the rate of the normal components of two fluxes at the moving ice–brine boundary: the salt-expulsion flux (2.7) and the diffusive flux (2.4). We let $\ell$ be a time-dependent characteristic domain size that can be seen as an average distance between ice boundaries. Observe that the entire process considered here is controlled by the critical salinity $u_{crit}$ and this can change rapidly in time, inducing different ice-formation regimes. Consequently, the characteristic domain size $\ell$ may change significantly in time.

We interpret $K$ as a characteristic normal velocity of the ice–brine boundary (to be defined precisely below), which scales as

(2.11) $$\begin{eqnarray}K\sim \frac{\text{d}\boldsymbol{s}}{\text{d}t}\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$

The ratio of mass-transfer rate to the diffusion rate is essential for this problem and is characterized by the Sherwood number (Cussler Reference Cussler2009), defined as (recall that the non-dimensionalized diffusion constant is unity)

(2.12) $$\begin{eqnarray}Sh=\frac{K}{\ell ^{-1}}.\end{eqnarray}$$

We emphasize that $Sh$ may change in time as it depends on both $\ell$ and $K$ , thus it is in principle a function. By a slight abuse of language, we will nevertheless refer to $Sh$ as a ‘number’. The case where $Sh$ is constant in time will be of particular importance later.

To make concrete the choice of $K$ , recall that the evolution of the ice–water interface is not known a priori but represents an unknown in the model. This evolution is controlled by $u_{crit}$ , which is a time-dependent input parameter of the model. As mentioned, $K$ itself depends on time through $\ell$ , which in its turn depends on $u_{crit}$ . Therefore it makes sense to express $K$ in terms of critical salinity. To do so we first refer to figure 1 sketching the change in the ice–brine system encountered within an infinitesimal time $\text{d}t$ . More precisely, as suggested by the area of fresh ice, in the normal direction the ice front is moving into the brine over the infinitesimal distance $\text{d}\boldsymbol{s}\boldsymbol{\cdot }\boldsymbol{n}$ , causing salt expulsion into the brine. The corresponding normal salt flux $\boldsymbol{J}_{bc}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}t$ can be obtained from (2.7). Consequently, the salt is distributed into the brine domain sized $\ell$ (the left part of 1), causing an increase in the salinity by $\text{d}\,u_{crit}$ (the grey dashed curve). By mass conservation, one gets

(2.13) $$\begin{eqnarray}\boldsymbol{J}_{bc}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}t=u_{crit}\,\text{d}\boldsymbol{s}\boldsymbol{\cdot }\boldsymbol{n}\sim \ell \left(\frac{\text{d}}{\text{d}t}u_{crit}\right)\text{d}t.\end{eqnarray}$$

Dividing in the above by $u_{crit}\,\text{d}t$ , from (2.11) we obtain the definition of the characteristic velocity of the system

(2.14) $$\begin{eqnarray}\displaystyle K=\frac{\ell }{u_{crit}}\frac{\text{d}}{\text{d}t}u_{crit}=\ell \frac{\text{d}}{\text{d}t}\ln u_{crit}. & & \displaystyle\end{eqnarray}$$

Combining the definition of the characteristic normal velocity with the Sherwood number in (2.12), we obtain

(2.15) $$\begin{eqnarray}Sh=\ell ^{2}\frac{\text{d}}{\text{d}t}\ln u_{crit}.\end{eqnarray}$$

The Sherwood number as defined above will play a critical role in our subsequent analysis.

In addition to the Sherwood number, the model under discussion is also influenced by the nucleation threshold. This is accounted for through the nucleation multiplier $\unicode[STIX]{x1D707}$ defined by equation (2.9), which is already a dimensionless parameter. In § 4 we show that the qualitative behaviour of the system can be characterized based exclusively on two dimensionless parameters: the Sherwood number $Sh$ and the nucleation multiplier  $\unicode[STIX]{x1D707}$ .

2.6 The dimensionless variables in 1D

While the mathematical model defined above is stated for arbitrary dimensions, the analysis in the remainder of this manuscript is restricted to the case of one spatial dimension. For concreteness, we reformulate the variables of (2.10) accordingly, by employing the characteristic quantities discussed in § 2.5.

Since in 1D the brine domain is a union of sub-intervals, we fix $t>0$ and for any brine interval we let $x_{0}(t)$ denote its left boundary and $\ell (t)$ its length. This allows us to rescale the spatial coordinate $x$ and introduce a local coordinate $y$ on each brine interval

(2.16) $$\begin{eqnarray}y=(x-x_{0}(t))/\ell (t).\end{eqnarray}$$

In this way, for each sub-interval of brine the solution is defined on the unit interval $y\in [0,1]$ .

Now note that the freezing and nucleation given in (2.10) ensures that for the salinity $u$ , one has

(2.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}u_{crit}(t)<u\leqslant u_{crit}(t), & & \displaystyle\end{eqnarray}$$

for any $t>0$ and $x$ in the brine domain. It therefore makes sense to introduce a time-dependent normalization by the salinity $u_{crit}$ , such that the dimensionless salinity is constrained to the time-independent interval $(\unicode[STIX]{x1D707},1]$ . This allows us to define a dimensionless salinity $\unicode[STIX]{x1D708}$ for a sub-domain through

(2.18a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D708}:(0,\infty )\times (0,1)\rightarrow (\unicode[STIX]{x1D707},1],\quad \unicode[STIX]{x1D708}(t,y)=\unicode[STIX]{x1D708}(t,(x-x_{0}(t))/\ell (t))=\frac{u(t,x)}{u_{crit}(t)}.\end{eqnarray}$$

Note that the time is dimensionless, as rescaled in § 2.4. Also, since in the dimensionless form $y=0$ and $y=1$ are ice–brine boundaries and due to the scaling of the salinity one has $\unicode[STIX]{x1D708}(0)=\unicode[STIX]{x1D708}(1)=1$ . By (2.9), $\unicode[STIX]{x1D708}$ cannot decrease below $\unicode[STIX]{x1D707}$ . In other words, whenever the minimum of $\unicode[STIX]{x1D708}$ is close to $\unicode[STIX]{x1D707}$ , the system is close to nucleation.

3.1 The critical case

We define the critical case to be the state where the salinity profile is self-similar and thus represents a case where no nucleation events occur. In other words, diffusion is in balance with the temperature decrease, representing a bifurcation point for the system. Recall that the dimensionless solution function $\unicode[STIX]{x1D708}$ , as defined in (2.18), is in this case constant in time

(3.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}(t,y)}{\unicode[STIX]{x2202}t}=0.\end{eqnarray}$$

To simplify the notation, since $\unicode[STIX]{x1D708}$ does not depend on time, we will in this section use the notation

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{1}(y)\equiv \unicode[STIX]{x1D708}(t,y).\end{eqnarray}$$

To find the boundary conditions leading to the critical case we consider $\unicode[STIX]{x1D708}_{1}$ as a self-similar solution of the original problem, in the spirit of Barenblatt (Reference Barenblatt1996). Using (2.16) and (2.5) together with the chain rule, we obtain

(3.3) $$\begin{eqnarray}\displaystyle 0=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}_{1}(y)}{\unicode[STIX]{x2202}t}\equiv \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}(t,y)}{\unicode[STIX]{x2202}t} & = & \displaystyle \frac{1}{u_{crit}^{2}}\left[\frac{\unicode[STIX]{x2202}u(t,x)}{\unicode[STIX]{x2202}t}u_{crit}-\frac{\text{d}u_{crit}(t)}{\text{d}t}u\right]\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{1}{u_{crit}^{2}}\left[\frac{\unicode[STIX]{x2202}^{2}u(t,x)}{\unicode[STIX]{x2202}x^{2}}u_{crit}(t)-\frac{\text{d}u_{crit}(t)}{\text{d}t}u(t,x)\right]\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{1}{\ell ^{2}(t)}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D708}_{1}(y)}{\unicode[STIX]{x2202}y^{2}}-\frac{1}{u_{crit}(t)}\frac{\text{d}u_{crit}(t)}{\text{d}t}\unicode[STIX]{x1D708}_{1}(y).\end{eqnarray}$$

This expression, together with the definition of the Sherwood number given in equation (2.15) provides an equation for the self-similar solution $\unicode[STIX]{x1D708}_{1}$

(3.6) $$\begin{eqnarray}\ell ^{2}(t)\frac{\text{d}}{\text{d}t}\ln u_{crit}(t)=Sh=\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D708}_{1}(y)}{\unicode[STIX]{x2202}y^{2}}\frac{1}{\unicode[STIX]{x1D708}_{1}(y)}.\end{eqnarray}$$

Since the self-similar solution $\unicode[STIX]{x1D708}_{1}$ does not include an explicit time dependency, we make the essential observation that for the critical case the Sherwood number is constant. Thus equation (3.6) provides two results. First, it is an ordinary differential equation for the boundary salinity $u_{crit}$ leading to a constant Sherwood number. Second, it is an equation providing the self-similar solution $\unicode[STIX]{x1D708}_{1}$ for this external control.

3.1.1 The external control required for the critical case

We start with interpreting equation (3.6) in terms of $u_{crit}$ . Observe that it involves two unknowns, $u_{crit}$ and $\ell$ . These unknowns are related by mass conservation of salt in the whole fluid domain, as there are no nucleation events.

In other words, mass conservation for salt gives us in terms of non-scaled variables

(3.7) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\int _{x_{0}(t)}^{x_{0}(t)+\ell (t)}u(t,x)\,\text{d}x=0.\end{eqnarray}$$

Using the definition of $\unicode[STIX]{x1D708}_{1}$ and transforming the above to the domain $[0,1]$ yields

(3.8) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left[\ell (t)u_{crit}(t)\int _{0}^{1}\unicode[STIX]{x1D708}_{1}(y)\,\text{d}y\right]=0.\end{eqnarray}$$

Since the integral is independent of time, we can conclude that the product of domain size and critical concentration is independent of time, i.e.

(3.9) $$\begin{eqnarray}\ell (t)u_{crit}(t)=u_{crit,0}\ell =1.\end{eqnarray}$$

Combining equations (3.6) and (3.9) we obtain

(3.10) $$\begin{eqnarray}u_{crit}(t)=\left(1-\frac{t}{t_{\infty }}\right)^{-1/2},\end{eqnarray}$$

where the integration constant is identified as the time before the brine salinity blows up, which is related to the Sherwood number by

(3.11) $$\begin{eqnarray}t_{\infty }=\frac{1}{2Sh}.\end{eqnarray}$$

We note that from equations (3.9) and (3.11) that the size of the brine region satisfies $\ell (t)=(1-t/t_{\infty })^{1/2}$ . Thus consistent with the notion of mass conservation, as the salinity blows up when $t$ approaches $t_{\infty }$ , the brine domain vanishes, $\ell (t_{\infty })=0$ .

3.1.2 The salinity in the critical case

We now solve equation (3.6) to obtain $\unicode[STIX]{x1D708}_{1}$ and from this we find the salinity. From the definition of $\unicode[STIX]{x1D708}$ in (2.18) the boundary conditions are

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{1}(0)=\unicode[STIX]{x1D708}_{1}(1)=1.\end{eqnarray}$$

The solution to (3.6) satisfying these conditions is

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{1}(y)=\frac{\cosh \left(Sh^{1/2}\left(y-\frac{1}{2}\right)\right)}{\cosh \left(\frac{1}{2}Sh^{1/2}\right)},\end{eqnarray}$$

where $\cosh$ is the hyperbolic cosine function. Equation (3.13) defines a family of solutions depending on the Sherwood number.

The above solutions are obtained in the absence of nucleation, in the sense that the constraint $\unicode[STIX]{x1D708}\geqslant \unicode[STIX]{x1D707}$ has not been enforced. The self-similar solution $\unicode[STIX]{x1D708}_{1}$ has a minimum at $y=1/2$ , thus the nucleation threshold is given as

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{1}(1/2)=\frac{1}{\cosh \left(\frac{1}{2}Sh^{1/2}\right)}>\unicode[STIX]{x1D707}.\end{eqnarray}$$

Since the self-similar solution is not valid across nucleation events, we interpret inequality (3.14) as a constraint on admissible Sherwood numbers for the critical solution, i.e. the freezing must be sufficiently slow such that

(3.15) $$\begin{eqnarray}Sh\leqslant \left(2\text{arcCosh}\frac{1}{\unicode[STIX]{x1D707}}\right)^{2}.\end{eqnarray}$$

This bound on the Sherwood number is essential, as it defines the largest possible constant Sherwood number for this system. When inequality (3.15) is violated over an extended time period, i.e. for a constant Sherwood number, nucleation will happen, and thus the length scale $\ell$ will be reduced, thus also reducing the Sherwood number. Conversely, we can consider the equality (3.15) as a characteristic property of nucleation events. We explore these concepts in § 4.

3.2 A super-critical, fractal case

Figure 3. A schematic plot illustrating the change in the solution $u$ due to spontaneous nucleation in the middle of the domain.

Super-critical freezing appears when cooling dominates diffusion. Due to this, an increase of the critical salinity is encountered and $\unicode[STIX]{x1D708}$ , the ratio between the critical salinity and the minimal one, decreases until it reaches $\unicode[STIX]{x1D707}$ , at which time nucleation occurs. For symmetry reasons we expect that nucleation happens near the middle of a brine domain, we approximate the splitting of the domain into two equal parts – see figure 3. This leads to two domains of halved lengths and, by (2.15), to a decrease in the Sherwood number. A lower Sherwood number together with a steeper salinity gradient near the newly created boundaries lead to a local diffusion-dominated solution profile until freezing starts to dominate again and the process repeats itself. The period where the solution is diffusion dominated is essential in order to transition from a non-symmetric profile after a nucleation event (solid line, figure 3), towards an asymptotically symmetric profile before the next nucleation event happens. Following this scenario, one may expect a fractal structure similar to a Cantor set for the fluid domain.

In the spirit of the critical case described in § 3.1 we seek a self-similar solution in the super-critical case, which will provide insight as to the structures arising during rapid freezing. However, in this case nucleation, and hence brine interval splitting, is encountered at certain times, which is the underlying mechanism for generating fractal solutions. Since $\ell$ is discontinuous in time we understand that, unlike the critical case, we cannot have a constant Sherwood number and thus a freezing-rate scaling (3.10) will not work. Indeed, in this section, we will show that for an external control of the critical salinity which is strictly faster than (3.10), persistent nucleation can occur. The explicit solution derived herein thus provides a connection between the external control on the critical salinity and the relative freezing and nucleation processes.

As the domain splits into smaller sub-domains, the typical length scale becomes smaller and the same should happen with the time scale characterizing the process between two nucleation events. To reflect this, a nonlinear transformation of the time variable is introduced, with respect to which the dimensionless solution $\unicode[STIX]{x1D708}$ has a periodic behaviour

(3.16) $$\begin{eqnarray}t=f(\unicode[STIX]{x1D702}).\end{eqnarray}$$

The function $f$ is a nonlinear, continuous and increasing transformation of time, which will be specified later. At present we require that the transformation ensures that $\unicode[STIX]{x1D708}$ is periodic with period 1 in the transformed time variable $\unicode[STIX]{x1D702}$ .

Figure 4. A conceptual illustration of the changes in the brine domain size during fractal-forming freezing. Regions of brine are represented as grey, with the darker overlay representing the smoothed, $\unicode[STIX]{x1D702}$ -independent approximation of a continuous brine sub-domain size evolving through time.

In the critical case discussed in § 3.1 the Sherwood number was constant in time. Here we use the ansatz that the Sherwood number has the same periodic structure as the domain, i.e. that it is one-periodic with respect to the time variable $\unicode[STIX]{x1D702}$ :

(3.17) $$\begin{eqnarray}Sh(\unicode[STIX]{x1D702})=Sh(\unicode[STIX]{x1D702}+1).\end{eqnarray}$$

In the super-critical case the brine domain $\unicode[STIX]{x1D6FA}$ is complex, consisting of several sub-domains (intervals). Since the Sherwood number is no longer constant in the nucleation case, we introduce a new parameter to characterize the system in this regime. Thus let $\unicode[STIX]{x1D706}>1$ be a similarity parameter, such that between two nucleation events each interval decreases $\unicode[STIX]{x1D706}$ times as the ice boundaries propagate inwards. This non-dimensionless parameter thus gives the balance between freezing and nucleation. In particular, small values of $\unicode[STIX]{x1D706}$ close to $1$ represent nucleation-dominated regimes, while large values of $\unicode[STIX]{x1D706}$ represent freezing-dominated regimes.

The definition of the time variable $\unicode[STIX]{x1D702}$ is such that nucleation events occur at integer values, which we will denote for clarity by $i$ . The definition of the similarity variable $\unicode[STIX]{x1D706}$ thus implies that

(3.18) $$\begin{eqnarray}\ell ((i+1)^{-})=\frac{1}{\unicode[STIX]{x1D706}}\ell (i^{+}).\end{eqnarray}$$

The superscripts $\pm$ should be interpreted as the right and left limits, see figure 4. Further, once a spontaneous nucleation event is encountered in an interval each fluid sub-domain splits into two halves, thus $2\ell (i^{+})=\ell (i^{-})$ and the total decrease of the domain over one self-similarity period is given by:

(3.19) $$\begin{eqnarray}\ell (i^{+})=\frac{1}{2\unicode[STIX]{x1D706}}\ell ((i-1)^{+})=\frac{1}{(2\unicode[STIX]{x1D706})^{i}}.\end{eqnarray}$$

Additionally, we can express mass conservation of salt in the system by accounting for the number of brine regions $n(i^{+})$

(3.20) $$\begin{eqnarray}n(i^{+})\ell (i^{+})\bar{u}(i^{+})=n(0)\ell (0)\bar{u}(0)=\bar{u}(0),\end{eqnarray}$$

where

(3.21) $$\begin{eqnarray}\displaystyle \bar{u}(i^{+})=\frac{1}{|\unicode[STIX]{x1D6FA}(i^{+})|}\int _{\unicode[STIX]{x1D6FA}(i^{+})}u(i^{+},x)\,\text{d}x & & \displaystyle\end{eqnarray}$$

stands for the average salinity in the brine at time $i^{+}$ . The quantity of $\bar{u}(i^{+})$ is unknown but can be approximated by $u_{crit}(i^{+})$ . Combining this with the initial conditions assuming a single brine region of length one allows us to deduce the characteristic domain size

(3.22) $$\begin{eqnarray}\ell (i^{+})=\frac{\bar{u}(0)}{n(i^{+})\bar{u}(i^{+})}\approx \frac{1}{2^{i}u_{crit}(i^{+})}.\end{eqnarray}$$

Here we also used that by definition $n(i^{+})=2^{i}$ . Eliminating $i$ from equations (3.19) and (3.22) gives

(3.23) $$\begin{eqnarray}\frac{\ln \ell (i^{+})}{\ln 2\unicode[STIX]{x1D706}}=\frac{\ln \ell (i^{+})+\ln u_{crit}(i^{+})}{\ln 2}.\end{eqnarray}$$

From equation (3.23) we can now express $\ell$ as an explicit function of $u_{crit}$

(3.24) $$\begin{eqnarray}\ell (i^{+})=u_{crit}(i)^{-1/\unicode[STIX]{x1D6FE}},\end{eqnarray}$$

where the exponent is defined by the similarity parameter as

(3.25) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}=\frac{\ln \unicode[STIX]{x1D706}}{\ln 2\unicode[STIX]{x1D706}}.\end{eqnarray}$$

This exponent represents an essential scaling parameter of the super-critical case and will appear throughout the following results. Note that for large $\unicode[STIX]{x1D706}$ , we obtain $\unicode[STIX]{x1D6FE}{\nearrow}1$ . In the limit case, when $\unicode[STIX]{x1D6FE}=1$ , we recover the results from § 3.1 since by (3.10) the brine region in that case satisfies $\ell (t)=1/u_{crit}(t)$ . On the other hand, when $\unicode[STIX]{x1D706}{\searrow}1$ one gets $\unicode[STIX]{x1D6FE}{\searrow}0$ implying the other extreme dominated by nucleation. Under these conditions the domain freezes all at once from the inside.

Inserting the above into the definition of the Sherwood number, together with the ansatz (3.17) yields

(3.26) $$\begin{eqnarray}u_{crit}(i)^{-(1+2/\unicode[STIX]{x1D6FE})}\frac{\text{d}}{\text{d}t}u_{crit}(t)=Sh(0^{+})\end{eqnarray}$$

at any time $t=f(i)$ when nucleation occurs.

It is natural to consider that the external control $u_{crit}$ evolves continuously. As an example we observe that equation (3.26) is satisfied by the following equation, similar in structure to (3.10)

(3.27) $$\begin{eqnarray}u_{crit}(t)=\left(1-\frac{t}{t_{\infty }}\right)^{-(\unicode[STIX]{x1D6FE}/2)}.\end{eqnarray}$$

As in the critical case, $t_{\infty }$ is a time when the entire domain is completely frozen, which can be calculated analogously as before as

(3.28) $$\begin{eqnarray}t_{\infty }=\frac{\unicode[STIX]{x1D6FE}}{2Sh(0^{+})}.\end{eqnarray}$$

We close the section by noting that we can obtain an explicit expression for nonlinear time transformation by combining equations (3.16), (3.24) and (3.27), which gives the expression

(3.29) $$\begin{eqnarray}t=f(\unicode[STIX]{x1D702})=t_{\infty }[1-(2\unicode[STIX]{x1D706})^{-2\unicode[STIX]{x1D702}}].\end{eqnarray}$$

We summarize the results of this section as follows. Equation (3.29) provides the temporal structure of nucleation events given a freezing-rate control of the form of equation (3.27). This freezing-rate control represents a two-parameter family of freezing rates parameterized by $t_{\infty }$ and $\unicode[STIX]{x1D6FE}<1$ , where the limiting case of $\unicode[STIX]{x1D6FE}=1$ corresponds to the self-similar solution without nucleation obtained in the preceding section. The parameter $t_{\infty }$ simply measures the time scale of the process, while all spatial structure is derived from the exponent $\unicode[STIX]{x1D6FE}$ . The freezing process thus defined has a time-periodic Sherwood number, which at nucleation events satisfies equation (3.28). Similarly, we can extract the reduction of the amount of brine between two nucleation events as $\unicode[STIX]{x1D706}$ , related to $\unicode[STIX]{x1D6FE}$ by equation (3.25). The characteristic length scale of a brine domain becomes related to the critical salinity through equation (3.24).

5.1 Convergent solutions in the critical case

Figure 5. (a) The function $\unicode[STIX]{x1D708}$ at various times (darker colour corresponding to later time); (b) the solution surface plotted against dimensionless time.

We start by considering the self-similar salinity profile for the critical case derived in (3.10). Figure 5 presents the numerical results for the critical case. More precisely, it displays the scaled, dimensionless salinity $\unicode[STIX]{x1D708}$ introduced in (2.18). The initial condition in the numerical simulation is

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D708}(x,0)=0.01\cdot |x-.5|+.995,\quad \text{for }x\in [0,1],\end{eqnarray}$$

and the computation is for $Sh=0.5$ . As seen in figure 5(a), the numerical approximation eventually converges to the analytical solution (3.13).

Figure 5(b) presents the solution in $\unicode[STIX]{x1D702}\times \unicode[STIX]{x1D708}$ axes. The lowest surface line in this projection illustrates the time evolution of $\min (u)/\max (u)=\min (\unicode[STIX]{x1D708})$ , that is, the quantity which will trigger secondary nucleation if it drops below the nucleation multiplier $\unicode[STIX]{x1D707}$ defined in (2.9). Thus the simulations reported here are valid for any $\unicode[STIX]{x1D707}\leqslant 0.993$ .

The convergence of the solution to the self-similar profile predicted by theory indicates that the self-similar solution $\unicode[STIX]{x1D708}$ is robust with respect to perturbations, an observation which is supported by other calculations not shown herein.

5.2 The numerical validation of the fractal-forming behaviour

The numerical experiments in this subsection verify that the condition for the external control on the critical salinity given in equation (3.27) indeed leads to binary, Cantor-set-type, fractals. This numerical approach is particularly important since in the super-critical fractal case it is not feasible to obtain an explicit self-similar solution as for the critical case. For this example we use $Sh(0+)=0.02$ , which is consistent with a nucleation threshold of $0.9925$ according to inequality (3.15). The fractal formation is then characterized by the balance between freezing and nucleation, which we set to $\unicode[STIX]{x1D706}=1.03$ for the sake of illustration.

Figure 6. Demonstration of the fractal-forming behaviour for the super-critical conditions given in equation (3.27). Panel (a) shows minimum concentration (lower bound) and concentration profiles (colours), while (b) shows the concentration (colours) plotted on the domains as they split, thus white represents ice. For both panels, time is given with respect to the transformed time $\unicode[STIX]{x1D702}$ .

The salinity profile immediately after nucleation events converges towards an asymmetric self-similar solution after a few nucleation events and thus the initial condition is immaterial for the asymptotic behaviour of the system. To reflect this, we only show the solution after it has converged to the self-similar profile, typically within two to three splitting steps. The self-similar profile obtained in this example is also used as the initial condition of salinity when we return to exploring the non-asymptotic regimes in § 5.3.

The solution is shown in figure 6 for a time span covering several nucleation events. As the numerical simulation evidences, each sub-domain splits into two equal sub-domains with a period 1 in the nonlinear time $\unicode[STIX]{x1D702}$ . This provides a post hoc justification of the assumptions stated in § 3.2.

The self similarity of the solution implies that the process will continue as $\unicode[STIX]{x1D702}\rightarrow \infty$ , leading to a countable infinity of splittings within a finite time (recall that $t_{\infty }$ is finite). This results in a bipartite tree of sub-domains that, at each time, is an approximation of a Cantor-like set. We have verified this behaviour for freezing rates of $\unicode[STIX]{x1D706}=1.01$ to $\unicode[STIX]{x1D706}=1.50$ , the results of which are summarized in figure 7.

Figure 7. This log–log plot summarizes the relationship between remaining fluid fraction and characteristic fluid domain size. The linear trends (in log–log space) confirm the fractal structure of the solution over the full range of freezing rates.

5.3 Applying the approximate method to determine the properties of sea ice

In this subsection we validate the heuristically derived Sherwood number $Sh$ and mean brine-sub-domain length scale $\ell$ derived in § 4. We use this opportunity to consider the pre-asymptotic regime alluded to in the previous section.

5.3.1 Fractal solutions including the pre-asymptotic regime

Figure 8. Comparison of simulated and estimated Sherwood number and average sub-domain size for a solution including the pre-asymptotic regime.

We consider the asymmetric initial condition (from § 5.2) leading to initial asymmetrical splitting and thus the different brine sub-domains have different widths. In this case, the fractal does not have a pure binary structure anymore. This can be seen in figure 8, where each individual branch of the solution approaches an asymptotic regime after about four splittings. Note that for these figures we choose to use dimensionless time $t$ rather than $\unicode[STIX]{x1D702}$ and, to better visualize the results, we use $t-t_{\infty }$ as the $x$ -axis with logarithmic scale.

As the numerical solution is driven by a prescribed $u_{crit}(t)$ , as detailed in § 5.2, we can apply the estimates for the mean length scale as well as the Sherwood numbers, as given by equations (3.22), (4.2) and (4.4). The estimates are compared against the true values obtained from the simulation in the lower part of figure 8. Since the estimate from § 4 does not know the initial salinity distribution, it cannot capture the details at the earliest stage of the simulation. However, for later times we find a reasonable agreement between the estimated and simulated values.

5.3.2 A generic rapid freezing scenario

Figure 9. Comparison of simulated and estimated Sherwood number and average sub-domain size for logistic function boundary conditions.

As a final comparison, we consider a case of rapid freezing where the boundary salinity controlling the freezing process does not follow the monomial scaling given in equation (3.27). In particular, we consider the logistic function

(5.2) $$\begin{eqnarray}u_{crit}(t)=1+\frac{5}{(1+\exp (10(-t+1)))},\quad \text{for }t\in [0,1.5].\end{eqnarray}$$

This choice is motivated by the transition from a slow, diffusion-dominated regime via a nucleation-dominated regime and finally a diffusion-dominated regime. Paired with the non-symmetric initial conditions from § 5.2, we consider this a challenging test for the heuristic estimates derived in § 4.

The numerical solution and the comparison of the characteristic parameters are depicted on figure 9. We first note the rapid establishment of an asymptotic-like fractal structure, in the sense of the division of the original brine domain into eight nearly equal-sized partitions. This provides additional support to the applicability of the analysis of § 3 for general freezing regimes.

Second, we note from the lower figures that the estimates for characteristic brine-domain length and Sherwood number are in general very close to the computed values. Again the discrepancy is related to the salinity distribution at the onset of the freezing process, which is not captured, leading to a prediction of premature nucleation. As the nucleation process is established, the estimates approach the calculated values both of Sherwood number and of domain-size, with the correct prediction of eight fluid domains. As expected, the heuristic algorithm gives a very close match for the final smooth regime (starting from $t\approx 0.7$ ) where the solution is close to the steady state.