2. The ‘extended’ Stefan problem

Moving phase boundary problems such as the one of borehole closure are typically solved mathematically as a Stefan problem. The classical Stefan problem describes two-phase thermal diffusion and movement of their interlaying phase boundary. The rate at which the phase boundary moves is dependent on thermal diffusion in the surrounding media, either supplying heat for melting or removing heat for freezing. As an extension to the classical Stefan problem, Worster (Reference Worster, Batchelor, Moffat and Worster2000) added a molecular diffusion component for the case of a binary solution (e.g. antifreeze–water mixture in a borehole). Now, the motion of the phase boundary is dictated not only by the temperature gradients on either side, but also by the concentration gradient which more fully addresses the case of slush formation in a glacier borehole filled with antifreeze solution. Much less work has been done for this coupled problem (Huppert and Worster, Reference Huppert and Worster1985; Worster, Reference Worster, Batchelor, Moffat and Worster2000) with most applications to sea-ice formation (Feltham and others, Reference Feltham, Untersteiner, Wettlaufer and Worster2006). For the specific phase-boundary problem of borehole melting/freezing, the Stefan problem is solved in cylindrical coordinates (Humphrey and Echelmeyer, Reference Humphrey and Echelmeyer1990), representing an infinite cylinder with the phase boundary moving only in the radial direction. Zagorodnov and others (Reference Zagorodnov1994) extended the classical borehole thermodynamics problem to include a binary solution, but only treated specific cases for borehole dissolution. Below, we describe the physics for melting/freezing of an infinite cylinder with a binary solution and test several simplified cases.

The model for phase-boundary motion with a binary solution is defined by three coupled partial differential equations. The first two are for thermal diffusion

(1)$$\rho _{\rm s}c_{\rm s}\displaystyle{{\partial T_{\rm s}} \over {\partial t}} = \nabla \cdot ( {k_{\rm s}\nabla T_{\rm s}} ) , \;$$

(2)$$\rho _{\rm l}c_{\rm l}\displaystyle{{\partial T_{\rm l}} \over {\partial t}} = \nabla \cdot ( {k_{\rm l}\nabla T_{\rm l}} ) , \;$$

with Eqn (1) for the solid (x _s) and Eqn (2) for the liquid (x _l). Here, T is temperature, t is time, ρ is density, c is specific heat capacity and k is thermal conductivity. The thermal diffusivity is α = k/ρc. The third governing equation is for molecular diffusion in the solution which is defined by the same fundamental equation as above (i.e. Fickian diffusion)

(3)$$\displaystyle{{\partial C} \over {\partial t}} = \nabla \cdot ( {D\nabla C} ) , \;$$

where C is the concentration of solute (in kg m⁻³) and D is the molecular diffusivity.

The three governing equations above are coupled at the phase boundary. The coupling is described mathematically by the Stefan condition, which dictates wall motion based on the heat flux in both solid and liquid media

(4)$$\rho _{\rm s}L\displaystyle{{{\rm d}R} \over {{\rm d}t}} = k_{\rm s}\displaystyle{{\partial T_{\rm s}} \over {\partial r}}-k_{\rm l}\displaystyle{{\partial T_{\rm l}} \over {\partial r}}, \;$$

where L is the latent heat of fusion, R is the phase boundary location (i.e. the borehole radius for cylindrical coordinates) and r is the radial coordinate. Heat flows down the temperature gradient, so a negative gradient in the solid (temperature decreasing away from the hole) causes freezing. The act of freezing gives off thermal energy to the ice at the hole wall which will eventually diffuse out into the bulk ice. During active drilling, there is an imposed heat flux (Q) at the phase boundary, so an additional term can be included in Eqn (4). A secondary Stefan condition is added for the molecular diffusion problem

(5)$$\displaystyle{{\partial C} \over {\partial r}} = {-}\displaystyle{C \over D}\displaystyle{{{\rm d} R} \over {{\rm d} t}}.$$

This condition states that the solute flux at the hole wall is dependent on the rate at which the wall moves. Conceptually, Eqn (5) expresses either solute rejection and accumulation in the case of borehole freezing (∂R/∂t is negative) or solute dilution in the case of borehole melting (∂R/∂t is positive). Thus, even when no additional mass of solute is added to the solution, freezing of the hole wall acts as a source term for the solute concentration and melting of the hole wall as a sink term.

To close the model, one additional equation is necessary which describes the concentration-dependent freezing temperature and fixes the thermal boundary condition at the hole wall for both Eqns (1) and (2),

