<figure><img src="image_1.png" /></figure>

HAL
open science

# A phase field model for snow crystal growth in three dimensions

Gilles Demange, Helena Zapolsky, Renaud Patte, Marc Brunel

► To cite this version:

Gilles Demange, Helena Zapolsky, Renaud Patte, Marc Brunel. A phase field model for snow crystal growth in three dimensions. npj Computational Materials, 2017, 3 (1), 10.1038/s41524-017-0015-1. hal-01861344

HAL Id: hal-01861344

https://normandie-univ.hal.science/hal-01861344

Submitted on 29 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE OPEN

# A phase field model for snow crystal growth in three dimensions

Gilles Demange<fnref n="1" />, Helena Zapolsky<fnref n="1" />, Renaud Patte<fnref n="1" /> and Marc Brunel²

Snowflake growth provides a fascinating example of spontaneous pattern formation in nature. Attempts to understand this phenomenon have led to important insights in non-equilibrium dynamics observed in various active scientific fields, ranging from pattern formation in physical and chemical systems, to self-assembly problems in biology. Yet, very few models currently succeed in reproducing the diversity of snowflake forms in three dimensions, and the link between model parameters and thermodynamic quantities is not established. Here, we report a modified phase field model that describes the subtlety of the ice vapour phase transition, through anisotropic water molecules attachment and condensation, surface diffusion, and strong anisotropic surface tension, that guarantee the anisotropy, faceting and dendritic growth of snowflakes. We demonstrate that this model reproduces the growth dynamics of the most challenging morphologies of snowflakes from the Nakaya diagram. We find that the growth dynamics of snow crystals matches the selection theory, consistently with previous experimental observations.

npj Computational Materials (2017)3:15; doi:10.1038/s41524-017-0015-1

## INTRODUCTION

Snowflake growth in supersaturated atmosphere is one of the most familiar, and at the same time scientifically challenging physical phenomena.<fnref n="1" /> Beyond aesthetic fascination for their symmetric shape, snow crystals provide a unique example of self-patterning systems.² Early experiments on artificial snow crystals in cold chamber led by Nakaya,³ revealed this fundamental phase transition resulted in a wide manifold of patterns, exclusively determined by supersaturation and temperature. This empirical dependency was later formalized in the meteorological classification of Magono,⁴ and the snow crystals morphology diagram of Nakaya.⁵ Despite this clear experimental picture, the underlying physical rules remain thoroughly debatable,⁶ as evidenced by the numerous models proposed to understand the variety of snowflake shapes, such as the surface diffusion model of Mason et al.,⁷ the quasi-liquid layer approach of Lacmann et al.,⁸ and the layer nucleation rates theory by Nelson.⁹

Alternatively, many simulation methods were developed to reproduce the growth dynamics of snowflakes. First step toward the comprehension of ice crystal growth was provided by molecular dynamics simulations.¹⁰⁻¹² Unfortunately, such simulations are still confined to space and time scales by several orders inferior to snowflake characteristic scales.⁴ Significant results were also achieved using the mesoscopic approach, such as the cellular automata model of Gravner and Griffeath.¹³ Though this model remarkably describes the morphology of snowflakes, the numerous parameters involved can hardly be related to physical quantities. This shortcoming was partly corrected in ref. ¹⁴. However, it is not perfectly clear how additional physical mechanisms, such as fluid dynamics can be added to the procedure. This might indeed be required to simulate the impact of air flows on the morphology of snowflakes. From this perspective, a 3D sharp interface model was developed by Barrett et al.¹⁵. Different snowflake morphologies were simulated.

Nevertheless, the side branching,⁶, ¹⁶ surface markings,¹⁷⁻¹⁹ and coalescence²⁰ of ice crystals could not be reproduced in this framework. Besides, only small supersaturations were prospected. Consequently, only the bottom of the Nakaya diagram was explored. This limitation comes from the Laplacian approximation framework, and the numerical cost of interface parametrization.²¹ The phase field model has the decisive advantage to overcome explicit tracking of the sharp interface, by spreading it out over a small layer.²¹ As a matter of fact, it has become standard to simulate the dendritic growth in alloy solidification.²²⁻²⁴ However, except a first mathematical study by Barret et al. in ref. ²⁵, the phase field approach was almost never used to simulate snow crystal growth in three dimensions. Indeed, until recently, it was assumed that the phase field approach was unable to reproduce facet formation and destabilization.²⁶ Yet, Debierre and Karma suggested in ref. ²⁷, that phase field could mimic faceting using highly anisotropic surface tension.

In this paper, we report the simulations of snow crystals growth in three dimensions, using a modified phase field model. More generally, by addressing the arduous study-case of ice crystal growth, this work fulfils the need of material science for faceted dendrite modelling at mesoscopic and macroscopic scale. A new surface tension anisotropy function accounting for the 6-fold horizontal and 2-fold vertical symmetry of snowflakes was derived. Contrary to previous phase field models,²²⁻²⁴ a supplementary anisotropy function, and anisotropic diffusion terms were also included to simulate the vertical anisotropy of snowflakes.⁶ To mimic faceting in snowflakes, the 2D regularisation algorithm of Eggleston et al.²⁸ was extended to three dimensions. As a result, the model reproduces the growth of the main snowflake morphologies of the Nakaya diagram, varying only four phenomenological parameters. It is shown that these parameters can be related to physical quantities. Simulated snowflakes show excellent agreement with experimental observations.¹⁶⁻¹⁹, ²⁹ Their

<footnote n="1">¹GPM, UMR CNRS 6643, University of Rouen, Saint Étienne du Rouvray 76575, France and ²CORIA UMR 6614, University of Rouen, Saint Étienne du Rouvray 76575, France
Correspondence: Gilles Demange (gilles.demange@univ-rouen.fr)</footnote>
ARTICLE OPEN

# A phase field model for snow crystal growth in three dimensions

Gilles Demange<fnref n="1" />, Helena Zapolsky<fnref n="1" />, Renaud Patte<fnref n="1" /> and Marc Brunel²

Snowflake growth provides a fascinating example of spontaneous pattern formation in nature. Attempts to understand this phenomenon have led to important insights in non-equilibrium dynamics observed in various active scientific fields, ranging from pattern formation in physical and chemical systems, to self-assembly problems in biology. Yet, very few models currently succeed in reproducing the diversity of snowflake forms in three dimensions, and the link between model parameters and thermodynamic quantities is not established. Here, we report a modified phase field model that describes the subtlety of the ice vapour phase transition, through anisotropic water molecules attachment and condensation, surface diffusion, and strong anisotropic surface tension, that guarantee the anisotropy, faceting and dendritic growth of snowflakes. We demonstrate that this model reproduces the growth dynamics of the most challenging morphologies of snowflakes from the Nakaya diagram. We find that the growth dynamics of snow crystals matches the selection theory, consistently with previous experimental observations.

npj Computational Materials (2017)3:15; doi:10.1038/s41524-017-0015-1

## INTRODUCTION

Snowflake growth in supersaturated atmosphere is one of the most familiar, and at the same time scientifically challenging physical phenomena.<fnref n="1" /> Beyond aesthetic fascination for their symmetric shape, snow crystals provide a unique example of self-patterning systems.² Early experiments on artificial snow crystals in cold chamber led by Nakaya,³ revealed this fundamental phase transition resulted in a wide manifold of patterns, exclusively determined by supersaturation and temperature. This empirical dependency was later formalized in the meteorological classification of Magono,⁴ and the snow crystals morphology diagram of Nakaya.⁵ Despite this clear experimental picture, the underlying physical rules remain thoroughly debatable,⁶ as evidenced by the numerous models proposed to understand the variety of snowflake shapes, such as the surface diffusion model of Mason et al.,⁷ the quasi-liquid layer approach of Lacmann et al.,⁸ and the layer nucleation rates theory by Nelson.⁹

Alternatively, many simulation methods were developed to reproduce the growth dynamics of snowflakes. First step toward the comprehension of ice crystal growth was provided by molecular dynamics simulations.¹⁰⁻¹² Unfortunately, such simulations are still confined to space and time scales by several orders inferior to snowflake characteristic scales.⁴ Significant results were also achieved using the mesoscopic approach, such as the cellular automata model of Gravner and Griffeath.¹³ Though this model remarkably describes the morphology of snowflakes, the numerous parameters involved can hardly be related to physical quantities. This shortcoming was partly corrected in ref. ¹⁴. However, it is not perfectly clear how additional physical mechanisms, such as fluid dynamics can be added to the procedure. This might indeed be required to simulate the impact of air flows on the morphology of snowflakes. From this perspective, a 3D sharp interface model was developed by Barrett et al.¹⁵. Different snowflake morphologies were simulated.

Nevertheless, the side branching,⁶, ¹⁶ surface markings,¹⁷⁻¹⁹ and coalescence²⁰ of ice crystals could not be reproduced in this framework. Besides, only small supersaturations were prospected. Consequently, only the bottom of the Nakaya diagram was explored. This limitation comes from the Laplacian approximation framework, and the numerical cost of interface parametrization.²¹ The phase field model has the decisive advantage to overcome explicit tracking of the sharp interface, by spreading it out over a small layer.²¹ As a matter of fact, it has become standard to simulate the dendritic growth in alloy solidification.²²⁻²⁴ However, except a first mathematical study by Barret et al. in ref. ²⁵, the phase field approach was almost never used to simulate snow crystal growth in three dimensions. Indeed, until recently, it was assumed that the phase field approach was unable to reproduce facet formation and destabilization.²⁶ Yet, Debierre and Karma suggested in ref. ²⁷, that phase field could mimic faceting using highly anisotropic surface tension.

