Nonequilibrium scale selection mechanism for columnar jointing
See allHide authors and affiliations

Edited by H. Jay Melosh, University of Arizona, Tucson, AZ, and approved October 31, 2008 (received for review May 29, 2008)
Abstract
Crack patterns in laboratory experiments on thick samples of drying cornstarch are geometrically similar to columnar joints in cooling lava found at geological sites such as the Giant's Causeway. We present measurements of the crack spacing from both laboratory and geological investigations of columnar jointing, and show how these data can be collapsed onto a single master scaling curve. This is due to the underlying mathematical similarity between theories for the cracking of solids induced by differential drying or by cooling. We use this theory to give a simple quantitative explanation of how these geometrically similar crack patterns arise from a single dynamical law rooted in the nonequilibrium nature of the phenomena. We also give scaling relations for the characteristic crack spacing in other limits consistent with our experiments and observations, and discuss the implications of our results for the control of crack patterns in thin and thick solid films.
Drying solids lose moisture from their exposed surfaces and shrink as a consequence. Similarly, cooling solids lose heat from their exposed surfaces and shrink as a consequence. In either case, this differential shrinkage of one part of the solid relative to another leads to stresses that can eventually lead to cracking (1–3). Although much is known about the nucleation, growth, dynamics, and stability of a single crack in an elastic solid, most questions associated with the patterns of multiple cracks due to stresses that arise from nonequilibrium processes such as drying and cooling (4–7) remain wide open. The resulting polygonal planform patterns can arise in a variety of situations, from the mundane cracks in drying mud, to the deliberately artistic cracks in ceramics and pottery, to the spectacular columnar joint formations of the Giant's Causeway in Northern Ireland, Fingal's Cave on Staffa, in Scotland, and the Devil's Postpile in California. The latter formations have fascinated casual observers, artists, and scientists for centuries (8–10), but no comprehensive physical theory for their form or scale exists. Indeed, it is only in the past decade or so that careful laboratory experiments have started to address the dependence of any of these crack patterns on such quantities as the rate of drying or cooling, the thickness of the layers, and their mechanical properties (4–7, 11–13). For example, recent experiments show that crack formation and propagation in drying thin films leads to length scales and patterns that can be strongly timedependent; cracks in directionally drying films grow diffusively at short times, and can advance intermittently via sticksliplike motion over longer times (11, 12). The patterns formed by these cracks depend in detail on the spatiotemporal dynamics of drying, substrate adhesion, and thickness variations (4, 6, 13, 14). This immediately suggests a nonequilibrium origin to these crack patterns, one that couples the heterogeneous elastic stresses in the cracking solid to the dynamics of drying or cooling, that might be contrasted with the equilibrium crack patterns that are seen and studied in a variety of engineering applications (15).
For columnar joints like those shown in Fig. 1, which only occur in relatively thick layers, the similarity between crack patterns induced by drying (16–21) and cooling (16, 22–25) can be traced to the fact that the transport of water in a drying slurry and the extraction of heat from a hot solid are mathematically analogous (1–3). In a poroelastic medium, fluid flow is coupled to the elastic deformation of a porous solid, whereas in a thermoelastic medium, heat conduction is coupled to elastic deformation of a conducting solid; pressure in one case is analogous to temperature in the other. In each case, a shrinkage front propagates through the medium leading to the penetration of a crack front that follows slightly behind. Ordering of the crack network at this front carves out the regular columns. In this article, we will show how the scale of this shrinkage front sets the average size of the resulting columnar joints. To understand this quantitatively, we present observations on both laboratory and geological columnar joints, and then deduce the appropriate scaling behavior from general theoretical considerations. The resulting scaling law is also applicable to a wide class of fracture problems where the elastic screening length of the crack pattern is smaller than the sample thickness.
Experimental Observations on Starch
We studied laboratory examples of columnar jointing made by drying slurries of corn starch (see Fig. 1A and D), that have been known since at least Victorian times (9), although they have only recently been investigated quantitatively (17–21).
Slurries of corn starch were prepared by mixing equal weights of dry (Canada brand) corn starch and water. Traces of bleach were added to sterilize the experiments, and samples were dried in glass dishes under a pair of 250 W heat lamps. The sample mass was automatically recorded every minute, and these data were used in a feedback loop to control the evaporation rate. If the observed evaporation rate was lower than desired, the duty cycle of the heat lamps was proportionally increased, to deliver more heat to the starchcake, and vice versa. Using these methods, which are described in more detail in ref. 21, we could fix the evaporation rate to within 10% of a desired value throughout the experiment. As discussed below, the evaporation rate controls the rate of fracture advance and the scale of the columns that are formed. In runs with feedback control, the evaporation rate was held constant, and a stable average column size was reached within 1 cm of the drying surface, and maintained to the base of the starchcake (21). In runs without feedback control, the evaporation rate decreases decreased with time, the fracture front slowed, and the columns became larger with depth (19, 21, 25). We measured column scales in both controlled and uncontrolled runs and were thus able to study the column scale over over a wide range of fracture advance rates.
The fracture position as a function of time can be deduced from the sample mass as a function of time by using a data inversion technique that integrates the water concentration data, presented in Fig. 2A and B, to construct a lookup table for the fracture position as a function of sample mass. Although the exact shape of the water concentration front depends on the evaporation rate (both current and past), these details do not substantially alter the relationship of the total sample mass to fracture position, which reflects the almost steplike jump in water concentration at the drying front, evident in Fig. 2A and B. To test this technique, we dried a series of 17 samples under heat lamps, using identical drying conditions for each sample. These were broken open after different drying periods had elapsed, and the fracture position was measured directly in each starchcake. The sample mass, which had been measured every minute, was then converted into an estimated fracture position by using the inversion algorithm. As shown in Fig. 3A, there is excellent agreement between direct observations of the fracture front position, and the position inferred from continuously weighing the starchcake.
The fracture front velocity v is taken to be the numerical derivative of the fracture position with respect to time. As shown in Fig. 3B, we can use this technique to measure the velocity of the fracture front as it passed through the sample, after drying is complete. Typically, in our controlled experiments, a desired average fracture velocity was well achieved in the middle 6075% of the sample. The velocity v is simply related to the volumetric flux of evaporated water per unit area J_{0} and the sample porosity ϕ by v = J_{0}/ϕ. All experiments were consistent with a porosity of ϕ = 26 ± 1%.
Feedback control produced very regular starch joints within a certain narrow size range, by using drying fronts that moved through the sample at fixed v. An example is shown in Fig. 1D. To probe a larger range of v, we let thick starchcakes directionally dry without feedback control, resulting in columns that slowly increase in scale. In both cases, we measured how the fracture spacing, defined as the square root of the average column crosssectional area, depends on v. Fig. 4A shows that the fracture spacing is inversely proportional to the instantaneous fracture front velocity v, so that faster drying leads to smaller columns.
Field Observations on Lava
It is not possible to make direct, dynamical observations on columnar joints as they slowly form in cooling lava flows, although some important insights can be had from borehole temperature measurements that have been made over the decades following the formation of the Kilauea Iki lava lake (26). It is, however, possible to deduce the cooling rate, fracture advance velocity, and column scale from purely geometric characteristics of joints in ancient, exposed flows (25).
We measured columnar joints at field sites across the Columbia Plateau in Washington State, in southwestern British Columbia, and on the island of Staffa. The Columbia River Basalt Group covers central Washington state with up to a kilometer of columnar basalt, and extends into Oregon and Idaho (27). This Miocene flood basalt plateau contains several hundred extensive, chemically homogeneous flows (27, 28). We used this terrain as a natural laboratory with repeatable lava properties and deposition histories. Nearby, the volcanism of British Columbia allows, in contrast, data to be gathered from highly varied eruption and lava types. Sites in the Garibaldi volcanic belt and the Intermontane belt allowed us to study jointing in this more heterogeneous environment, including the diverse small recent flows of columnar dacite, trachydacite, and alkali basalt in the Vancouver–Whistler corridor (29–31). Finally, the lava on the island of Staffa in Scotland is part of the North Atlantic Tertiary Igneous Province, and is compositionally similar to that of the Columbia Plateau, but slightly more basic (32). We studied the Fingals cave lava unit at several places along the southern coast of Staffa.
At each field site, the average column face width and stria height were measured. Striae are chisellike marks (see Fig. 1E), which are left by the individual, sequential fracture advances of a growing crack as it intrudes into the slowly cooling lava (23, 24, 33). These features, which are also called joint increments (23), are related to the commonly seen striations characteristic of metal fatigue (33), and always form perpendicular to the direction of crack motion (33). More than a few meters away from any cooling surface, the reflux of water within cracks is the dominant cooling mechanism of lava (26), and establishes a selfsustaining cooling front moving at a constant speed, v.
The temperature gradients near the moving cooling front arise from a competition between the diffusion of heat and the advance of the front at speed v, and it can be shown that they scale inversely with v (25). The striae leave a permanent record of the temperature field. As each stria forms, the crack tip advances, on average, over a drop in thermal stress given by the tensile strength of the host rock (25, 33). This implies, in turn, that v and the average stria height are inversely proportional to each other. The coefficient of proportionality linking these terms can be calculated theoretically by using the more detailed methods presented in ref. 25. Here, we directly apply this theory, assuming that the physical properties of lavas in all locations were the same as those of the tholeiitic basalts of the Columbia Plateau. Although this is a valid assumption for the Staffa lava, which is compositionally similar to those of the Columbia Plateau, it may contribute to systematic uncertainties in the results from British Columbia.
In general, we found that large columns are due to slower cooling than small columns, consistent with qualitative inferences made previously by using crystal texture analysis (22). Fig. 4B shows that the fracture spacing of columnar joints in lava is inversely proportional to the speed of the cooling front over a wide range of speeds, and is consistent with the linear relationship between stria height and column face width (23, 25, 34). We see that the same general relationship between fracture spacing and front speed holds for both cooling lava and drying corn starch slurries. To understand this, we now turn to the mathematical similarity between these physical processes.
Scaling Theory
Our dynamic similarity argument rests on the observation that the processes of drying and cooling of a deformable solid are mathematically analogous. In either case, we consider an elastic solid for which the relative displacement of material is described by a strain tensor
When appropriate boundary conditions are specified for T or p, displacements and/or tractions, the force and energy balance equations completely determine the state of the solid. We see that for a poroelastic solid, pressure is mathematically analogous to temperature, up to factors of the coefficient of thermal expansion and an elastic modulus; indeed the underlying similarity has a thermodynamic origin (1). Given this, we see that the dynamic similarity of columnar joints in starch and lava arises because the transport of the relevant quantity, water or heat, is essentially diffusive in each case. Indeed, for a linearly elastic isotropic solid the Eq. 1, yield, in the limit of slow deformations, the uncoupled thermal diffusion equation ∂T/∂t = ∇ · (D_{T}(∇T)), with the thermal diffusivity D_{T} = k/ρC_{p}. For poroelasticity, we get an analogous equation of the form ∂(∇ · u + βp)/∂t = ∇ · (D_{P}(∇ · u + β∇p)), with the hydraulic or poroelastic diffusivity D_{P} = k_{h}/H. We see that the diffusion of heat in thermoelasticity is analogous to the diffusion of water mass content in poroelasticity. To complete the formulation of the problem we need to specify boundary conditions of zero normal traction and a prescribed flux of either heat or moisture on the exposed surface of the halfspace. Here, we do not solve this problem, but restrict ourselves to a simple testable scaling theory, starting with the case of drying starch, and then transposing our results to the geological problem.
Following an initial transient after a drying/cooling experiment is started, a propagating nonequilibrium steady state of drying/cooling is set up in the solid that leads to the formation of a fracture front that relieves stress in the solid behind it. During the transient drying process near the exposed free surface of a wet, thick starchcake of large lateral extent, the surface is exposed to low humidity air. After a time t following the beginning of the experiment, a surface layer of thickness L_{s} shrinks because of the loss of water from it. Shrinkage is easily accommodated in the direction normal to the surface, but leads to an isotropic tensile stress in the plane parallel to it. If the stresses induced by drying are sufficient to exceed some cracking threshold, then cracks will nucleate from the free surface and propagate into the bulk material. Each crack will relieve stress in a region given by the elastic screening or stress relaxation length perpendicular to the crack face, which is proportional to L_{s} (15). However, because the drying/cooling field inducing stresses is dynamic, the cracks move from the highly stressed drier surface into the wetter and less stressed bulk, and so slow down and eventually stop. Further drying leads to a repetition of this cycle. Although individual cracks move intermittently via a sticksliplike motion, on average the crack front advances diffusively, following the drying front, albeit slightly behind it (12).
During the initial phase of drying near the surface of a very thick sample, approximated as a halfspace, the only dynamical length scale that arises naturally is given by where D_{P} is the poroelastic diffusivity characterizing the equilibration of fluid mass content. If the evaporative flux J_{0} decreases with time, as when the external humidity is kept constant, the crack front slows as the transport of water is limited by diffusion out the sample. As time progresses, L_{s} increases and the crack pattern should coarsen, i.e., the crack spacing should increase with increasing distance from the exposed surface. This coarsening has been previously noticed in starch experiments (17, 19, 20). Our explanation for the coarsening follows from the nonequilibrium nature of drying/cooling as outlined above, in sharp contrast with equilibrium theories given hitherto for these crack patterns that are based on freeenergy minimization (36). When the evaporation rate is kept constant, so that the volumetric flux J_{0} is invariant in time, we find ourselves in a nonequilibrium steady state. This is achieved in the fluxcontrolled starch experiments presented earlier, where, after a transient near the surface, the crack spacing does not coarsen further. Achieving this steady state requires that the drying power increase with time. In this case, the drying front progresses into the sample at a constant speed v. Under these conditions, the elastic stress screening length is also invariant, and scales as This allows us to define a dimensionless Péclet number for a crack spacing L_{c}, where the drying front velocity v is related to the moisture flux J_{0} and the porosity ϕ by v = J_{0}/ϕ. Since Pe ∼ L_{c}/L_{s}, the observed inverse scaling of fracture spacing to fracture advance rate follows if jointing proceeds at a constant Pe ∼ O(1).
Similar considerations apply to the lava cooling case, with the poroelastic diffusivity D_{P} replaced by the thermal diffusivity D_{T} = k/ρC_{p}. In the case of cooling lava, coarsening has been observed near the flow margin (23), whereas the nonequilibrium steady state is maintained far away from the flow margin by the heat transport due to the reflux of water in the cracks (24–26). In this regime, a very regular column scale is observed as in Fig. 1C and D.
Dynamic Similarity of Starch and Lava Columns
Our scaling theory simplifies somewhat for the case of drying starch, even though the relationship between the local water concentration C and the pressure p is complicated by the fact that the pores are not fully saturated with water, but instead consist of water bridges, air, and water vapor. During the drying process, water is transported under the pressure gradient predominantly as liquid at the wet end of the sample, and predominantly as vapor near the dry end (20, 37). Because linear poroelasticity theory shows that the water mass content satisfies a linear diffusion equation, we expect that the moisture concentration C, which we measured directly in the experiments, also obeys a diffusion equation. However, given the large variations in water content that lead to fundamentally different mechanisms for moisture transport at early and late times (37), the diffusivity is a function of the moisture content in the solid so that C now satisfies a nonlinear diffusion equation of the form where D(C) is a concentrationdependent effective diffusivity (related to the poroelastic diffusivity D_{P} ∼ k_{h}/H) which can be deduced from the experimental data.
In our starch experiments, we controlled the volumetric flux of water J_{0} at the surface, rather than the vapor pressure p. We sampled the moisture concentration C(z,t) in the vertical direction z in many identical starchcakes dried for different times t, by sectioning the cakes into 1 to 2mmthick layers, and weighing, drying, and reweighing the sections. As shown in Fig. 2A, the cakes dried uniformly to a critical moisture concentration of 0.30 g/cm^{3}, after which time a sharp drying front was initiated at the upper drying surface. Below this front, C remained constant. As the drying front steadily advanced, it maintained a concave, selfsimilar shape, as shown in Fig. 2B. To understand this, we consider Eq. 6 that can be integrated to give the concentrationdependent diffusivity (37), Here, we have assumed 1dimensional transport, and that the sample lies in 0 ≤ z ≤ h. We have imposed a noflux boundary condition on the lower boundary z = h, which is in contact with the base of the container. Fig. 2C shows the effective diffusivity D(C) as calculated from Eq. 7 and the water concentration data presented in Fig. 2A. There is a broad minimum in D(C), reaching D_{0} = 1.1 × 10^{−9} m^{2}/s over approximately C = 0.1 − 0.3 g/cm^{3}. Above a concentration of 0.3 g/cm^{3} the diffusivity must continue to increase with C, as the starchcake dries almost uniformly from saturation down to this concentration. However, under these wet conditions, ∂C/∂x is very small and these data cannot be reliably used to find D(C). The minimum in D(C), which we call D_{0}, is due to the crossover between liquid and vapor transport (37) that occurs at the sharp dropoff in water concentration which forms the drying front. It is therefore this minimum D_{0} which is relevant to the scaling of the columns. By using Eq. 5 with D_{P} = D_{0}, the measured starch fracture spacing shown in Fig. 4A is thus found to correspond to a range of Péclet numbers Pe = 0.15 ± 0.05.
In the case of lava, the energy balance Eq. 1, reduces to a linear diffusion equation with a diffusivity D_{T} = k/ρC_{p}. Then, we transpose the result for drying starch to the case of cooling lava by reinterpreting D as the thermal diffusivity D_{T} and J_{0} as the heat flux that is strongly influenced by the continuous reflux of water in the cracks (24, 25), although this effect does not significantly influence the water transport in the case of drying starch (21). Because the thermal diffusivity of basalt (25) is approximately D_{T} = 6.5 × 10^{−7} m^{2}/s, and the fracture front velocities are comparable to those in the starch experiments, columnar jointing in lava should be 10100 times larger than in starch. We see that this is indeed the case, as shown in Fig. 4B; jointing in lava is well described by Pe = 0.3 ± 0.2, very similar to that found in the starch experiments.
Discussion
Our simple scaling theory embodied in Eqs. 3–5 provides a unified way of understanding the basic observations of columnar jointing induced by the nonequilibrium and inhomogeneous cooling in a geological context, and via a mathematical analogy, the columnar jointing in laboratory experiments of the drying of a fluid infiltrated porous solid such as starch. In particular, we are able to derive a simple scaling law for the crack spacing consistent with both previous and current experimental observations.
In addition to explaining an old geological mystery, and connecting it to a number of everyday observations of cracking, our work paves the way for the engineered cracking of solids to obtain patterns at will. To take one example, in the production of fine coke from crushed coal, it is the coolinginduced shrinkage cracks that control the size, and hence quality, of the product (38). Other applications range from the industrial drying of paint, concrete, and photographic paper, to mudcracks and the desiccation of soils, the thermal fracture of permafrost (39), and the jointing of metamorphic rocks (S. ShorelandBall, R. Sparks, M. Murphy, C. MacNicaill, and B. Daniel, unpublished work). In all of these examples, the crack patterns are controlled by the competition between nonequilibrium processes such as the relatively slow drying and cooling by which stresses develop, and rapid cracking that causes these stresses to be relieved. Thus, although elastic deformations are important, the dynamics of crack propagation are much faster than the very slow process of stress generation that ultimately determines the rate of cracking and thence the scale of jointing.
We conclude with a brief discussion of the drying or cooling of a solid of finite thickness h. For very thin or very slowly cooling (or drying) solids, the crack spacing will be ∼ O(h), if cracks form at all (4, 5, 13). However, if the characteristic scale given by Eq. 3 is smaller than h, the fracture spacing ∼ L_{s}, no longer depends on h, but is instead a function of the cooling (or drying) rate. More generally, if the characteristic elastic stress screening length is smaller than the thickness of the drying solid, then the crack spacing is just the stress relaxation length, so that the crack spacing should be proportional to min(h,L_{s}).
Acknowledgments
We thank Catherine Duquette, Robert Mehew, and Brian Goehring for their assistance in the field, National Trust for Scotland for allowing extended access to Staffa, and PierreYves Robin, Mark Jellinek, and Kelly Russell for their advice on geology.
Footnotes
 ^{2}To whom correspondence should be addressed. Email: lg352{at}cam.ac.uk, lm{at}seas.harvard.edu, or smorris{at}physics.utoronto.ca

