A single snowflake has over a quintillion (10¹⁸) water molecules! How do molecules arrange themselves into such beautiful and intricate shapes? Is it possible to model snowflakes in computer code, on a pretty basic laptop?

From a few simple rules emerges the complexity of a snowflake. Two 'forces' are at work on the molecules of water vapour: the pulling strength at the surface of the snowflake, and pushing pressure of the vapour around the snowflake. If the surface is rough and bumpy, as opposed to smooth, there is a stronger attractive force to pull molecules into the nooks and crannies, to smooth it out. And if there is more vapour pushing, then more molecules bump up against the surface.

The 6 points of a snowflake will grow faster - and in unison - because they stick further out into the space around the snowflake where the vapour has higher saturation. And the spaces in between the points will fill in to smooth out the bumps.

My aim is to faithfully emulate the natural physics as much as I can. I want to generate the hexagonal symmetry naturally, not by coying one point to the 5 others. I want to generate a unique snowflake every time. But I've had to cut a few corners ... because a full model would be too computationally complex on my PC.

My second aim is to explore how different conditions reliably grow different types of crystal, based on the scientify models. At this stage, I have a working model that can generate different shapes of snow crystal, very reminiscent of real snow crystal, however I haven't been able to model real-world conditions like temperature and humidity. That is for future work.

My starting point is the science of snow crystals, developed by Ken Libbrecht and others. Crystallization is actually an extremely complex physical process that involves chaotic - but statistical - molecular reactions at the interface between a frozen surface and a vapour surface (some of which are not well understood yet). Despite this, there's a simple formula for growth.

(Summarized from Libbrecht's page on Snowflake Science)

A note on terminology: 'Snowflake' is the everyday name for just about anything that falls from the sky in the winter, whereas 'snow crystal' is the more specific technical name for an individual crystal of ice, what we normally call a snowflake.

A snow crystal forms as the water vapour around a seed ice crystal high in the clouds freezes onto the seed.

Snow crystals can form as plates, dendrites, prisms, needles, and columns, but always with hexagonal symmetry. Variations in temperature and humidity causes different types of snow crystals in the wild.

The temperature and humidity surrounding the crystal determine if molecules will freeze onto it. As a snow crystal falls through the clouds, the conditions change unpredictably, making the arms grow differently. But at any one instant the conditions around all six arms are the same, so each arm grows the same. Since every path through the clouds is different, every snow crystal is unique.

Snow crystals have a hexagonal, or six-fold, symmetry because of the angular 3D geometry of water molecules (H₂O). This means new water molecules can only attach in certain way .

Two factors determine how water molecules freeze and attach to the growing crystal:

• The surface physics of the crystal itself, including the surface roughness.
• The vapour saturation near the surface of the crystal.

Variations in these factors affect the velocity of crystal growth normal to each point on the surface, through a simple formula (see below).

Below is a basic description of the algorithm. See the code on github for full details.

Setup

Growth phase

• During the growth phase, check every vapour cell on the crystal surface (ie next to already crystallized cells). Each such potential attachment point will have a growth velocity. If the velocity is high enough, then freeze the cell onto the crystal.
• Calculate the velocity using a simplified version of Libbrecht's key formula:

  velocity = alpha x sigma

  where alpha is the Condensation Coefficient and sigma is the Vapour Saturation, at that potential attachment point.
• The Condensation Coefficient (alpha) is based on the number of surrounding already condensed cells (more surrounding ice means a strongger attractive force).
• The Vapour Saturation (sigma) is the percentage saturation of the cell (a higher saturation means higher pressure, and more water molecules bumping onto the crystal).
• Initially, the Vapour Saturation of every vapour cell is set to a uniform initial value. When a cell condenses, 'capture' (ie remove) the vapour from the cell and adjoining cells, thus taking water out of the air and into the ice, reducing its saturation. If the saturation is too low, because it's all been captured in the snow crystal, then a cell cannot condense (because the velocity will be low).

Diffusion phase

• In the diffusion phase, move vapour saturation from highly saturated cells to lower saturated cells over the whole vapour array, simulating how the vapour pressure drives to equalize saturation.

Changing conditions

• Model the changing conditions of temperature and humidity using a single Conditions variable, and factor it into the sigma calculation of the velocity formula. After each iteration, modify Conditions by a random walk, up or down, simulating unpredicable temperature and humidity changes.
• There are several user-specified parameters:
  • Attachment_Power: A constant factor that affects the contribution of the Condensation Coefficient (alpha) to the velocity.
  • Velocity_Limit: A minimum threshold on velocity (expressed as a percentage of the maximum velocity in the iteration) for a cell to condense into ice.
  • Vapour_Reduction_Factor: The vapour saturation left in the surroundings cell, after a cell condenses into ice, capturing some of the vapour.