In this paper, we report the simulations of snow crystals growth in three dimensions, using a modified phase field model. More generally, by addressing the arduous study-case of ice crystal growth, this work fulfils the need of material science for faceted dendrite modelling at mesoscopic and macroscopic scale. A new surface tension anisotropy function accounting for the 6-fold horizontal and 2-fold vertical symmetry of snowflakes was derived. Contrary to previous phase field models,²²⁻²⁴ a supplementary anisotropy function, and anisotropic diffusion terms were also included to simulate the vertical anisotropy of snowflakes.⁶ To mimic faceting in snowflakes, the 2D regularisation algorithm of Eggleston et al.²⁸ was extended to three dimensions. As a result, the model reproduces the growth of the main snowflake morphologies of the Nakaya diagram, varying only four phenomenological parameters. It is shown that these parameters can be related to physical quantities. Simulated snowflakes show excellent agreement with experimental observations.¹⁶⁻¹⁹, ²⁹ Their

<footnote n="1">¹GPM, UMR CNRS 6643, University of Rouen, Saint Étienne du Rouvray 76575, France and ²CORIA UMR 6614, University of Rouen, Saint Étienne du Rouvray 76575, France
Correspondence: Gilles Demange (gilles.demange@univ-rouen.fr)</footnote>
growth satisfies the microscopic solvability theory,³⁰⁻³³ consistently with experiments.³⁴

## MODEL