Author contributions: L.G., L.M. and S.W.M. designed research; L.G. performed research; L.M. contributed analytic tools; L.G. analyzed data; and L.G., L.M., and S.W.M. wrote the paper.

↵^{1}Present address: BP Institute for Multiphase Flow, University of Cambridge, Madingley Road, Cambridge, UK CB3OEZ.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Freely available online through the PNAS open access option.
 © 2009 by The National Academy of Sciences of the USA
References
 ↵
 ↵
 ↵
 Boley BA,
 Weiner JH
 ↵
 ↵
 ↵
 ↵
 ↵
 Bulkeley R
 ↵
 Huxley TH
 ↵
 Kennedy A
 ↵
 ↵
 ↵
 ↵
 ↵
 Hutchinson JW,
 Suo Z
 ↵
 French JW
 ↵
 Müller G
 ↵
 Toramaru A,
 Matsumoto T
 ↵
 ↵
 ↵
 ↵
 Long PE,
 Wood BJ
 ↵
 ↵
 ↵
 Goehring L,
 Morris SW
 ↵
 ↵
 Hooper PR,
 Kawkesworth CJ
 ↵
 Reidel SP,
 et al.
 ↵
 Green NL
 ↵
 Bye A,
 Edwards BR,
 Hickson CJ
 ↵
 Stelling P,
 Tucker DS
 Russell JK,
 Hickson CJ,
 Andrews G
 ↵
 Kerr AC
 ↵
 Ryan MP,
 Sammis CG
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
Citation Manager Formats
Article Classifications
 Physical Sciences
 Physics