When cracks form in a thin contracting layer, they sequentially break the layer into smaller and smaller pieces. A rectilinear crack pattern encodes information about the order of crack formation, as later cracks tend to intersect with earlier cracks at right angles. In a hexagonal pattern, in contrast, the angles between all cracks at a vertex are near 120°. Hexagonal crack patterns are typically seen when a crack network opens and heals repeatedly, in a thin layer, or advances by many intermittent steps into a thick layer. Here, it is shown how both types of pattern can arise from identical forces, and how a rectilinear crack pattern can evolve towards a hexagonal one. Such an evolution is expected when cracks undergo many opening cycles, where the cracks in any cycle are guided by the positions of cracks in the previous cycle but when they can slightly vary their position and order of opening. The general features of this evolution are outlined and compared with a review of the specific patterns of contraction cracks in dried mud, polygonal terrain, columnar joints and eroding gypsum–sand cements.
Cracks in cooling or drying media can form captivating patterns of connected networks, such as the artistic craquelure patterns sometimes seen in pottery glazes to those found in dried mud, or the polygonal networks covering the polar regions of the Earth and Mars. Two types of patterns are common. The first is a rectilinear pattern, such as one that may form when a homogeneous slurry is dried or cooled uniformly [1–3]. The vertices of these patterns typically contain one crack path intersecting another straight crack at right angles, a T-junction. Second are crack networks with a more regular hexagonal pattern, such as the columnar joints of the Giant’s Causeway in Northern Ireland, UK [4,5]. These are dominated by cracks intersecting at 120° or Y-junctions. Curiously, examples of both extreme types of patterns can sometimes be found in the same systems, such as garden-variety mud cracks. These, and other examples, are shown in figure 1. Hexagonally ordered patterns may form as the result of either desiccation or thermal contraction cracks, and there is a deep physical symmetry between these two mechanisms [8,9]. Furthermore, they may be classified into crack patterns that have evolved in both time and space (such as a network of crack tips that delimit the prismatic forms of columnar joints as they advance ), and those that have recurred in the same place but evolved over time as the pattern repeatedly healed and recracked [3,11].
Here, it is shown that both hexagonal and rectilinear crack patterns can arise naturally from the same assumptions of fracture behaviour, and how a rectilinear, T-junction-dominated pattern can develop into a hexagonal pattern, with Y-junctions. This paper is roughly divided into two parts. In the first part (§§2 and 3), the relevant physics of linear elastic fracture mechanics in a thin brittle layer is summarized, and used to develop a model for the evolution of a rectilinear crack pattern into a hexagonal one, through the cyclic opening of cracks. The second part (§§4–7) reviews the ordering of crack patterns in a number of geophysical systems: desiccation cracks in clays that are repeatedly wetted and dried; the thermal contraction cracks in permafrost, known as polygonal terrain; columnar joints in lava and starch; and cracks in eroding gypsum sand. These patterns are each explored and shown to obey the same general physical mechanism of crack pattern evolution. In the case of desiccation cracks in clays, new experiments are also reported, which show that the length scale of the evolution of a rectilinear mud-crack pattern towards a hexagonal one scales linearly with the thickness of the cracking layer.
2. Contraction cracks in thin layers
When an elastic body cools or dries, it tends to shrink and may crack. Here, the case of a brittle film under tensile stress and firmly attached to a much thicker, non-cracking substrate, as shown in figure 2a, is reviewed. As a drying body loses water, it develops an internal pore pressure as a result of capillary effects. In the absence of cracks, the equations of stress equilibrium and poro-elasticity are easily solved and predict an in-plane stress, σ0, proportional to this capillary pressure. Similarly, if a film is cooled, it will have an in-plane stress proportional to the temperature change, and to the mismatch in thermal expansivities of the film and substrate. The similarity between the patterns of cracks that develop during drying and cooling reflects the underlying similarity between the physics of poro-elasticity and thermo-elasticity [8,9].
Tension is released when a crack opens in the brittle layer, as the normal stresses in the film must vanish on the new crack surfaces. The substrate will resist deformation, however, and exert a restoring force on the film–substrate interface, approximately proportional to the local in-plane strain [12,13]. For a long straight crack in a thin film, as shown in figure 2b, these conditions give rise to an exponential dependence of the height-averaged film stress on the distance x away from the crack, 2.1where σ⊥ and σ∥ are the stresses normal and parallel to the crack face, and ν1 is the Poisson ratio of the film . The length scale over which stresses relax, l, is proportional to the film thickness, h, 2.2The order-one term g is given numerically in  and depends on the Dundurs parameters α and β, ratios involving the plane-strain elastic moduli E1 and E2 and Poisson ratios ν1 and ν2 of the film and substrate. If there is no elastic mismatch between the cracking layer and its substrate, l=2.0 h, while for a compliant film on an unyielding substrate l is between 1.2 and 1.3 h, depending on ν1. In the limit of very soft (or fluid) substrates, l diverges, and the film behaves as an unconstrained layer.
Once a crack has formed, the Griffith criterion compares the energy required to create new fracture surface, Gc, with the rate of change of elastic energy in the cracked body for an infinitesimal extension of the crack, G, and predicts crack growth when more energy would be released than consumed . While the first term, the critical energy release rate, is a material parameter, the second term depends on the geometry of the crack and its local surroundings. For an isolated straight film crack that is long compared with the film thickness [14,12], the strain energy release rate 2.3Thus, when the film stress reaches some critical value , the film can crack. As Gc is a material constant, thinner layers will require higher stresses to crack than thicker layers, a fact well established in the engineering literature [16–18].
As stress develops, it can be relieved by the sequential fragmentation of the film. This may be illustrated by considering an array of parallel cracks, each separated by a distance λ that is large compared with l. If each crack relaxes stress as per equation (2.1), then the energy release rate, per crack, of the array is 2.4where σc is the stress required to open a single isolated crack or the critical cracking stress . As drying or cooling proceeds, the background stress σ0 increases, and the allowed crack spacing λ will shrink. Counterintuitively, this process does not continue indefinitely, no matter how high the film stress gets. When the cracks are closer than a few film thicknesses to each other, the stresses are no longer well approximated by a through-thickness-average stress. Numerical simulations of the stress state between two parallel cracks have shown that, for cracks closer than some limiting spacing λc∼h, the surface stress in the region between the two cracks is compressive and will tend to close any further flaws or potential cracks [19,20]. For cracks in a compliant layer on a hard substrate, experiments have confirmed that the crack spacing is limited to about three times the film thickness [1,3,21–23]. It is likely that, in other cases, the saturation spacing scales with the stress relaxation length l, although this has not been rigorously demonstrated. If so, then where there is no difference in the elastic properties of the cracking layer and its substrate, the limiting crack spacing should be approximately five times the depth of the cracks.
When multiple cracks form, they will interact with each other and pattern the film. Rather than the simple parallel cracks just considered, cracks will normally nucleate randomly on material defects or along the faces of pre-existing cracks and grow until they hit something, such as a boundary or another crack. Each crack will tend to grow in the direction that maximizes G−Gc, i.e. the difference between the strain energy released and the fracture energy consumed during growth . This is a local principle, as the energy release rates are calculated at the point and time of the crack growth. As in equation (2.1), any pre-existing crack is more efficient at relieving stress perpendicular to its face than that parallel to it. Therefore, a growing crack will tend to curve to intersect another crack at right angles to it. As shown in figure 2c, the resulting T-junctions thus encode time information that can be used to partially reconstruct the order in which cracks form .
3. A model for evolving crack patterns
Cracks in a homogeneous film will tend to pattern the surface into roughly equal-sized rectilinear pieces, bounded with right-angled T-junctions. However, many natural crack patterns, such as columnar joints, show hexagonal planforms and equi-angled Y-junctions. In both cases, Euler’s formula implies that the average number of neighbours for any polygon is six, as long as cracks do not cross each other, or vertices form from more than three cracks. It is possible to continuously deform a rectilinear crack network into a hexagonal one (allowing vertices to pass each other, if necessary ), if the pattern is allowed to evolve. Development towards a hexagonal pattern might thus be expected by noting that a hexagonal tiling minimizes the ratio of the crack-length (perimeter) to area of the polygons . However, as discussed in the previous section, a growing crack tip can only respond to the local strain energy release rate at the position and moment where it is opening. There is no general way to link this local energy maximization principle to a global constraint. Here, instead, it is shown that evolution from a rectilinear to hexagonal pattern can be explained by entirely local responses of advancing cracks, as the result of meeting three conditions: (i) that cracks recur or advance repeatedly, (ii) that the previous positions of cracks act as lines of weakness, guiding the next iteration of cracking, and (iii) that the order of opening of cracks may change in each iteration.
In order to develop from a rectilinear network of cracks into a hexagonal one, the pattern must have some means of changing with time. Polygonal terrain consists of networks of thermal contraction cracks that open during the cold polar winters and which thaw, and heal, during the summers; the crack pattern reforms each year . For columnar joints, the pattern of growing cracks at any particular moment is confined to a thin region, approximating a two-dimensional film . Each crack advances stepwise, and the columns are the forms delimited by the planar crack network as it advances through space.
Even if a layer cracks and heals repeatedly, the pattern will not necessarily evolve. If the system has no memory of previous patterns, each cycle will lead to a novel, but statistically similar, rectilinear crack pattern. If there is memory of vertex positions only, perhaps through fixed defects (e.g. a pebble in a mud puddle), then cracks may still nucleate along the same points in each generation, but then grow to intersect with each other at right angles, again creating a statistically similar pattern each time they open. If any memory is too strong, however, the crack pattern may simply repeat in each cycle without change. To evolve away from a rectilinear pattern, new cracks must be guided, but imperfectly, by previous cracks.
A nascent crack can be guided by a crack in an earlier cycle in a number of ways. During columnar joint formation, each crack advance follows the crack above it, which concentrates stress along its tip, but can grow downwards in a slightly different direction . In a thin layer, in contrast, the previous cycle’s crack may not have fully healed or may have damaged their vicinity. In either case, the critical energy release rate Gc is lower near old crack paths . Alternatively, if the previously cracked region is more compliant (lower E1) than uncracked material, then near a previous crack, the strain energy release rate G will be higher, as suggested by the form of equation (2.3). Finally, as with polygonal terrain, the crack may leave a depression after healing [11,7], locally lowering the average height h by δ. In this case, the total energy spent by a unit length of crack extension would be reduced to Gc(h−δ). However, the rate of energy release, G, remains integrated throughout the region surrounding the crack of width ∼h and will be relatively unchanged. It will thus be energetically favourable for a crack to be confined within the depression of the former crack. The action of any, or a combination, of these factors means that a crack tip advancing to maximize G−Gc will be guided to follow a previously existing crack path as it grows.
Changes in the order in which cracks approach and leave vertices allow for the above conditions to evolve a crack pattern. The case of a crack guided by a network of partially healed cracks is demonstrated in figure 3a–f. In this example, consider a film that cracks and heals regularly, but where the regions near previous cracks are more compliant than their surroundings. In the first generation, as shown in figure 3a, one crack has advanced from left to right, while a later crack has curved in from below to intersect it at right angles. After healing, in the next cycle of cracking, a crack is guided along the path of a previous crack. As it approaches the area near the previous vertex, the stress field around its tip becomes asymmetric (i.e. the shear stress intensity factor is non-zero). As has been recently demonstrated, a crack path is attracted by an inclusion of more compliant material . If this attraction is weak, the crack will follow close to its original path and the pattern will have no opportunity to evolve. However, as shown in figure 3b, the asymmetry may cause the crack to be attracted downwards and the curve to follow a different path from its first incarnation. A further crack approaching the vertex in this cycle, as in figure 3c, will also be deflected downwards slightly, as it curves to intersect the first crack. As variations in film stresses can only be felt over distances of a few film thicknesses (equation (2.1)), these changes in crack paths will initially only be expected in a small region near the vertex, whose size should vary in proportion to l, and thus h. Within this region, the new crack paths are deflected, and the angles in which the cracks approach the vertex will be more equi-angled than the original T-junction.
If there are repeated changes in the order in which cracks approach a vertex, then changes in the vertex position will accumulate over many cycles and lead to a directed drift of the vertex in the direction of the crack that originally approached it at right angles. For example, if, in any cycle, a crack approaches first from the lower branch of the vertex shown in figure 3d, the path near the vertex will be unstable, in that any slight perturbation would cause it to deflect to either the left or right branch. If the right branch is chosen, the crack will further round the vertex and move it further down. Over many cycles, if all branches near a vertex crack first with equal frequency, then, by symmetry, the vertex will stabilize with three equal angles of 120°, as shown in figure 3e. However, as shown in figure 3f, the order of cracking in any particular cycle will give rise to deviations in the immediate vicinity of vertices, where cracks curve slightly to intersect each other. Alternatively, this might be notable by a damage zone or through crack tip splittings near vertices.
A pattern which evolves as it passes through space, such as columnar joints, can be thought of as the structure formed by stacking evolving crack networks over many cycles, as shown in figure 3g.
This model of crack ordering makes several predictions. First, that any such hexagonal pattern must have cracks recurring repeatedly, near older cracks, but not always opening in the same direction—and, equivalently, that the order in which cracks arrive around any vertex must change with time. It implies that any environment with a hexagonal crack pattern should also be able to support, and have previously shown, a rectilinear one. Both types of pattern should scale in the same manner, with crack spacings a few times the stress relaxation length l∼h, as long as sufficient stress exists to saturate the spacing. The hexagonal patterns must have had a chance to evolve over many crack opening cycles. The reverse, that rectilinear patterns are ‘young’, is not necessarily true, as a crack pattern may heal sufficiently between openings to erase any memory and reset the pattern. Finally, intermediate patterns should link the rectilinear and hexagonal forms. The effects of an evolving pattern should be first notable by small crack deflections in the vicinities of vertices, which will tend to equalize the angles between cracks. Further from the vertices, in the early stages of evolution, cracks should still appear to be approaching as a T-junction. During this time, the vertices themselves should also show a directed drift in the direction of the cracks that originally approached each vertex at two right angles.
The remaining sections of this article will show how the model conditions and predictions outlined here are realized for desiccation cracks in soils, polygonal terrain, columnar joints and eroding gypsum–sand cements.
4. Mud cracks
Desiccation cracking in thin elastic layers, such as clays or other granular media, is a simple model system that is frequently used to study crack patterns. Experiments in thin starch layers have been used to show how the timing information of a crack pattern is encoded in the shape of T-junctions , for example, while the vibrating of pastes prior to drying has been used to show how crack patterns are influenced by anisotropy in film properties [28,29]. In Nature, mud cracks can be found with both rectilinear and hexagonal (figure 1h) cells. It has been demonstrated experimentally that mud-crack patterns evolve when they are repeatedly moistened and dried . Over only a few cycles of cracking and healing, the crack angles approach 120° and vertices move. In this section, these experiments are further developed in an attempt to look at how the evolution of mud-crack patterns scales with film thickness.
Slurries of bentonite clay (Acros Organics, Bentonite K-10) were prepared by mixing 1 part bentonite with 1.9 parts deionized water (Millipore) by weight, stirring the mixture to be homogeneous and then immediately pouring into large flat-bottomed glass Petri dishes (diameter 185 or 195 mm and height 15 mm). The slurries were dried 60 cm under a 500 W halogen heat lamp in a well-ventilated room, and were fully dry within 48 h. Previous work  shows that, in this system, the crack spacing is about three times the layer thickness, except in layers thinner than approximately 2 mm, when the critical layer thickness is approached. After each drying phase was complete, the dishes were imaged by a digital camera (Nikon D5100). Markings made around the circumference of the dishes were used as fixed points to adjust the magnification, position and rotation of each image to match the configuration of the first drying’s image. Typically images could be matched to within 2 pixels or 0.1 mm accuracy. Each dish was then re-wetted, by spraying with a light mist of deionized water, until it was visibly damp throughout the layer (a mass of water between 1.1 and 1.3 times that of the dry clay was added), and then set back to re-dry.
For the first two drying cycles, the positions of vertices were measured by hand from digital images and compared with 10 dishes between 2.6 and 6.9 mm thick. Two types of vertices were investigated: those where cracks had deflected around a vertex and those where cracks had formed as in the first cycle. Examples of both types are shown in figure 4a,b. The morphology of deflected vertices was sometimes variable. In the thicker layers, especially, there were cases of branching or en passant  cracks near vertices (figure 4c–d) alongside more smoothly curving intersections (figure 4e). For vertices where cracks were deflected off their original path in the second drying cycle, the deflection scaled linearly with the film thickness, as shown in figure 4f. Although these displacements showed large standard deviations, the mean displacement over many such vertices was 14±1% of h. For the remaining vertices, there was a smaller average deflection. However, this was above the limiting resolution of 0.1 mm and scaled with film thickness, at 5±1% of h.
For four dishes, between 2.6 and 6.3 mm thick, the evolution of the crack pattern was monitored over 10 wetting and drying cycles. Although it was originally intended to observe these for longer, the spray bottle used for re-wetting was replaced in the 11th cycle, and the new spray was sufficiently hard to dramatically modify the pattern. For the first 10 cycles of drying, images of the crack patterns were processed automatically (figure 5a–d) by thresholding greyscale images based on local intensity and gradient (edge detection). These images were then skeletonized, and vertices were found as the branch-points of the thresholded image skeleton. Using a Matlab-based particle tracking code , the motion of the vertices was followed. As with previous experiments , cracks could sometimes disappear for several cycles, only to reappear close to their original location. For this analysis, only vertices that were present in all 10 crack patterns were analysed. After this time, the average vertex had moved between about 0.5 and 1 mm, depending on the layer thickness. As shown in figure 5e, when scaled with the layer thickness h, the motion from all patterns collapses onto a single curve, which also overlaps data published previously  for bentonite clay 3.1 mm thick. In all cases, there was a significant change in the average vertex positions within the first three or four cycles, followed by a much more gradual vertex motion. Previously , it was observed that 25–35% of vertices changed the order of formation in any particular generation. The slowdown in vertex motion thus appears to be identifiable with the time it takes for all vertices to change order of opening once. Furthermore, the average displacement after four cycles, 15±1% of h, agrees with the vertex shift, shown in figure 4f, for the subset of cracks which changed the order of arrival at a vertex between the first and second drying cycles.
The angle at which cracks approached vertices was also measured, automatically. An annulus was drawn around each vertex of the skeletonized image, and the points at which the cracks crossed the boundaries of this annulus were found. An inner radius of 5 pixels was used to minimize the effects of any errors in determining the vertex position precisely. An outer radius of half the layer thickness was used to observe the local changes in crack angles near the vertices, which should scale with h. As shown in figure 6, there is a gradual but systematic equalization of the crack angles towards 120° over many cracking cycles. This trend agrees in all experiments.
Finally, figure 5f shows that the direction of the displacement of vertices between the first and tenth drying cycles is well correlated with the direction of the crack that originally approached the junction at the smallest angles to its neighbours. In other words, as predicted in §3, the vertex is moving preferentially along the direction required to stretch a T-junction into a Y-junction, the direction of the crack that originally approached the vertex at two near-right angles.
5. Polygonal terrain
Polygonal terrain consists of networks of thermal contraction cracks, typically a few tens of metres apart, which form in ice-cemented soils . It is common in terrestrial circumpolar regions. The thermal stresses arising from the intense winters of these regions are sufficient to open cracks a few centimetres wide and several metres deep [11,32,33]. While open, loose soil, snow and wind-born detritus may fall into the crack. Thus, the crack does not heal completely when thawed and is likely to reopen in the following winter. The resulting patterns are similar to evolving mud-crack patterns, where each seasonal thermal cycle is equivalent to one cycle of drying and wetting. However, the net addition of material into cracks creates wedges of ice or sand that grow at rates of the order of a millimetre per year and lead to a slow turnover of the near-surface soil [7,34,35].
Recently, perhaps as a result of the landings of Phoenix  and Curiosity  on Martian polygonal terrain, much work has gone into analysing polygonal crack patterns via satellite imaging. Large datasets have been collected for sites on the Earth and Mars [38–41]. These data have been used to demonstrate that polygonal terrain shows some statistical trends that are also found in foam networks  and that the vertex positions contain spatial correlations . Work characterizing these networks is ongoing and interesting. Comparisons with ground observations have shown, however, that very-high-resolution images (approx. 0.25 m pixel−1) are required to correctly measure the features of polygon vertices [39,41], making detailed interpretation of satellite images difficult.
More precise results have been gathered through a few long-term field experiments into the dynamics of polygonal terrain. At Illisarvik, in Arctic Canada, Mackay drained a small permafrost lake and studied the development of new ice wedges over the subsequent 20 years [33,35,42]. The initial crack pattern was strongly influenced by local features, but the most common vertex type was a rectilinear T-junction . There were rapid changes in the ice wedges within the first few years, before ceasing owing to vegetation growth .
In contrast to this, Mackay also established field sites on Garry Island to study well-developed crack patterns . At both locations, he installed cables across the positions of ice wedges and measured the time of crack opening as the time at which the cables broke, during several consecutive years [32,33,35,42,43]. Not all cracks opened each year, yet cracks could reappear in the same spot after several-year absences . This is similar to drying experiments in clays, where cracks sometimes reappeared after a hiatus of several drying cycles . The cracks themselves appear to initiate between the position of the buried ice wedge and the soil surface [35,42]. The interpretation of this is that the wedge itself acts as a memory in this system, as it is preserved over many generations, and is composed of a foreign material that is weaker, or more compliant, than the surrounding soil [11,35,42]. The cracks advanced slowly, in many cases taking several days to travel from one vertex to the next . It was noted that the cracks could shift positions by several centimetres, in different years . The cracks often showed en passant-style doubled paths , like the mud cracks shown in figure 4d, and other branching structures. It was also seen that cracks initiated at different points, and opened in differing orders in different years .
Another long-term experiment was initiated by Black and Berg in the Dry Valleys of Antarctica [34,44] and developed later by Sletten et al. . At several sites, pairs of metal rods were driven into the permafrost, bracketing the cracks of polygonal terrain, and the separation of these rods was carefully measured. An intermittent series of measurements, spanning over 40 years, has been made on the separation of these rods [7,34,44], which shows that the intrusive wedges associated with these cracks are growing at the rate of approximately 0.5 mm yr−1. The observed sites include a range of ages at New Harbour/Taylor Valley (1–2 ka; figures 1d and 7a–c), Victoria Valley (10–12 ka; figure 7d) and Beacon Valley ( Ma; figure 1e), which also display different patterns . In older locations, the crack pattern is dominated by Y-junctions, while at New Harbour it is closer to a rectilinear patten . Based on the relative ages of these sites and their degree of ordering, it was suggested that an initially rectilinear pattern of polygonal terrain cracks evolves towards a more hexagonal pattern over thousands of years .
In the vicinity of the New Harbour site, cracks show additional features that agree with the means of evolution proposed here. As shown in figure 7a, there are regions around many vertices where the cracks curve suddenly to meet the vertex positions, as was shown in figure 2c,d or shown in mud cracks in figure 4b,e. The deviations are consistent with a gradual vertex motion in the direction of the crack that had originally approached at two right angles, similar to that shown for drying mud in figure 5f. For polygonal terrain, this motion should be recorded in an observable asymmetry of the wedge material beneath the crack pattern; the deposits on the side that was originally the primary branches of a T-junction should be thicker than the deposits on the side that had approached at two right angles. This prediction could be checked by excavating a statistically significant number of vertices.
Finally, an opportunity was given to revisit the Black and Berg sites in December 2007. At New Harbour, some vertices showed a complex inner structure, which appears to record how active crack branches curve slightly to intersect each other during different orders of opening in different years, as shown in figure 7b. Here, it was also noted that cracks at two of the 44 pairs of rods installed in 1962 had shifted position to now lie outside the pairs of marker rods (figure 7c). As with Mackay’s more direct measurements of crack positions in Arctic polygons, this indicates that the crack positions can shift.
Taken together, the observations of polygonal terrain show that the fracture network is persistent but that the pattern does change over time. The persistence is related to the fact that cracks leave permanent records of their existence within the soil, which preferentially nucleate and guide later generations of cracks. The order of crack opening also varies, as, for example, not all cracks are active in each year. All these features were seen in mud-crack patterns, and were essential to the model of crack pattern evolution developed here.
6. Columnar joints
Columnar joints in lava have attracted scientific curiosity for centuries [4,5]. The same pattern can be observed in other media, such as quenched glass , thermally shocked sandstone  and vitrified ice . It also appears in desiccated starch, where it has been particularly well studied [5,24,25,45,48–55]. The geometry of the columns in starch is statistically identical to that in lava [24,56]. Columns in starch are, however, usually 0.1–1 cm in diameter, while columns in lava are typically of the order of 10–100 cm across.
In lava, columnar joints form by the gradual extension of thermal contraction cracks into a cooling body, with the cracks growing perpendicular to the isothermal surfaces at the time of cracking . At any particular time, the crack tips represent a network in a thin region near the surface defined by the glass transition temperature isotherm  and the columnar forms can be interpreted as a record of this network as it evolved in time, while propagating through space. Lava may cool, and thus columns may grow, in any direction; in many cases, both an upper and lower colonnade are present in a single flow unit, as the lava cools from both its top and bottom surfaces [58,59]. In starch, columns grow normal to the surfaces of constant moisture content [25,60] and track the position of the so-called funicular–pendular transition, where the water transport mechanisms of capillary flow and vapour diffusion cross over in efficiency .
As cracks advance in lava, they do not do so continuously, but intermittently, leaving striations on the sides of columns, as shown in figure 8e,h or sketched in figure 9a. Each stria represents a single growth increment: once a crack tip is activated (from some weak point or vertex), it advances rapidly forwards through the brittle lava and also laterally following near the path of the edge of the crack face’s previous increment. In the forward direction of column growth, however, the tip blunts and stops as it intrudes into a warmer, more plastic zone. As a result, each crack advance leaves a broad thin mark on the side of a column . These striae record the cooling conditions of the lava at the time of their formation and can be interpreted accordingly [26,56,62,63]. In particular, the laterally averaged height of the striae represents the distance between two isotherms in the cooling lava body and can be used to measure the cooling rate at the time of fracture .
Striae also record the direction of lateral crack propagation during each small forward advance, in their plumose structure . From this, one can tell whether a crack moved towards, or away from, any vertex in any particular growth increment. As has been seen elsewhere [26,63,64], and as shown in figure 8h, this direction can change along subsequent striae. At one colonnade near Bridge Lake, a site described in , the directions of growth of 835 striae were measured. Some example growth direction sequences are shown in figure 8i. It was found that 36±2% of the time a stria differed from the direction of its predecessor. While each stria direction is not randomly chosen, the sequence at which cracks are activated at any vertex changes frequently during a column’s formation.
The crack pattern at the surface of lava lakes observed in Hawaii are rectilinear in appearance and dominated by T-junctions [26,65]. However, the evolution of column vertices was observed at the Boiling Pots, Hawaii, where columns with an average size of approximately 20 cm  were seen to develop from T-junctions to Y-junctions, over a distance of 1 m . Under careful inspection, most vertices were classified there as ‘pseudo-Y’ junctions, as the cracks curved slightly in the immediate vicinity of the vertex , similar to the sketch shown in figure 3f. This small deflection was used to argue that cracks did not all interact at a vertex synchronously . In other words, during each incremental crack advance, there was one crack that arrived at a vertex after the others, causing this slight asymmetry. Similar vertices can be observed in many colonnades; an example from Staffa is shown in figure 8a–d. In starch colonnades, an evolution from a disordered surface pattern to a hexagonal mature pattern is also observed . For columns 1–2 mm across, the pattern orders over a distance of approximately 1 cm from the drying surface ; this ordering region can be seen in figure 1b. Although these observations are sparse, it seems that, in both starch and lava, the hexagonal columnar pattern is well ordered after 5–10 column diameters of growth.
Thus, columns can evolve by slight variations of the crack plane in each crack increment , or by individual cracks overshooting or rounding a vertex formed in the previous cracking cycle . In general, each crack increment is guided by its predecessor, which concentrates stress along its edge, but can advance in a slightly different direction. In starch, this evolution has been followed through X-ray tomography [24,48,53,55] and magnetic resonance imaging . In both cases, the evolution can be followed for very long times. Over such times, the topology of the network also evolves, and by only three distinct mechanisms: two vertices may collide and exchange position; crack faces can fail to propagate, merging two (or more) columns; and new columns can initiate from existing vertices . For neighbouring columns of different sizes, there is also a strong tendency to equalize the column area . Taken together, these observations suggest long-term dynamics that are not covered by the model presented here.
Next, to consider the scaling of columns, we now need to evaluate the effective thickness of the elastic layer in which they form. During the solidification of a lava flow, water or moisture in the environment can act as an efficient transporter of heat, and the thermal transport problem is approximately that of an arbitrarily thick layer, cooled entirely to 100°C, bound to a purely diffusive region over which the temperature transits to that of the internal melt [56,61,66]. Cooling proceeds by the propagation of this solidification front into the lava, with all isotherms travelling together at some velocity v, which can be inferred from the average striae heights . In the reference frame of the solidification front, the linear advection–diffusion problem has a steady-state temperature profile, T, that solves 6.1where Pe=vL/D is the Péclet number. The ratio of diffusion D to advection speed v defines a length scale which is similar to the film thickness h from §2. As shown in figure 9a, the thermal gradients are confined to this advective–diffusive layer, as are the elastic stresses. Therefore, according to the arguments presented in §2, the crack spacing should be proportional to this effective layer thickness and inversely proportional to v . More formally, if the length L in the Péclet number is identified with the column diameter, then all columns should have the same Pe, which must be of order 1. As shown in figure 8f, this is the case: large columns cooled slower than small ones.
For drying starch, the moisture concentration ϕ is functionally equivalent to the temperature, in generating stress. There is a bottleneck in water transport as the two mechanisms of vapour diffusion and capillary flow cross over in efficiency at some critical moisture concentration . Most of the stress gradients in a drying starch body are confined to a thin layer, near this bottleneck, where an advection–diffusion problem can be posed for the moisture concentration, 6.2This is demonstrated in figure 9b,c. Here, Dh is the minimum in the moisture diffusivity, and Peh=vL/Dh [53,25]. As shown in figure 8g, the size of columnar joints in starch is inversely proportional to their cracking front velocity v, and when scaled by Dh are found at the same Péclet number as those in lava .
In summary, the spatially extended forms of columnar joints develop in a very similar way to the other hexagonal crack patterns described and modelled here, where the cycles of cracking and healing have become a stacked set of crack advances. In each advance, the network of cracks is guided by the cracks that precede it but is allowed to take a slightly different path. The sequence of cracking at any vertex is variable, however, and the vertices themselves interact and move over time. The crack pattern rapidly evolves from a rectilinear, T-junction-dominated pattern to a hexagonal one but retains a slight asymmetry in the immediate vicinity of vertices that is indicative of the order of crack opening in any particular growth increment. Finally, as with thin-film crack patterns, the length scale of the polygonal columns scales with the effective layer thickness.
7. Eroding crack patterns
Recently, a crack pattern that forms on cemented gypsum sands, on the windward (and therefore eroding) faces of dunes, has been described [6,67] and is briefly reviewed here. These crack patterns, shown in figure 1c, are found in the dunes of the White Sands National Monument, in New Mexico. Genetically, they have been argued to be similar to columnar joints, but where the cracks are advancing into the dunes as the cemented sand surface is removed . Crack depths of 0.1–45 cm were reported, and polygon sizes ranged from 5 to greater than 70 cm. As with naturally dried desiccation cracks in starch [48,24] and thermal contraction cracks in coke , the polygon size was noted to increase with depth.
It was found that cracks changed directions abruptly as vertices were approached . In this, they appear similar to the mud-crack experiments reported here and in  and to the developing polygonal terrain at New Harbour . The angles between cracks at the junctions were measured at four sites and found to deviate significantly from those of orthogonal T-junctions . Only in very thin (1 mm) cements, or cases where the cracking directions were strongly influenced by the dune slope, were T-junctions typical.
A convincing model for the origin of these cracks was also proposed . The cycling of moisture in the dunes allows for gypsum to dissolve and re-precipitate near the dune surfaces, cementing the sand. The same water, as it evaporates, causes stresses which fragment the cemented layer. Finally, the wind erodes the stoss face of the dune, allowing the wetting (and cracking) layer to slowly intrude into the dune.
Several features of these patterns, such as the local changes in crack orientation near vertices and the opportunity for dynamics as the face erodes, are consistent with an evolution similar to that outlined here. However, direct evidence for this is limited. At one site, 13 junctions were excavated and measured at the surface 2 and 10 cm depth. Although it was stated that cracks became straighter and more ordered with depth, the distribution of angles from this sample did not significantly change . As it is an interesting pattern, potentially related to ordered cracks and water cycling on Mars , additional work is to be encouraged.
The crack patterns discussed here, such as those that are commonly seen in dried mud, frozen soils or old lava flows, can take the appearance of a rectilinear network, with cracks intersecting preferentially at near right-angled T-junctions or a hexagonal network of equi-angled Y-junctions or sometimes an intermediate state. A genetic model was sketched which could account for these observations, among others, as the result of simple interactions. A crack curves to maximize the difference between the strain energy release rate, at the moment of cracking, and the cost of creating the new crack surfaces—an established principle of fracture mechanics . This is a local energy consideration: the crack tip only senses its immediate environment and responds accordingly. When a layer repeatedly cracks and heals, the positions of the previous cracks may be weaker, or more compliant, than the surrounding material. Similarly, in a crack network like that of columnar joints, which is growing through space, the stress state in an uncracked layer is strongly influenced by the prior pattern. In this way, cracks are guided along the paths taken in previous cracking cycles. At vertices, however, this argument introduces asymmetries and a growing crack can be deflected. If the order in which cracks arrive at a vertex changes with time, in different cracking cycles, these deflections can accumulate and will distort the vertex, moving it in a direction to equalize the angles between all cracks, creating a Y-junction. The requirements and predictions of this process were outlined in §§2 and 3.
This model was compared with and elaborated by presentations of evolving crack patterns in desiccated bentonite clays, polygonal terrain in permafrost, columnar joints in lava and starch, and a pattern of desiccation cracks that evolves as it erodes. These patterns were found to share similar features, such as a scaling of the crack spacing with the layer thickness, changes in the order of crack opening and distortions of the crack pattern in the immediate vicinity of vertices. Furthermore, it was shown in bentonite clay layers that were repeatedly dried and re-wetted that the drift of vertices away from their original positions, during this evolution, scales with the thickness of the cracking layer. This rate may also depend on the materials involved—as permafrost apparently takes much longer to order than a few tens of cracking cycles, measured there by years. The direction of this motion was seen to be precisely that predicted to deform vertices towards a hexagonal crack pattern.
Note added in proof
After acceptance of this article, I was made aware of a recent publication  which demonstrates that the head scales of crocodiles form by crack growth and evolution during embryonic development. This finding is connected to the ideas developed here, and the report is to be recommended.
I thank C. Duquette for her assistance in the field studying columnar joints, and R. Sletten, B. Hagedorn and D. Mann for providing the opportunity to participate in fieldwork in Antarctica and visit the Black and Berg sites. I also particularly thank B. Hallet for discussions on this topic over many years and P. Nandakishore for suggesting an improvement to the algorithms to measure crack angles.
One contribution of 13 to a Theme Issue ‘Pattern formation in the geosciences’.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.