Snowflakes growth in supersaturated water vapour was simulated using a phase field approach in three dimensions, based on a methodology developed in.²² In this model, two coupled variables $\phi$ and $u$ are considered. $\phi$ is an order parameter referring to the ice (+1) and vapour (-1) phases. The ice/vapour interface is described by a continuous variation of $\phi$, connecting -1 and +1. $u = (c - c'_{sat})/c'_{sat}$ is the reduced supersaturation of water vapour, where $c'_{sat}(T)$ is the saturation number density of vapour above ice, at temperature $T$. At initial time, the reduced supersaturation $u_0$ is homogeneous. The growth kinetics of snowflakes is governed by two non conservative phase field equations. Their adimensionalized form is given by:

$$A(n)^2 \partial_t \phi = -f'(\phi) + \lambda B(n) g'(\phi) u + \frac{1}{2} \nabla_r \cdot \left( |\nabla \phi|^2 \frac{\partial [A(n)^2]}{\partial \nabla \phi} + A(n)^2 \nabla_r \phi \right) \quad (1)$$

$$\partial_t u = \tilde{D} \nabla_r \cdot (q(\phi) \nabla_r u) - \frac{L_{sat}}{2} B(n) \partial_t \phi, \quad (2)$$

where space and time are scaled by the interface width $W_0$, and the relaxation time $\tau_0$, respectively. In the Eq. (1), the double well potential $f(\phi) = -\phi^2/2 + \phi^4/4$ is the free energy density $f_{iv}(c, T', \phi)$ of the ice/vapour system, at temperature $T$, and saturation concentration $c = c'_{sat}$: $f(\phi) = f_{iv}(c = c'_{sat}, T' = T, \phi)$. The second term in the Eq. (1) accounts for the coupling between $u$ and $\phi$, and promotes ice phase growth in supersaturated atmosphere, where $c > c'_{sat}$. Formally, it corresponds to the first order term in the Taylor expansion of the bulk potential, for $c$ in the neighbourhood of $c'_{sat}$. The coupling constant $\lambda$ can thus be computed by $\lambda = c_{sat}/(30H) \partial_c f_{iv}|_{c_{sat}}$, where $H$ is the free energy barrier. $g'(\phi) = (1 - \phi^2)^2$ is an interpolation function introduced in.²² Its form allows to keep the bulk potential minima at $\phi = \pm 1$ for any $u$. To describe the strong anisotropy of snowflakes growth along the vertical axis,⁶ we propose to introduce the kinetic anisotropy function $B(n) = \sqrt{n_x^2 + n_y^2 + \Gamma^2 n_z^2}$, where $n = -\nabla \phi / |\nabla \phi|$ is the unit normal vector of $\phi$. This anisotropic contribution may be justified by the premelted layer growth theory, stating that the primary habit results from different attachment rates of water molecules between prismatic and basal facets.⁸ The parameter $\Gamma > 0$ governs the preference between horizontal and vertical growth, called the primary habit of snowflakes. It can be empirically related to the temperature.³⁵ The addition of the first two terms in the Eq. (1) corresponds to the anisotropic thermodynamic driving force.

The last two terms in the Eq. (1) describe the ice/vapour interface formation and propagation,²² where $A(n)$ is the surface tension anisotropy function. In this work, a new expression of the anisotropy function was derived: $A(n) = 1 + \varepsilon_{xy} \cos(6\theta) + \varepsilon_z \cos(2\psi)$, where $\theta = \arctan(n_y/n_x)$ and $\psi = \arctan(\sqrt{n_x^2 + n_y^2}/n_z)$ are the polar and azimuthal angles, respectively. It accounts for both the horizontal 6-fold symmetry, and the vertical planar symmetry of snowflakes. $\varepsilon_{xy}$ and $\varepsilon_z$ are anisotropy constants.

The Eq. (2) describes one-sided diffusion of water molecules in vapour, and vapour condensation on ice. Diffusion is controlled by two quantities: the function $q(\phi) = 1 - \phi$, which prohibits diffusion within ice, and the reduced diffusion coefficient $\tilde{D} = D\tau_0/W_0^2$, where $D$ is the diffusion coefficient. The second term in the Eq. (2) accounts for conversion of vapour into ice. $L_{sat}$ can thus be interpreted as the rate of water depletion in vapour, via molecule attachment at the interface. In this study, $L_{sat}$ was treated as a

numerical parameter, and its setting was conditioned by the choice of $\lambda$. Here, $L_{sat} \sim 1$, to allow a sufficiently fast crystal growth. Obviously, attachment kinetics limited growth enabling snowflakes faceting is lost. Faceting is yet recovered through the highly anisotropic interface, as suggested by Debierre and Karma in ref. 27. One may argue that no experimental argument alleges for the predominance of surface tension over attachment kinetics in the faceting process.³⁶ However, it was numerically evidenced in ref. 15, that attachment kinetics played only a minor role in ice crystal faceting compared to surface tension, at least in the continuum model. To account for the anisotropy of the water molecule attachment kinetics on snow crystals,⁸,⁹ the anisotropy function $B(n)$ is introduced in the second term of Eq. (2) as well. The vertical/horizontal preference is also considered through anisotropic surface diffusion,¹¹ $\nabla_r = (\partial_{x_i} \partial_{y_j} \Gamma_{\partial_z})$.

Constants $\lambda, W_0, \tau_0$ and $\tilde{D}$ are entangled by the asymptotic analysis mapping the phase field model to Stefan sharp interface model, as shown in ref. 22. This sets $\tilde{D} = 0.6267\lambda$, and $W_0 = d_0\lambda/0.8839$, where $d_0$ is the isotropic capillarity length. This analysis also links the anisotropic interface width $W_0A(n)$ to the characteristic time of interface propagation $\tau_0A(n)^2$, where $\tau_0 = 0.6267\lambda W_0^2$.

The coupling constant was set to $\lambda = 3.0$ as in.²³ The horizontal anisotropy constant was set to $\varepsilon_{xy} = 0.1$ for vertical growth ($\Gamma \ge 1$), and $\varepsilon_{xy} = 0.2$ for horizontal growth ($\Gamma < 1$). To reproduce different snowflakes morphologies, parameters $\Gamma, u_0, L_{sat}$ and $\varepsilon_z$ were varied. All study case parameters are gathered in Table 1.

## RESULTS

Two limit cases of snowflake growth are displayed in Fig. 1. Figure 1a shows the different stages of the growth of a fernlike dendrite 1 (p. 59 of ref. 16), and Fig. 1b details the formation of a hollow prism 1 (p. 64–66 of ref. 16).

For both cases, the growth stages are very similar to the experimental growth kinetics obtained by Libbrecht:³⁷ the snowflake first aligns on the equilibrium Wulff shape, until the Mullins–Sekerka instability³⁸ occurs for a critical crystal radius. This instability is related to the Berg effect,³⁹ stating that the supersaturation field around a faceted snow crystal is largest at facet edges. Destabilization thus occurs when the snowflake reaches a critical radius,⁴⁰ and kinetic effects at corners overcome the highly anisotropic surface tension.⁶ Two examples of snowflake growth are provided at npj computational material’s website (supplementary video 1 and 2).

<table><caption>**Table 1.** Parameters for different morphologies of snowflakes</caption><thead><tr><th></th><th>Γ</th><th>u<sub>0</sub></th><th>L<sub>sat</sub></th><th>ε<sub>z</sub></th></tr></thead><tbody><tr><td>stellar dendrite</td><td>0.5</td><td>0.7</td><td>1.0</td><td>0.05</td></tr><tr><td>fern dendrite</td><td>0.5</td><td>0.8</td><td>1.6</td><td>0.05</td></tr><tr><td>dendritic arms plate</td><td>0.4</td><td>0.6</td><td>1.0</td><td>0.1</td></tr><tr><td>stellar plate</td><td>0.5</td><td>0.5</td><td>1.8</td><td>0.3</td></tr><tr><td>& stars</td><td>0.5</td><td>0.5</td><td>1.0</td><td>0.3</td></tr><tr><td>double plate</td><td>0.4</td><td>0.5</td><td>1.6</td><td>0.1</td></tr><tr><td>sectored plate</td><td>0.4</td><td>0.5(0.8)</td><td>1.0</td><td>0.25(0.1)</td></tr><tr><td>solid plate</td><td>0.25</td><td>0.4</td><td>2.0</td><td>0.2</td></tr><tr><td>rolls on plate</td><td>0.4(0.2)</td><td>0.5(0.8)</td><td>1.0</td><td>0.25(0.4)</td></tr><tr><td>& needles</td><td>3.0</td><td>0.8</td><td>1.0</td><td>0.5</td></tr><tr><td>hollow prism</td><td>3.0</td><td>0.3</td><td>3.0</td><td>0.5</td></tr><tr><td>- capped column</td><td>5.0(0.2)</td><td>0.8</td><td>1.0(2.0/1.0)</td><td>0.5(0.4)</td></tr></tbody></table>

Note: Parenthesis correspond to the case where parameter were changed during simulations
growth satisfies the microscopic solvability theory,³⁰⁻³³ consistently with experiments.³⁴

## MODEL

Snowflakes growth in supersaturated water vapour was simulated using a phase field approach in three dimensions, based on a methodology developed in.²² In this model, two coupled variables $\phi$ and $u$ are considered. $\phi$ is an order parameter referring to the ice (+1) and vapour (-1) phases. The ice/vapour interface is described by a continuous variation of $\phi$, connecting -1 and +1. $u = (c - c'_{sat})/c'_{sat}$ is the reduced supersaturation of water vapour, where $c'_{sat}(T)$ is the saturation number density of vapour above ice, at temperature $T$. At initial time, the reduced supersaturation $u_0$ is homogeneous. The growth kinetics of snowflakes is governed by two non conservative phase field equations. Their adimensionalized form is given by:

$$A(n)^2 \partial_t \phi = -f'(\phi) + \lambda B(n) g'(\phi) u + \frac{1}{2} \nabla_r \cdot \left( |\nabla \phi|^2 \frac{\partial [A(n)^2]}{\partial \nabla \phi} + A(n)^2 \nabla_r \phi \right) \quad (1)$$

$$\partial_t u = \tilde{D} \nabla_r \cdot (q(\phi) \nabla_r u) - \frac{L_{sat}}{2} B(n) \partial_t \phi, \quad (2)$$

where space and time are scaled by the interface width $W_0$, and the relaxation time $\tau_0$, respectively. In the Eq. (1), the double well potential $f(\phi) = -\phi^2/2 + \phi^4/4$ is the free energy density $f_{iv}(c, T', \phi)$ of the ice/vapour system, at temperature $T$, and saturation concentration $c = c'_{sat}$: $f(\phi) = f_{iv}(c = c'_{sat}, T' = T, \phi)$. The second term in the Eq. (1) accounts for the coupling between $u$ and $\phi$, and promotes ice phase growth in supersaturated atmosphere, where $c > c'_{sat}$. Formally, it corresponds to the first order term in the Taylor expansion of the bulk potential, for $c$ in the neighbourhood of $c'_{sat}$. The coupling constant $\lambda$ can thus be computed by $\lambda = c_{sat}/(30H) \partial_c f_{iv}|_{c_{sat}}$, where $H$ is the free energy barrier. $g'(\phi) = (1 - \phi^2)^2$ is an interpolation function introduced in.²² Its form allows to keep the bulk potential minima at $\phi = \pm 1$ for any $u$. To describe the strong anisotropy of snowflakes growth along the vertical axis,⁶ we propose to introduce the kinetic anisotropy function $B(n) = \sqrt{n_x^2 + n_y^2 + \Gamma^2 n_z^2}$, where $n = -\nabla \phi / |\nabla \phi|$ is the unit normal vector of $\phi$. This anisotropic contribution may be justified by the premelted layer growth theory, stating that the primary habit results from different attachment rates of water molecules between prismatic and basal facets.⁸ The parameter $\Gamma > 0$ governs the preference between horizontal and vertical growth, called the primary habit of snowflakes. It can be empirically related to the temperature.³⁵ The addition of the first two terms in the Eq. (1) corresponds to the anisotropic thermodynamic driving force.

The last two terms in the Eq. (1) describe the ice/vapour interface formation and propagation,²² where $A(n)$ is the surface tension anisotropy function. In this work, a new expression of the anisotropy function was derived: $A(n) = 1 + \varepsilon_{xy} \cos(6\theta) + \varepsilon_z \cos(2\psi)$, where $\theta = \arctan(n_y/n_x)$ and $\psi = \arctan(\sqrt{n_x^2 + n_y^2}/n_z)$ are the polar and azimuthal angles, respectively. It accounts for both the horizontal 6-fold symmetry, and the vertical planar symmetry of snowflakes. $\varepsilon_{xy}$ and $\varepsilon_z$ are anisotropy constants.

The Eq. (2) describes one-sided diffusion of water molecules in vapour, and vapour condensation on ice. Diffusion is controlled by two quantities: the function $q(\phi) = 1 - \phi$, which prohibits diffusion within ice, and the reduced diffusion coefficient $\tilde{D} = D\tau_0/W_0^2$, where $D$ is the diffusion coefficient. The second term in the Eq. (2) accounts for conversion of vapour into ice. $L_{sat}$ can thus be interpreted as the rate of water depletion in vapour, via molecule attachment at the interface. In this study, $L_{sat}$ was treated as a

numerical parameter, and its setting was conditioned by the choice of $\lambda$. Here, $L_{sat} \sim 1$, to allow a sufficiently fast crystal growth. Obviously, attachment kinetics limited growth enabling snowflakes faceting is lost. Faceting is yet recovered through the highly anisotropic interface, as suggested by Debierre and Karma in ref. 27. One may argue that no experimental argument alleges for the predominance of surface tension over attachment kinetics in the faceting process.³⁶ However, it was numerically evidenced in ref. 15, that attachment kinetics played only a minor role in ice crystal faceting compared to surface tension, at least in the continuum model. To account for the anisotropy of the water molecule attachment kinetics on snow crystals,⁸,⁹ the anisotropy function $B(n)$ is introduced in the second term of Eq. (2) as well. The vertical/horizontal preference is also considered through anisotropic surface diffusion,¹¹ $\nabla_r = (\partial_{xy}\partial_{yz}\Gamma\partial_z)$.

Constants $\lambda, W_0, \tau_0$ and $\tilde{D}$ are entangled by the asymptotic analysis mapping the phase field model to Stefan sharp interface model, as shown in ref. 22. This sets $\tilde{D} = 0.6267\lambda$, and $W_0 = d_0\lambda/0.8839$, where $d_0$ is the isotropic capillarity length. This analysis also links the anisotropic interface width $W_0A(n)$ to the characteristic time of interface propagation $\tau_0A(n)^2$, where $\tau_0 = 0.6267\lambda W_0^2$.

The coupling constant was set to $\lambda = 3.0$ as in.²³ The horizontal anisotropy constant was set to $\varepsilon_{xy} = 0.1$ for vertical growth ($\Gamma \ge 1$), and $\varepsilon_{xy} = 0.2$ for horizontal growth ($\Gamma < 1$). To reproduce different snowflakes morphologies, parameters $\Gamma, u_0, L_{sat}$ and $\varepsilon_z$ were varied. All study case parameters are gathered in Table 1.

## RESULTS

Two limit cases of snowflake growth are displayed in Fig. 1. Figure 1a shows the different stages of the growth of a fernlike dendrite 1 (p. 59 of ref. 16), and Fig. 1b details the formation of a hollow prism 1 (p. 64–66 of ref. 16).

For both cases, the growth stages are very similar to the experimental growth kinetics obtained by Libbrecht:³⁷ the snowflake first aligns on the equilibrium Wulff shape, until the Mullins–Sekerka instability³⁸ occurs for a critical crystal radius. This instability is related to the Berg effect,³⁹ stating that the supersaturation field around a faceted snow crystal is largest at facet edges. Destabilization thus occurs when the snowflake reaches a critical radius,⁴⁰ and kinetic effects at corners overcome the highly anisotropic surface tension.⁶ Two examples of snowflake growth are provided at npj computational material’s website (supplementary video 1 and 2).

<table><caption>**Table 1.** Parameters for different morphologies of snowflakes</caption><thead><tr><th></th><th>Γ</th><th>u<sub>0</sub></th><th>L<sub>sat</sub></th><th>ε<sub>z</sub></th></tr></thead><tbody><tr><td>stellar dendrite</td><td>0.5</td><td>0.7</td><td>1.0</td><td>0.05</td></tr><tr><td>fern dendrite</td><td>0.5</td><td>0.8</td><td>1.6</td><td>0.05</td></tr><tr><td>dendritic arms plate</td><td>0.4</td><td>0.6</td><td>1.0</td><td>0.1</td></tr><tr><td>stellar plate</td><td>0.5</td><td>0.5</td><td>1.8</td><td>0.3</td></tr><tr><td>& stars</td><td>0.5</td><td>0.5</td><td>1.0</td><td>0.3</td></tr><tr><td>double plate</td><td>0.4</td><td>0.5</td><td>1.6</td><td>0.1</td></tr><tr><td>sectored plate</td><td>0.4</td><td>0.5(0.8)</td><td>1.0</td><td>0.25(0.1)</td></tr><tr><td>solid plate</td><td>0.25</td><td>0.4</td><td>2.0</td><td>0.2</td></tr><tr><td>rolls on plate</td><td>0.4(0.2)</td><td>0.5(0.8)</td><td>1.0</td><td>0.25(0.4)</td></tr><tr><td>& needles</td><td>3.0</td><td>0.8</td><td>1.0</td><td>0.5</td></tr><tr><td>hollow prism</td><td>3.0</td><td>0.3</td><td>3.0</td><td>0.5</td></tr><tr><td>- capped column</td><td>5.0(0.2)</td><td>0.8</td><td>1.0(2.0/1.0)</td><td>0.5(0.4)</td></tr></tbody></table>

Note: Parenthesis correspond to the case where parameter were changed during simulations
<figure><img src="image_1.png" /><figcaption>(a) Fern dendrite 2</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>(b) Hollow prism 13</figcaption></figure>

**Fig. 1** Isosurface representation ($\phi=0$) of snowflakes growth, using the software Blender for visual rendering. Top: fern dendrite 2. The growth mechanism is essentially two dimensional.⁶ First, the initial disc grows into a transient flat faceted hexagon. Then, the branching instability occurs in the horizontal plane ($t=50$), for a critical crystal radius. The resulting faceted dendrites are equipped with developed faceted side branches. Bottom: prism 13. The seed grows vertically into a prism aligned on the Wulff shape. Then the basal facet breaks, and a hollow is formed at $t=10$. Time in $\tau_0$ unit

In Fig. 1a, the simulated fernlike dendrite morphology differs from the non axisymmetric shape predicted by the theory of Brener,³³ usually used for dendrites with cubic symmetry.²⁴, ⁴¹ We suggest this is due to the strong vertical anisotropy flattening, and to faceting, which limits the formation of a horizontal lacuna at the tail of dendrites. The obtained feathery shape with transversally sharp tip rather resembles Furukawa experimental ice crystals.¹⁷ Besides, experimental snowflakes display characteristic surface patterns such as ridges,¹³ and flat basal planes forming steps.¹⁷ Such patterns are also reproduced by our model in Fig. 1a. At low temperatures, surface nucleation and spiral growth are the leading

growth processes on the basal faces for horizontal growth.⁹ Though such mechanisms are not explicitly included in our model, the presence of surface steps in our simulations evokes the terrace growth resulting from nucleation, experimentally observed on ice crystals.¹⁹ It can be underlined that the succession of flat basal planes on both real and simulated snowflakes, reflects the main stages of snowflake history. For instance, the hexagonal Wulff shape before branching instability is clearly memorised in Fig. 1b. Comparison with real snowflakes is provided in the supplementary material 1 at npj computational material's website.
<figure><img src="image_3.png" /><figcaption>(a) Fern dendrite 2</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Hollow prism 13</figcaption></figure>

**Fig. 1** Isosurface representation ($\phi=0$) of snowflakes growth, using the software Blender for visual rendering. Top: fern dendrite 2. The growth mechanism is essentially two dimensional.⁶ First, the initial disc grows into a transient flat faceted hexagon. Then, the branching instability occurs in the horizontal plane ($t=50$), for a critical crystal radius. The resulting faceted dendrites are equipped with developed faceted side branches. Bottom: prism 13. The seed grows vertically into a prism aligned on the Wulff shape. Then the basal facet breaks, and a hollow is formed at $t=10$. Time in $\tau_0$ unit

In Fig. 1a, the simulated fernlike dendrite morphology differs from the non axisymmetric shape predicted by the theory of Brener,³³ usually used for dendrites with cubic symmetry.²⁴, ⁴¹ We suggest this is due to the strong vertical anisotropy flattening, and to faceting, which limits the formation of a horizontal lacuna at the tail of dendrites. The obtained feathery shape with transversally sharp tip rather resembles Furukawa experimental ice crystals.¹⁷ Besides, experimental snowflakes display characteristic surface patterns such as ridges,¹³ and flat basal planes forming steps.¹⁷ Such patterns are also reproduced by our model in Fig. 1a. At low temperatures, surface nucleation and spiral growth are the leading

growth processes on the basal faces for horizontal growth.⁹ Though such mechanisms are not explicitly included in our model, the presence of surface steps in our simulations evokes the terrace growth resulting from nucleation, experimentally observed on ice crystals.¹⁹ It can be underlined that the succession of flat basal planes on both real and simulated snowflakes, reflects the main stages of snowflake history. For instance, the hexagonal Wulff shape before branching instability is clearly memorised in Fig. 1b. Comparison with real snowflakes is provided in the supplementary material 1 at npj computational material's website.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Snowflakes size *L* (*W*₀ unit) vs. time *t* (*τ*₀ unit) (marks), for *u*₀ = 0.6. Two growth regimes appear, with a transition at *t* = 50 for vertical growth, and *t* = 80 for horizontal growth (dots). First regime: *L*(*t*) ∝ *t*ª (*ρ* = 0.99), with *a* = 0.62 for horizontal growth (**line**), and *a* = 0.63 for vertical growth (**dashes**). Second regime: *L*(*t*) ∝ *At* (*ρ* = 0.99), with *A* = 0.12 for horizontal growth (**line**), and *A* = 0.41 for vertical growth (**dashes**). Insert: tip velocity *V* (*W*₀/*τ*₀ unit) vs. reduced supersaturation *u*₀ (marks). *V* = *A*′*u*₀ (*ρ* = 0.99) with *A*′ = 0.62 for horizontal growth (**line**), and *A*′ = 0.62 for vertical growth (**dashes**). **Blue**: horizontal growth (parameters 2 except *u*₀ varied). **Red**: vertical growth (parameters 11 except *u*₀ varied).</figcaption></figure>

Figure 2 displays the time evolution of snowflake size *L*, defined as arm length for horizontal growth (*blue*), and column length for vertical growth (*red*).

It can be seen that two growth regimes appear. After a first growth regime, when both dynamics are similar, the needle growth accelerates, while plate-like growth slows down. In the case of horizontal growth, Libbrecht experimentally noted in ref. 42, that this change of regime coincides with full faceting occurrence on the prismatic face, causing growth to decelerate. In the case of vertical growth, it was also observed in ref. 42, that the growth acceleration was related to the formation of a vicinal surface at the needle tip, which enhances water attachment at the tip. During the first regime, the crystal growth is slower than diffusion, and the system satisfies the Laplacian approximation. Within this framework, Algrem et al. found in ref. 43, that arm growth should display the self similar scaling behaviour $L(t)\propto t^{3/5}$. This is confirmed by our simulations. It is interesting to note that the disparity of our simulations with the $t^{1/2}$ diffusive law of Zener et al. for spherical precipitates,⁴⁴ is due to surface tension.

In the second regime, the growth kinetics is linear (**line** and **dashes** in Fig. 2), consistently with the selection theory.¹⁸, ³² The associated tip velocity *V* is thus the slope of this linear function. Using the same procedure for different *u*₀, *V* vs. *u*₀ could be plotted in the insert. It appears that *V* is proportional to *u*₀. Therefore, the growth velocity of both faceted and non-faceted dendrites satisfies the universal law.³⁰ It is consistent with previous experimental,³⁴ and numerical observations.¹⁵, ²⁷ A similar relation was observed during faceted growth in other materials, including silicon.⁴⁵ For further discussion, we refer to ref. 46. Langer et al. argued in ref. 30. that the frontmost tip of dendrites was molecularly rough for both non-faceted and faceted dendrites. The attachment kinetics can thus be considered isotropic with circular symmetry.

DISCUSSION

Figure 3 compiles equilibrium shapes of simulated snowflakes, using the parameters in Table 1. An achievement of our model is

to recover the principal snowflakes morphologies of the Nakaya diagram,⁵ depending exclusively on four parameters: $\Gamma$, $u_0$, $L_{sat}$ and $\varepsilon_z$.

The first parameter $\Gamma$ fully determines the primary habit in our model. We suggest this parameter may be fitted on the dimensionless attachment coefficient $a^{\text{kin}}(T, u_{\text{surf}})$, encompassing the molecular kinetics at the solid/vapour interface.⁴⁷ Here, $u_{\text{surf}}$ is the supersaturation at the surface of the snowflake. Recent developments suggest $a^{\text{kin}}$ can be parametrised by $A(T) \exp(-u_{\infty}(T)/u_{\text{surf}})$, where the relation between $A$ and $u_{\text{surf}}$ and the temperature is empirically known.³⁵ The coefficient $a^{\text{kin}}$ accounts for the alternation of dominating growth mechanisms with the temperature. However, the physical interpretation of this dependence is still unclear.⁶ To explain this behaviour, two theories are proposed: surface nucleation rates for low temperatures,⁹ and quasi-liquid layer growth with different wettability for basal and prismatic facets⁴⁸ near the melting point. In this work, the link between $\Gamma$ and atmospheric parameters was simplified. First, the variation of $a^{\text{kin}}$ with the surface supersaturation was neglected, as fitting $\Gamma$ to the temperature only was sufficient to match the Nakaya diagram in a first approximation. Second, only two values of $\Gamma$ were used: $\Gamma \sim 0.4$ for horizontal growth, and $\Gamma = 3.0$ for vertical growth. This rough phenomenological setting of $\Gamma$ reduces to $\Gamma < 1$ or $\Gamma > 1$ depending on temperature. Then, Hertz-Knudsen relation for facet normal velocity¹⁸ makes supersaturation $u_0$ the priming parameter in facet destabilization and growth velocity.⁹ In the phase field model, it is reinterpreted as the thermodynamic driving force, and higher values also favour branching instability and faster growth.²¹ The parameter $L_{sat}$ controls the compactness of branching for non faceted growth.²⁶ Here, low values foster branching instability, while raising $L_{sat}$ drives the system toward quasistatic diffusion, and fosters a stronger faceting. As for the parameter $\varepsilon_z$, it mainly influences the formation of surface patterns for horizontal growth.

Especially interesting are two model examples. First, the double plate 1 (p. 75 in ref. 16), results from the sandwich instability of a flat snowflake, and the twelve arms 1 was simulated by the aggregation of two simple stars, with a 30° tilt.²⁰

Finally, simulated snowflakes morphology dependence on temperature and density excess over vapour/water equilibrium $\rho_{\text{super}}^W = \rho - \rho_{\text{sat}}^W$ (g m⁻³), is quantitatively consistent with the Nakaya diagram (see ref. 49). Using $u_0 = (\rho_{\text{super}}^W + \rho_{\text{sat}}^W - \rho_{\text{sat}})/\rho_{\text{sat}}$, where is the vapour/ice equilibrium density, and fitting $\Gamma$ on ref. 47, (Libbrecht, K. G. Physical dynamics of ice crystal growth. Annu. Rev. Mater. Res. (2017), Unpublished), our simulations predict plate formation near ($T, \rho_{\text{super}}^W$) = (−19, 0.15) and (−12, 0.4), vs. (−19, 0.15), (−12, 0.4), but also (−10, 0.18) in the Nakaya diagram. Simulated dendrites occur near (−16, 0.4), vs. (−15, 0.4) and (−15, 0.25) for experiments. Finally, hollow columns are formed at (−8, 0.6) in our simulations, vs. (−8, 0.6) and (−8, 0.3) in the Nakaya diagram. However, small supersaturations $\rho_{\text{sat}}^l < \rho < \rho_{\text{sat}}^W$ are beyond reach for our model. This shortcoming is complementary to the model of Barrett et al.¹⁵. Indeed, the latter is bound to $\rho \le \rho_{\text{sat}}^W$ or $T \to 0$ °C, due to the Laplacian approximation for small supersaturations.⁶ Our model on the contrary requires larger values for $L_{sat}$ and greater computational resources to reach such supersaturations.

CONCLUSION

In this paper, we have shown that the proposed modified phase field model is able to reproduce the complex dynamics of snowflake growth. Benefiting from the modularity and versatility of the phase field model, this approach can be extended to describe the growth of snowflakes in real atmospheric conditions. It can notably be equipped with stochastic effects, for the asymmetric growth of more exotic morphologies of snowflakes in the Nakaya diagram. Continuum fluid dynamics can also be easily
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Snowflakes size *L* (*W*₀ unit) vs. time *t* (*τ*₀ unit) (marks), for *u*₀ = 0.6. Two growth regimes appear, with a transition at *t* = 50 for vertical growth, and *t* = 80 for horizontal growth (dots). First regime: *L*(*t*) ∝ *t*ª (*ρ* = 0.99), with *a* = 0.62 for horizontal growth (**line**), and *a* = 0.63 for vertical growth (**dashes**). Second regime: *L*(*t*) ∝ *At* (*ρ* = 0.99), with *A* = 0.12 for horizontal growth (**line**), and *A* = 0.41 for vertical growth (**dashes**). Insert: tip velocity *V* (*W*₀/*τ*₀ unit) vs. reduced supersaturation *u*₀ (marks). *V* = *A*′*u*₀ (*ρ* = 0.99) with *A*′ = 0.62 for horizontal growth (**line**), and *A*′ = 0.62 for vertical growth (**dashes**). **Blue**: horizontal growth (parameters 2 except *u*₀ varied). **Red**: vertical growth (parameters 11 except *u*₀ varied)</figcaption></figure>

Figure 2 displays the time evolution of snowflake size *L*, defined as arm length for horizontal growth (*blue*), and column length for vertical growth (*red*).

It can be seen that two growth regimes appear. After a first growth regime, when both dynamics are similar, the needle growth accelerates, while plate-like growth slows down. In the case of horizontal growth, Libbrecht experimentally noted in ref. 42, that this change of regime coincides with full faceting occurrence on the prismatic face, causing growth to decelerate. In the case of vertical growth, it was also observed in ref. 42, that the growth acceleration was related to the formation of a vicinal surface at the needle tip, which enhances water attachment at the tip. During the first regime, the crystal growth is slower than diffusion, and the system satisfies the Laplacian approximation. Within this framework, Algrem et al. found in ref. 43, that arm growth should display the self similar scaling behaviour $L(t)\propto t^{3/5}$. This is confirmed by our simulations. It is interesting to note that the disparity of our simulations with the $t^{1/2}$ diffusive law of Zener et al. for spherical precipitates,⁴⁴ is due to surface tension.

In the second regime, the growth kinetics is linear (**line** and **dashes** in Fig. 2), consistently with the selection theory.¹⁸, ³² The associated tip velocity *V* is thus the slope of this linear function. Using the same procedure for different *u*₀, *V* vs. *u*₀ could be plotted in the insert. It appears that *V* is proportional to *u*₀. Therefore, the growth velocity of both faceted and non-faceted dendrites satisfies the universal law.³⁰ It is consistent with previous experimental,³⁴ and numerical observations.¹⁵, ²⁷ A similar relation was observed during faceted growth in other materials, including silicon.⁴⁵ For further discussion, we refer to ref. 46. Langer et al. argued in ref. 30. that the frontmost tip of dendrites was molecularly rough for both non-faceted and faceted dendrites. The attachment kinetics can thus be considered isotropic with circular symmetry.

DISCUSSION

Figure 3 compiles equilibrium shapes of simulated snowflakes, using the parameters in Table 1. An achievement of our model is

to recover the principal snowflakes morphologies of the Nakaya diagram,⁵ depending exclusively on four parameters: $\Gamma$, $u_0$, $L_{sat}$ and $\varepsilon_z$.

The first parameter $\Gamma$ fully determines the primary habit in our model. We suggest this parameter may be fitted on the dimensionless attachment coefficient $a^{\text{kin}}(T, u_{\text{surf}})$, encompassing the molecular kinetics at the solid/vapour interface.⁴⁷ Here, $u_{\text{surf}}$ is the supersaturation at the surface of the snowflake. Recent developments suggest $a^{\text{kin}}$ can be parametrised by $A(T) \exp(-u_{\infty}(T)/u_{\text{surf}})$, where the relation between $A$ and $u_{\text{surf}}$ and the temperature is empirically known.³⁵ The coefficient $a^{\text{kin}}$ accounts for the alternation of dominating growth mechanisms with the temperature. However, the physical interpretation of this dependence is still unclear.⁶ To explain this behaviour, two theories are proposed: surface nucleation rates for low temperatures,⁹ and quasi-liquid layer growth with different wettability for basal and prismatic facets⁴⁸ near the melting point. In this work, the link between $\Gamma$ and atmospheric parameters was simplified. First, the variation of $a^{\text{kin}}$ with the surface supersaturation was neglected, as fitting $\Gamma$ to the temperature only was sufficient to match the Nakaya diagram in a first approximation. Second, only two values of $\Gamma$ were used: $\Gamma \sim 0.4$ for horizontal growth, and $\Gamma = 3.0$ for vertical growth. This rough phenomenological setting of $\Gamma$ reduces to $\Gamma < 1$ or $\Gamma > 1$ depending on temperature. Then, Hertz-Knudsen relation for facet normal velocity¹⁸ makes supersaturation $u_0$ the priming parameter in facet destabilization and growth velocity.⁹ In the phase field model, it is reinterpreted as the thermodynamic driving force, and higher values also favour branching instability and faster growth.²¹ The parameter $L_{sat}$ controls the compactness of branching for non faceted growth.²⁶ Here, low values foster branching instability, while raising $L_{sat}$ drives the system toward quasistatic diffusion, and fosters a stronger faceting. As for the parameter $\varepsilon_z$, it mainly influences the formation of surface patterns for horizontal growth.

Especially interesting are two model examples. First, the double plate 1 (p. 75 in ref. 16), results from the sandwich instability of a flat snowflake, and the twelve arms 1 was simulated by the aggregation of two simple stars, with a 30° tilt.²⁰

Finally, simulated snowflakes morphology dependence on temperature and density excess over vapour/water equilibrium $\rho_{\text{super}}^W = \rho - \rho_{\text{sat}}^W$ (g m⁻³), is quantitatively consistent with the Nakaya diagram (see ref. 49). Using $u_0 = (\rho_{\text{super}}^W + \rho_{\text{sat}}^W - \rho_{\text{sat}})/\rho_{\text{sat}}$, where is the vapour/ice equilibrium density, and fitting $\Gamma$ on ref. 47, (Libbrecht, K. G. Physical dynamics of ice crystal growth. Annu. Rev. Mater. Res. (2017), Unpublished), our simulations predict plate formation near ($T, \rho_{\text{super}}^W$) = (−19, 0.15) and (−12, 0.4), vs. (−19, 0.15), (−12, 0.4), but also (−10, 0.18) in the Nakaya diagram. Simulated dendrites occur near (−16, 0.4), vs. (−15, 0.4) and (−15, 0.25) for experiments. Finally, hollow columns are formed at (−8, 0.6) in our simulations, vs. (−8, 0.6) and (−8, 0.3) in the Nakaya diagram. However, small supersaturations $\rho_{\text{sat}}^l < \rho < \rho_{\text{sat}}^W$ are beyond reach for our model. This shortcoming is complementary to the model of Barrett et al.¹⁵. Indeed, the latter is bound to $\rho \le \rho_{\text{sat}}^W$ or $T \to 0$ °C, due to the Laplacian approximation for small supersaturations.⁶ Our model on the contrary requires larger values for $L_{sat}$ and greater computational resources to reach such supersaturations.

CONCLUSION

In this paper, we have shown that the proposed modified phase field model is able to reproduce the complex dynamics of snowflake growth. Benefiting from the modularity and versatility of the phase field model, this approach can be extended to describe the growth of snowflakes in real atmospheric conditions. It can notably be equipped with stochastic effects, for the asymmetric growth of more exotic morphologies of snowflakes in the Nakaya diagram. Continuum fluid dynamics can also be easily
$\Gamma < 1$

<figure><img src="image_4.png" /><figcaption>1 - stellar dendrite</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>2 - fern dendrite</figcaption></figure>

$\Gamma \ge 1$

<figure><img src="image_6.png" /><figcaption>11 - needle</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>3 - dendrite plate</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>4 - stellar plate</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>12 - needles cluster</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>5 - simple star*</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>6 - 12-star*</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>13 - hollow prism</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>7 - double plate</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>8 - sectored plate</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>14 - capped col. a</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>9 - solid plate</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>10 - scrolls on plate</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>15 - capped col. b</figcaption></figure>

**Fig. 3** Simulated snowflakes classification, using the software Blender for visual rendering. $\Gamma < 1$ gives the flat snow crystals 1 to 9. $\Gamma \ge 1$ gives columns 11 & 13. Primary habit switch is achieved by shifting $\Gamma$ during simulations. $\Gamma$ from 3.0 to 0.2 at $t = 70$, gives capped columns 14 and 15. $\Gamma$ from 0.4 to 3.0 gives scrolls on plate 10. For $\Gamma < 1$, no facet breaking happens when $u_0 \le 0.4$ (solid plate 9). Branching instability appears for $u_0 = 0.5$ (snowflakes 3 to 8), and side branching requires $u_0 \ge 0.7$ (dendrites 1 and 2). For $\Gamma > 1$, $u_0 = 0.3$ gives a prism 13, and $u_0 = 0.8$ induces vertical branching with needles 11 and 12. $L_{sat} = 1.0$ gives a thick side branching for the stellar dendrite 1. Rising $L_{sat}$ to $L_{sat} = 1.6$ refines side branches with the fern dendrite 1. Higher values of $L_{sat}$ destroy side branching, and lead to branches maximum faceting with the stellar plate 4 ($L_{sat} = 1.8$), (*): see Method
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Simulated snowflakes classification, using the software Blender for visual rendering. $\Gamma < 1$ gives the flat snow crystals 1 to 9. $\Gamma \ge 1$ gives columns 11 & 13. Primary habit switch is achieved by shifting $\Gamma$ during simulations. $\Gamma$ from 3.0 to 0.2 at $t = 70$, gives capped columns 14 and 15. $\Gamma$ from 0.4 to 3.0 gives scrolls on plate 10. For $\Gamma < 1$, no facet breaking happens when $u_0 \le 0.4$ (solid plate 9). Branching instability appears for $u_0 = 0.5$ (snowflakes 3 to 8), and side branching requires $u_0 \ge 0.7$ (dendrites 1 and 2). For $\Gamma > 1$, $u_0 = 0.3$ gives a prism 13, and $u_0 = 0.8$ induces vertical branching with needles 11 and 12. $L_{sat} = 1.0$ gives a thick side branching for the stellar dendrite 1. Rising $L_{sat}$ to $L_{sat} = 1.6$ refines side branches with the fern dendrite 1. Higher values of $L_{sat}$ destroy side branching, and lead to branches maximum faceting with the stellar plate 4 ($L_{sat} = 1.8$), (*): see Method</figcaption></figure>
introduced in Eq. (1) and Eq. (2), to simulate the air flow around the snowflake during its fall from clouds.⁵⁰ This could notably answer the need expressed by meteorology and aircraft safety for new metrologies, able to estimate ice water content in clouds.⁵¹, ⁵² Besides, the developed approach is not restricted to snowflake growth. It opens a new way to answer numerous outstanding questions concerning the general issue of faceted dendritic growth in a wide range of materials. We can cite, for instance, the faceted dendritic growth observed in pure isotactic polystyrene films, or the snowflake-like growth of graphene monocrystal.⁵³

## METHODS

Phase field simulations were performed using the Fourier-spectral semi-implicit scheme with periodic boundary conditions. For horizontal growth, simulations were performed with grid spacing $Δx = Δy = Δz = 0.8$, and time step $Δt = 0.05$,²³ on a $400 × 400 × 64$ simulation box ($800 × 800 × 96$ for (*) in Fig. 3). Vertical growth simulations required a more precise discretization. Grid spacing was thus reduced to $Δx = Δy = Δz = 0.4$, and $Δt = 0.01$ as in ref. 23, on a $128 × 128 × 256$ simulation box. Snow crystal growth was initiated by a circular-disc shape germ ($φ = 1$) of radius $R = 8Δx$, within water vapour ($φ = -1$) of homogeneous reduced supersaturation $u_0 > 0$. A new regularisation method for high interfacial energy leading to crystal missing orientations was also adapted from two dimensions,²⁸ to three dimensions. This allowed to overcome the restriction for the critical values of the anisotropy constants and choose $ε_{xy} > 1/35$ and $ε_z > 1/3$, as required to achieve faceting.²¹ Additional information on the method is provided in the supplementary material 2 at npj computational material's website.

## ACKNOWLEDGEMENTS

The authors would like to thank warmly Prof. K.G. Libbrecht of the California Institute of Technology, for his interest in this study, and for providing us with his instructive review (Libbrecht, K. G. Physical dynamics of ice crystal growth. Annu. Rev. Mater. Res. (2017), Unpublished). The simulations were performed at the Centre de Ressources Informatiques de Haute-Normandie (CRIHAN) and at the IDRIS of CNRS. The authors thank the French National Agency for financial support through the CISTIC project of programme Investissements d'Avenir LabEx EMC3 (ANR-10-LABX-09-01).

## COMPETING INTERESTS

The authors declare that they have no competing interests.

## REFERENCES

1. Pamuk, B. et al. Anomalous nuclear quantum effects in ice. *Phys. Rev. Lett.* **108**, 193003 (2012).

2. Ball, P. Material witness: close to the edge. *Nat. Mater.* **15**, 1060–1060 (2016).

3. Nakaya, U., Satô, I. & Sekido, Y. Preliminary experiments on the artificial production of snow crystals. *J. Fac Sci. Hokkaido Imp. Uni. Ser. 2*, 1–11 (1938).

4. Magono, C. & Chung, W. Meteorological classification of natural snow crystals. *J. Fac Sci. Hokkaido Uni. Ser. 7 Geophys* **2**, 321–335 (1966).

5. Nakaya, U. *The Formation of Ice Crystals*, 207–220 (American Meteorological Society, 1951).

6. Libbrecht, K. G. The physics of snow crystals. *Rep. Prog. Phys.* **68**, 855 (2005).

7. Mason, B. J. Snow crystals, natural and man made. *Contemp. Phys.* **33**, 227–243 (1992).

8. Kuroda, T. & Lacmann, R. Growth kinetics of ice from the vapour phase and its growth forms. *J. Cryst. Growth* **56**, 189–205 (1982).

9. Nelson, J. Growth mechanisms to explain the primary and secondary habits of snow crystals. *Philos. Mag.* **A 81**, 2337–2373 (2001).

10. Nada, H. & Furukawa, Y. Anisotropic growth kinetics of ice crystals from water studied by molecular dynamics simulation. *J. Cryst. Growth* **169**, 587–597 (1996).

11. Furukawa, Y. & Nada, H. Anisotropic surface melting of an ice crystal and its relationship to growth forms. *J. Phys. Chem. B* **101**, 6167–6170 (1997).

12. Benet, J., Llombart, P., Sanz, E. & MacDowell, L. G. Premelting-induced smoothening of the ice-vapor interface. *Phys. Rev. Lett.* **117**, 096101 (2016).

13. Gravner, J. & Griffeath, D. Modeling snow-crystal growth: A three-dimensional mesoscopic approach. *Phys. Rev. E* **79**, 011601 (2009).

14. Kelly, J. G. & Boyer, E. C. Physical improvements to a mesoscopic cellular automaton model for three-dimensional snow crystal growth. *Cryst. Growth Des.* **14**, 1392–1405 (2014).

15. Barrett, J. W., Garcke, H. & Nürnberg, R. Numerical computations of faceted pattern formation in snow crystal growth. *Phys. Rev. E* **86**, 011604 (2012).

16. Libbrecht, K. G. *Ken Libbrecht's Field Guide to Snowflakes* (Voyageur Press, 2006).

17. Furukawa, Y. & Shimada, W. Three-dimensional pattern formation during growth of ice dendrites relation to universal law of dendritic growth. *J. Cryst. Growth* **128**, 234–239 (1993).

18. Saito, Y. *Statistical Physics of Crystal Growth*, vol. 2 (World Scientific, 1996).

19. Sazaki, G., Zepeda, S., Nakatsubo, S., Yokoyama, E. & Furukawa, Y. Elementary steps at the surface of ice crystals visualized by advanced optical microscopy. *Proc. Natl. Acad. Sci* **107**, 19702–19707 (2010).

20. Kikuchi, K. & Uyeda, H. Formation mechanisms of multibranched snow crystals (twelve-, eighteen-, twenty-four-branched crystals). *Atmos. Res.* **474**, 169–179 (1998).

21. Singer-Loginova, I. & Singer, H. The phase field technique for modeling multi-phase materials. *Rep. Prog. Phys.* **71**, 106501 (2008).

22. Karma, A. & Rappel, W.-J. Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics. *Phys. Rev. E* **53**, R3017 (1996).

23. Ramirez, J., Beckermann, C., Karma, A. & Diepers, H.-J. Phase-field modeling of binary alloy solidification with coupled heat and solute diffusion. *Phys. Rev. E* **69**, 051607 (2004).

24. Cartalade, A., Younsi, A. & Plapp, M. Lattice boltzmann simulations of 3d crystal growth: Numerical schemes for a phase-field model with anti-trapping current. *Comput. Math. Appl.* **71**, 1784–1798 (2016).

25. Barrett, J. W., Garcke, H. & Nürnberg, R. On the stable discretization of strongly anisotropic phase field models with applications to crystal growth. *ZAMM Z. Angew. Math. Mech.* **93**, 719–732 (2013).

26. Kobayashi, R. Modeling and numerical simulations of dendritic crystal growth. *Physica D* **63**, 410–423 (1993).

27. Debierre, J.-M., Karma, A., Celestini, F. & Guérin, R. Phase-field approach for faceted solidification. *Phys. Rev. E* **68**, 041604 (2003).

28. Eggleston, J. J., McFadden, G. B. & Voorhees, P. W. A phase-field model for highly anisotropic interfacial energy. *Physica D* **150**, 91–103 (2001).

29. Libbrecht, K. G. Toward a comprehensive model of snow crystal growth dynamics: 1. overarching features and physical origins. Preprint at arXiv:1211.5555 (2012).

30. Langer, J., Sekerka, R. & Fujioka, T. Evidence for a universal law of dendritic growth rates. *J. Cryst. Growth* **44**, 414–418 (1978).

31. Langer, J. & Müller-Krumbhaar, H. Theory of dendritic growth. elements of a stability analysis. *Acta Metall* **26**, 1681–1687 (1978).

32. Amar, M. B. & Brener, E. Theory of pattern selection in three-dimensional non-axisymmetric dendritic growth. *Phys. Rev. Lett.* **71**, 589 (1993).

33. Brener, E. Three-dimensional dendritic growth. *J. Cryst. Growth* **166**, 339–346 (1996).

34. Libbrecht, K. G., Crosby, T. & Swanson, M. Electrically enhanced free dendrite growth in polar and non-polar systems. *J. Cryst. Growth* **240**, 241–254 (2002).

35. Libbrecht, K. G. & Rickerby, M. E. Measurements of surface attachment kinetics for faceted ice crystal growth. *J. Cryst. Growth* **377**, 1–8 (2013).

36. Libbrecht, K. G. On the equilibrium shape of an ice crystal. Preprint at arXiv:1205.1452 (2012).

37. Libbrecht, K. G. Snowcrystals.com. http://www.snowcrystals.com.

38. Mullins, W. W. & Sekerka, R. Stability of a planar interface during solidification of a dilute binary alloy. *J. Appl. Phys.* **35**, 444–451 (1964).

39. Berg, W. Crystal growth from solutions. *Proc. R. Soc. A* **916**, 79–95 (1938).

40. Fujioka, T. & Sekerka, R. Morphological stability of disc crystals. *J. Cryst. Growth* **24**, 84–93 (1974).

41. Singer, H. & Bilgram, J. Three-dimensional reconstruction of xenon dendrites. *Europhys. Lett.* **68**, 240 (2004).

42. Libbrecht, K. G. Toward a comprehensive model of snow crystal growth: 5. measurements of changes in attachment kinetics from background gas interactions at -5 c. Preprint at arXiv:1602.08528 (2016).

43. Almgren, R., Dai, W.-S. & Hakim, V. Scaling behavior in anisotropic hele-shaw flow. *Phys. Rev. Lett.* **71**, 3461 (1993).

44. Zener, C. Theory of growth of spherical precipitates from solid solution. *J. Appl. Phys.* **20**, 950–953 (1949).

45. Yang, X. et al. Crystal and faceted dendrite growth of silicon near (100). *Acta Mater.* **60**, 3259–3267 (2012).

46. Wang, Y. & Khachaturyan, A. G. Shape instability during precipitate growth in coherent solids. *Acta Metall.* **43**, 1837–1857 (1995).

47. Libbrecht, K. G. Crystal growth in the presence of surface melting and impurities: an explanation of snow crystal growth morphologies. Preprint at arXiv:0810.0689 (2008).
introduced in Eq. (1) and Eq. (2), to simulate the air flow around
the snowflake during its fall from clouds.⁵⁰ This could notably
answer the need expressed by meteorology and aircraft safety for
new metrologies, able to estimate ice water content in clouds.⁵¹, ⁵²
Besides, the developed approach is not restricted to snowflake
growth. It opens a new way to answer numerous outstanding
questions concerning the general issue of faceted dendritic
growth in a wide range of materials. We can cite, for instance,
the faceted dendritic growth observed in pure isotactic polystyr-
ene films, or the snowflake-like growth of graphene
monocrystal.⁵³

METHODS

Phase field simulations were performed using the Fourier-spectral semi-
implicit scheme with periodic boundary conditions. For horizontal growth,
simulations were performed with grid spacing Δx = Δy = Δz = 0.8, and time
step Δt = 0.05,²³ on a 400 × 400 × 64 simulation box (800 × 800 × 96
for (*) in Fig. 3). Vertical growth simulations required a more precise
discretization. Grid spacing was thus reduced to Δx = Δy = Δz = 0.4, and
Δt = 0.01 as in ref. 23, on a 128 × 128 × 256 simulation box. Snow crystal
growth was initiated by a circular-disc shape germ (φ = 1) of radius R = 8Δx,
within water vapour (φ = −1) of homogeneous reduced supersaturation
u₀ > 0. A new regularisation method for high interfacial energy leading to
crystal missing orientations was also adapted from two dimensions,²⁸ to three
dimensions. This allowed to overcome the restriction for the critical values of
the anisotropy constants and choose εₓᵧ > 1/35 and εₓₓ > 1/3, as required to
achieve faceting.²¹ Additional information on the method is provided in the
supplementary material 2 at npj computational material's website.

ACKNOWLEDGEMENTS

The authors would like to thank warmly Prof. K.G. Libbrecht of the California Institute of Technology, for his interest in this study, and for providing us with his instructive review (Libbrecht, K. G. Physical dynamics of ice crystal growth. Annu. Rev. Mater. Res. (2017), Unpublished). The simulations were performed at the Centre de Ressources Informatiques de Haute-Normandie (CRIHAN) and at the IDRIS of CNRS. The authors thank the French National Agency for financial support through the CISTIC project of programme Investissements d'Avenir LabEx EMC3 (ANR-10-LABX-09-01).

COMPETING INTERESTS

The authors declare that they have no competing interests.

REFERENCES

1. Pamuk, B. et al. Anomalous nuclear quantum effects in ice. *Phys. Rev. Lett.* **108**, 193003 (2012).

2. Ball, P. Material witness: close to the edge. *Nat. Mater.* **15**, 1060–1060 (2016).

3. Nakaya, U., Satô, I. & Sekido, Y. Preliminary experiments on the artificial pro-
duction of snow crystals. *J. Fac Sci. Hokkaido Imp. Uni. Ser. 2*, 1–11 (1938).

4. Magono, C. & Chung, W. Meteorological classification of natural snow crystals. *J. Fac Sci. Hokkaido Uni. Ser. 7 Geophys* 2, 321–335 (1966).

5. Nakaya, U. *The Formation of Ice Crystals*, 207–220 (American Meteorological Society, 1951).

6. Libbrecht, K. G. The physics of snow crystals. *Rep. Prog. Phys.* **68**, 855 (2005).

7. Mason, B. J. Snow crystals, natural and man made. *Contemp. Phys.* **33**, 227–243 (1992).

8. Kuroda, T. & Lacmann, R. Growth kinetics of ice from the vapour phase and its growth forms. *J. Cryst. Growth* **56**, 189–205 (1982).

9. Nelson, J. Growth mechanisms to explain the primary and secondary habits of snow crystals. *Philos. Mag.* **A 81**, 2337–2373 (2001).

10. Nada, H. & Furukawa, Y. Anisotropic growth kinetics of ice crystals from water studied by molecular dynamics simulation. *J. Cryst. Growth* **169**, 587–597 (1996).

11. Furukawa, Y. & Nada, H. Anisotropic surface melting of an ice crystal and its relationship to growth forms. *J. Phys. Chem. B* **101**, 6167–6170 (1997).

12. Benet, J., Llombart, P., Sanz, E. & MacDowell, L. G. Premelting-induced smoothening of the ice-vapor interface. *Phys. Rev. Lett.* **117**, 096101 (2016).

13. Gravner, J. & Griffeath, D. Modeling snow-crystal growth: A three-dimensional mesoscopic approach. *Phys. Rev. E* **79**, 011601 (2009).

14. Kelly, J. G. & Boyer, E. C. Physical improvements to a mesoscopic cellular automaton model for three-dimensional snow crystal growth. *Cryst. Growth Des.* **14**, 1392–1405 (2014).

15. Barrett, J. W., Garcke, H. & Nürnberg, R. Numerical computations of faceted pattern formation in snow crystal growth. *Phys. Rev. E* **86**, 011604 (2012).

16. Libbrecht, K. G. *Ken Libbrecht's Field Guide to Snowflakes* (Voyageur Press, 2006).

17. Furukawa, Y. & Shimada, W. Three-dimensional pattern formation during growth of ice dendrites relation to universal law of dendritic growth. *J. Cryst. Growth* **128**, 234–239 (1993).

18. Saito, Y. *Statistical Physics of Crystal Growth*, vol. 2 (World Scientific, 1996).

19. Sazaki, G., Zepeda, S., Nakatsubo, S., Yokoyama, E. & Furukawa, Y. Elementary steps at the surface of ice crystals visualized by advanced optical microscopy. *Proc. Natl. Acad. Sci* **107**, 19702–19707 (2010).

20. Kikuchi, K. & Uyeda, H. Formation mechanisms of multibranched snow crystals (twelve-, eighteen-, twenty-four-branched crystals). *Atmos. Res.* **474**, 169–179 (1998).

21. Singer-Loginova, I. & Singer, H. The phase field technique for modeling multi-phase materials. *Rep. Prog. Phys.* **71**, 106501 (2008).

22. Karma, A. & Rappel, W.-J. Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics. *Phys. Rev. E* **53**, R3017 (1996).

23. Ramirez, J., Beckermann, C., Karma, A. & Diepers, H.-J. Phase-field modeling of binary alloy solidification with coupled heat and solute diffusion. *Phys. Rev. E* **69**, 051607 (2004).

24. Cartalade, A., Younsi, A. & Plapp, M. Lattice boltzmann simulations of 3d crystal growth: Numerical schemes for a phase-field model with anti-trapping current. *Comput. Math. Appl.* **71**, 1784–1798 (2016).

25. Barrett, J. W., Garcke, H. & Nürnberg, R. On the stable discretization of strongly anisotropic phase field models with applications to crystal growth. *ZAMM Z. Angew. Math. Mech.* **93**, 719–732 (2013).

26. Kobayashi, R. Modeling and numerical simulations of dendritic crystal growth. *Physica D* **63**, 410–423 (1993).

27. Debierre, J.-M., Karma, A., Celestini, F. & Guérin, R. Phase-field approach for faceted solidification. *Phys. Rev. E* **68**, 041604 (2003).

28. Eggleston, J. J., McFadden, G. B. & Voorhees, P. W. A phase-field model for highly anisotropic interfacial energy. *Physica D* **150**, 91–103 (2001).

29. Libbrecht, K. G. Toward a comprehensive model of snow crystal growth dynamics: 1. overarching features and physical origins. Preprint at arXiv:1211.5555 (2012).

30. Langer, J., Sekerka, R. & Fujioka, T. Evidence for a universal law of dendritic growth rates. *J. Cryst. Growth* **44**, 414–418 (1978).

31. Langer, J. & Müller-Krumbhaar, H. Theory of dendritic growth. elements of a stability analysis. *Acta Metall* **26**, 1681–1687 (1978).

32. Amar, M. B. & Brener, E. Theory of pattern selection in three-dimensional non-axisymmetric dendritic growth. *Phys. Rev. Lett.* **71**, 589 (1993).

33. Brener, E. Three-dimensional dendritic growth. *J. Cryst. Growth* **166**, 339–346 (1996).

34. Libbrecht, K. G., Crosby, T. & Swanson, M. Electrically enhanced free dendrite growth in polar and non-polar systems. *J. Cryst. Growth* **240**, 241–254 (2002).

35. Libbrecht, K. G. & Rickerby, M. E. Measurements of surface attachment kinetics for faceted ice crystal growth. *J. Cryst. Growth* **377**, 1–8 (2013).

36. Libbrecht, K. G. On the equilibrium shape of an ice crystal. Preprint at arXiv:1205.1452 (2012).

37. Libbrecht, K. G. Snowcrystals.com. http://www.snowcrystals.com.

38. Mullins, W. W. & Sekerka, R. Stability of a planar interface during solidification of a dilute binary alloy. *J. Appl. Phys.* **35**, 444–451 (1964).

39. Berg, W. Crystal growth from solutions. *Proc. R. Soc. A* **916**, 79–95 (1938).

40. Fujioka, T. & Sekerka, R. Morphological stability of disc crystals. *J. Cryst. Growth* **24**, 84–93 (1974).

41. Singer, H. & Bilgram, J. Three-dimensional reconstruction of xenon dendrites. *Europhys. Lett.* **68**, 240 (2004).

42. Libbrecht, K. G. Toward a comprehensive model of snow crystal growth: 5. measurements of changes in attachment kinetics from background gas interactions at-5 c. Preprint at arXiv:1602.08528 (2016).

43. Almgren, R., Dai, W.-S. & Hakim, V. Scaling behavior in anisotropic hele-shaw flow. *Phys. Rev. Lett.* **71**, 3461 (1993).

44. Zener, C. Theory of growth of spherical precipitates from solid solution. *J. Appl. Phys.* **20**, 950–953 (1949).

45. Yang, X. et al. Crystal and faceted dendrite growth of silicon near (100). *Acta Mater.* **60**, 3259–3267 (2012).

46. Wang, Y. & Khachaturyan, A. G. Shape instability during precipitate growth in coherent solids. *Acta Metall.* **43**, 1837–1857 (1995).

47. Libbrecht, K. G. Crystal growth in the presence of surface melting and impurities: an explanation of snow crystal growth morphologies. Preprint at arXiv:0810.0689 (2008).
48. Meuler, A. J. et al. Relationships between water wettability and ice adhesion. *ACS Appl. Mater. Interfaces* **2**, 3100–3110 (2010).

49. Marshall, J. S. & Langleben, M. P. A theory of snow-crystal habit and growth. *J. Meteorol* **11**, 104–120 (1954).

50. Jeong, J.-H., Goldenfeld, N. & Dantzig, J. A. Phase field model for three-dimensional dendritic growth with fluid flow. *Phys. Rev. E* **64**, 041602. doi:10.1103/PhysRevE.64.041602 (2001).

51. Porcheron, E. et al. Development of a spectrometer for airborne measurement of droplet sizes in clouds. *J. Eur. Opt. Soc. Rapid. Publ.* **10**, 15030-1-15030-6 (2015).

52. Kielar, J. J. et al. Simultaneous interferometric in-focus and out-of-focus imaging of ice crystals. *Opt. Commun.* **372**, 185-195 (2016).

53. Whiteway, E., Yang, W., Yu, V. & Hilke, M. Time evolution of the growth of single graphene crystals and high resolution isotope labeling. Preprint at arXiv:1509.01579 (2015).

<figure><img src="image_19.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017

Supplementary Information accompanies the paper on the npj Computational Materials website (doi:10.1038/s41524-017-0015-1).
48. Meuler, A. J. et al. Relationships between water wettability and ice adhesion. *ACS Appl. Mater. Interfaces* **2**, 3100–3110 (2010).

49. Marshall, J. S. & Langleben, M. P. A theory of snow-crystal habit and growth. *J. Meteorol* **11**, 104–120 (1954).

50. Jeong, J.-H., Goldenfeld, N. & Dantzig, J. A. Phase field model for three-dimensional dendritic growth with fluid flow. *Phys. Rev. E* **64**, 041602. doi:10.1103/PhysRevE.64.041602 (2001).

51. Porcheron, E. et al. Development of a spectrometer for airborne measurement of droplet sizes in clouds. *J. Eur. Opt. Soc. Rapid. Publ.* **10**, 15030-1-15030-6 (2015).

52. Kielar, J. J. et al. Simultaneous interferometric in-focus and out-of-focus imaging of ice crystals. *Opt. Commun.* **372**, 185-195 (2016).

53. Whiteway, E., Yang, W., Yu, V. & Hilke, M. Time evolution of the growth of single graphene crystals and high resolution isotope labeling. Preprint at arXiv:1509.01579 (2015).

<figure><img src="image_7.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017

Supplementary Information accompanies the paper on the npj Computational Materials website (doi:10.1038/s41524-017-0015-1).