Finite-Element Approximation of One-Sided Stefan Problems with Anisotropic, Approximately Crystalline, Gibbs--Thomson Law

We present a finite-element approximation for the one-sided Stefan problem and the one-sided Mullins–Sekerka problem, respectively. The problems feature a fully anisotropic Gibbs–Thomson law, as well as kinetic undercooling. Our approximation, which couples a parametric approximation of the moving boundary with a finite element approximation of the bulk quantities, can be shown to satisfy a stability bound, and it enjoys very good mesh properties which means that no mesh smoothing is necessary in practice. In our numerical computations we concentrate on the simulation of snow crystal growth. On choosing realistic physical parameters, we are able to produce several distinctive types of snow crystal morphologies. In particular, facet breaking in approximately crystalline evolutions can be observed.

💡 Research Summary

The paper introduces a novel finite‑element/parametric scheme for solving one‑sided Stefan and one‑sided Mullins‑Sekerka problems in which the interfacial condition is governed by a fully anisotropic Gibbs–Thomson law together with kinetic undercooling. In classical solid–liquid phase‑change models the Gibbs–Thomson relation is either isotropic or only weakly anisotropic, which limits the ability to reproduce the highly faceted morphologies observed in natural snow crystals. The authors therefore adopt a “approximately crystalline” description of the surface energy: the anisotropy is allowed to be strong enough to generate flat facets separated by sharp edges, yet the energy function remains smooth enough for variational treatment.

The computational framework couples a parametric representation of the moving interface with a bulk finite‑element discretisation of the temperature (or concentration) field. The interface is approximated by a piecewise linear (in 2‑D) or piecewise planar (in 3‑D) mesh whose nodes move with the normal velocity prescribed by the Stefan/Mullins‑Sekerka law. Because the mesh points are advected by the same velocity that defines the interface motion, the mesh automatically adapts to the evolving geometry; no post‑processing mesh‑smoothing or remeshing is required. The bulk equations are solved on a fixed background mesh, while the coupling terms are evaluated on the moving interface mesh, leading to a symmetric, positive‑definite linear system at each time step.

A key theoretical contribution is the derivation of a discrete energy stability estimate. By constructing a discrete analogue of the total free‑energy (bulk diffusion energy plus interfacial energy) and showing that the numerical scheme dissipates this quantity monotonically, the authors prove unconditional stability with respect to the time step size. This result holds for arbitrary anisotropy strength and for any admissible kinetic coefficient, thereby guaranteeing robustness even in regimes where the interface develops sharp corners or facets.

Implementation details include: (i) the computation of the anisotropic curvature term via the surface‑energy tensor evaluated on the local normal; (ii) a consistent treatment of kinetic undercooling, which adds a linear term proportional to the normal velocity in the Gibbs–Thomson condition; and (iii) an efficient assembly of the coupled system using standard finite‑element libraries. The authors also discuss how to choose realistic physical parameters (diffusivity, surface tension, kinetic coefficient, undercooling) for water/ice, enabling quantitative comparison with laboratory snow‑crystal growth.

Numerical experiments focus on snow‑crystal morphogenesis. Starting from a small seed, the algorithm reproduces a wide variety of crystal types—plates, dendrites, columns, and the so‑called “stalk‑prism” shapes—by varying the supersaturation and temperature. Importantly, the simulations capture facet‑breaking events: as the crystal grows, initially flat basal or prism facets develop instabilities that lead to the emergence of sidebranches, a phenomenon that is difficult to obtain with isotropic or weakly anisotropic models. The authors present both 2‑D cross‑sections and full 3‑D visualisations, demonstrating that the mesh quality remains high throughout the evolution, with no need for ad‑hoc smoothing.

In summary, the paper makes four major contributions: (1) a unified finite‑element/parametric discretisation that simultaneously resolves bulk diffusion and interface motion; (2) a rigorous discrete energy stability proof that holds for fully anisotropic Gibbs–Thomson laws and kinetic undercooling; (3) an algorithmic design that preserves mesh quality automatically, eliminating the need for mesh‑regularisation; and (4) a set of realistic snow‑crystal simulations that reveal complex morphological transitions, including facet breaking. The work opens the door to high‑fidelity simulations of anisotropic solidification in other materials, to coupling with additional physics (e.g., solute transport, elastic stresses), and to quantitative validation against experimental data.