(6)$$T_{\rm f} = T_0 + f( C ) .$$

Here, T ₀ is the pure liquid freezing temperature and f(C) is a function for freezing-point depression. We assign the freezing-point depression empirically (Supplementary Material S.1).

All model results below are reported in non-dimensional units: r* = r/R ₀, R* = R/R ₀ where R ₀ is the drilled radius (0.04 m in our case), T* = T/T _∞ where T _∞ is the bulk ice temperature (−20 °C in our case), t* = t/t ₀ where $t_0 = R_0^2 L/k_{\rm s}\vert {T_\infty } \vert$. The model domain is 1-dimensional (radial), and the imposed heat flux for melting out the initial hole is constant (Q = 2.5  kW). Adjusting the problem-specific parameters (i.e. R ₀, T _∞, and Q) would change the precise timing and magnitude of borehole evolution, but the qualitative conclusions presented below are consistent for any choice of parameters. For cases in which the numerical model can be tested analytically, those tests were done (Supplementary Material S.2).

3 Sequence of borehole stabilization

We use the equations described above to simulate borehole melting/freezing in the presence of an antifreeze solution for an ensemble of injection timings and volumes. Unlike prior studies, we emphasize the hot-point case where downward drilling motion is relatively slow but antifreeze can be injected soon after the drill passes through the borehole. For every simulation, the antifreeze solution is aqueous methanol and is considered ‘well mixed’ immediately after injection. Hence, the solution concentration and temperature are constant throughout the borehole while the solution dissolves the hole wall, reducing the model to Eqns (1) and (4). During this instantaneous-mixing phase, our model matches Humphrey and Echelmeyer (Reference Humphrey and Echelmeyer1990) exactly, including the numerics of the logarithmic transform and non-dimensionalization (Jarvis and Clarke, Reference Jarvis and Clarke1974; Humphrey and Echelmeyer, Reference Humphrey and Echelmeyer1990). The concentration within the hole is updated purely based on volume changes of the hole, and the solution temperature is assumed to match the depressed freezing temperature. The mixing process is exothermic, so thermal energy is added to the solution immediately after antifreeze injection.

We divide the simulation into four stages based on the dominating process for each stage. An example hole radius through time and a representative temperature profile from each stage are shown in Figure 1:

Fig. 1. A time series of borehole melting and refreezing. The labeled black line shows the hole radius through time from when drilling starts (bottom) to where the borehole stabilizes with slush forming in the solution (top). The horizontal scale changes at stage 4 to emphasize peak dissolution and subsequent slush formation (indicated by the expanding dashed lines). Representative temperature profiles are shown for each stage (1–4) with black the temperature in the ice, gray the temperature in a previous stage, blue the temperature in the solution which changes with the solution concentration, and the red area showing the warm thermal annulus in the ice. The inset for stage 4 is a cartoon representation which shows the competition between thermal and molecular diffusion that causes constitutional supercooling (slush) where the solution is colder than the liquidus line.

Stage 1. Melting – The melting stage is defined by a heat flux that follows the borehole wall as it moves outward. In addition to the energy for melting ice and moving the borehole wall outward, there is inevitably some amount of energy that goes toward warming ice in an annular zone outside the hole wall. The amount of warming in this thermal annulus is dependent on the rate of drilling with the fastest drilling being the most efficient to direct its energy toward melting while only minimally warming the ice (in the Supplementary Material S.3, we assess the influence of drilling speed on the two-dimensionality of an advection-diffusion problem for the down going drill).

Stage 2. Freezing – After the drill melts out the hole and the heat flux at the hole wall is removed, the energy in the hole begins to diffuse outward into the bulk ice; that is, the hole begins to freeze shut. During this process, the warm thermal annulus continues to grow as the latent energy within the hole is converted to thermal energy in the encroaching ice. The rate of freezing is rapid at first when the temperature gradient at the hole wall is greatest. Then, freezing slows down as the temperature gradient flattens out but the hole is still a substantial fraction of its original volume. Finally, when the hole approaches closure, the rate of freezing accelerates again because the total hole volume to freeze decreases non-linearly with the decreasing radius (i.e. volume is dependent on the square of the radius).

Stage 3. Dissolution – Antifreeze can be injected into the borehole to stop Stage 2 before the hole freezes shut. In all simulations, the injected antifreeze is at 0 °C, but directly after injection, the solution mixes and cools to the freezing temperature of the newly mixed solution. For some time, the solution is colder than the surrounding ice within the thermal annulus. Therefore, energy from the annulus initially moves back toward the hole center and causes the hole wall to re-melt (or dissolve). The total amount of dissolution depends on the size of the thermal annulus, as well as the total amount and rate of antifreeze injection. We assume that the dissolution process happens quickly after injection (all simulations shown here are on the scale of ~hours) and the hole therefore remains well mixed. If, on the other hand, the effect of molecular diffusion is considered important, the antifreeze solution would be diluted near the hole wall as ice dissolves, and the rate-limiting factor for dissolution would be the molecular diffusion of antifreeze toward the hole wall which is orders of magnitude slower than thermal diffusion for most solutions.

Stage 4. Refreezing (Slush) – Eventually, the dissolution in Stage 3 dissipates all the energy from the thermal annulus. If, at this stage, the solution temperature is colder than the bulk-ice temperature, the hole wall continues to dissolve and refreezing as slush is avoided entirely. If, on the other hand, the thermal annulus was large enough to dissolve the hole and warm the solution above the bulk-ice temperature, the energy that was originally directed toward dissolution moves back into the bulk ice and the hole refreezes. In the case of a well-mixed hole, refreezing happens as in Stage 2, by accretion on the hole wall. However, the timescale for this process (~days) is significantly longer than Stage 3 and unless there is some forced mixing within the hole, a solute gradient will be established. In this case, solute accumulates within a short distance (~millimeters) of the hole wall, with the exact characteristic length scale determined by the molecular diffusivity ($\sqrt {Dt}$). As solute accumulates, the solution temperature is depressed to the freezing temperature at the hole wall. Then, thermal and molecular diffusion compete within the hole. Thermal diffusion acts to distribute the latent energy of freezing throughout the solution, while molecular diffusion acts to distribute the accumulated solute that is rejected from the freezing ice lattice. For all cases in which the non-dimensional Lewis number (i.e. thermal over molecular diffusivity, Le = α/D) is greater than 1, the solution cools more quickly than solute moves. Then, slush forms in an area of ‘constitutional supercooling’ where the solution temperature is effectively below the freezing temperature (Fig. 1 stage 4 inset) (Worster, Reference Worster, Batchelor, Moffat and Worster2000). Considering aqueous methanol as the antifreeze agent Le ≫ 1, so we assume that all refreezing is in the form of slush.

To calculate the amount of dissolved and refrozen ice in our simulation results, we first assume that the borehole radius eventually converges to equilibrium, R _eq. This radius is dependent on the total mass of injected antifreeze, M, and the bulk-ice temperature, T _∞,

(7)$$T_{\rm f}\left({\displaystyle{M \over {\pi R_{{\rm eq}}^2 }}} \right) = T_\infty , \;$$

where $C_{{\rm eq}} = M/\pi R_{{\rm eq}}^2$ is the solution concentration at equilibrium volume, $V_{{\rm eq}} = \pi R_{{\rm eq}}^2$ (per unit depth). Hence, the borehole changes size until the concentration is such that the depressed freezing temperature equals the bulk-ice temperature. The fraction of dissolved ice (Fig. 2a) is the ratio between the maximum dissolved borehole volume (V _④; Figure 1 stage 4) and the volume at injection (V _②; Figure 1 stage 2). It gives a measure of the over-dissolution of the hole before settling to the equilibrium volume. The fraction of refrozen ice (Fig. 2b) is the ratio between the equilibrium volume (V _eq) and the maximum dissolved borehole volume (V _④). It gives an indication of the amount of slush formed. In some cases, enough antifreeze is injected that V _④ $ = V_{{\rm eq}}$ and the refreezing stage is avoided altogether (white space in Fig. 2b).

Fig. 2. Fraction dissolved (a) and fraction refrozen (b) results for an ensemble of simulations with variable injection timing and total mass of injected antifreeze. Each of the dots represents the result from one simulation, with its location on the plot being the timing of antifreeze injection and the equilibrium radius (proportional to the total mass of injected antifreeze through Eqn (7)). Fractional volumes shown by the colormaps are calculated after the entire simulation runs as described in the text, and colors are interpolated between individual simulations (i.e. individual dots). The radius at injection (dashed line) provides a reference for the state of the hole before any dissolution or refreezing. The pure antifreeze limit (solid line) is the equilibrium radius when the solution is pure antifreeze at the time of injection.