## Abstract

Progress towards the development of a better understanding of the formation of geological patterns in wet systems due to precipitation and dissolution is reviewed. Emphasis is placed on the formation of terraces, stalactites, stalagmites and other carbonate patterns due to precipitation from flowing supersaturated solutions and the formation of scallops by dissolution in undersaturated turbulent fluids. In addition, the formation of spherulites, dendrites and very large, essentially euhedral, crystals is discussed. In most cases, the formation of very similar patterns as a result of the freezing/melting of ice and the precipitation/dissolution of minerals strongly suggests that complexity associated with aqueous chemistry, interfacial chemistry and biological processes has only a secondary effect on these pattern formation processes.

## 1. Introduction

Progress towards understanding and simulating some striking and often beautiful geological patterns formed by mineral precipitation and dissolution is reviewed. At present, the origins of the vast majority of macroscopic geological patterns are not well understood, and much work remains to be done. We start by discussing systems in which pattern formation is dominated by the precipitation of dissolved minerals from flowing fluids at or near the Earth’s surface, with a focus on the growth of terraces and steps, and precipitation from slowly flowing thin liquid films to form stalactites, stalagmites, icicles, travertine domes and other formations. Terraces and steps are typically formed where there is a plentiful flow of water and where at least part of the flow is turbulent. This adds substantial complexity to attempts to accurately simulate these patterns. When the supply of water is less plentiful, dissolved minerals and other substances are transported in thin films under laminar flow conditions, and the hydrodynamics is much simpler (Oron *et al.* 1997).

Pattern formation in aqueous systems occurs under a wide variety of biogeochemical conditions, but the most familiar and extensively investigated patterns are brought about by the precipitation of carbonate minerals controlled primarily by the loss of CO_{2} from the aqueous phase due to transport to the atmosphere. In cave waters, the primary source of CO_{2} is usually the respiration of micro-organisms, plants and animals, in soil (Raich & Schlessinger 1992), which raises the CO_{2} chemical potential far above that in typical cave atmospheres. Water carrying dissolved CO_{2} percolates through the ground and dissolves carbonate minerals. When water with a calcium concentration that has been elevated in this manner reaches a cave, the process is reversed, and carbonate minerals are precipitated, as dissolved CO_{2} is transported from the aqueous phase to the vapour phase resulting in a pH increase and associated calcium carbonate supersaturation. In hot springs, the source of CO_{2} is mainly abiotic (Crossey *et al*. 2009), and the CO_{2} potential (concentration) is often much larger than biological processes can produce. Consequently, the rate of precipitation is generally higher near hot springs than that in caves.

Many pattern formation processes that are dominated by mineral precipitation and/or dissolution are paralleled by pattern formation by the freezing and/or melting of ice (e.g. Short *et al.*2005, 2006). This provides important model systems that can be conveniently used in laboratory experiments. Ice has the advantage that freezing and melting is very fast relative to precipitation and dissolution in natural systems, there are no geochemical complications and micro-organisms usually play no role in the pattern formation process (Lock 1990). The observation that patterns formed by freezing and melting are strikingly similar to those formed by precipitation and dissolution strongly suggests that pattern formation can often be explained without including all of the inherent biogeochemical and geological complexity (Goldenfeld *et al.* 2006; Hammer *et al.* 2007).

The next topic that we review is the formation of patterns consisting of a dense array of scallop-shaped depressions. These patterns are formed by the coupling between fluid flow and material removal processes such as dissolution, melting, corrasion, corrosion and ablation. The origin of these patterns is understood only in qualitative terms. The most familiar example of this class of patterns is scallops in the walls of phreatic cave passages, and the size of the scallops is often used to estimate the flow conditions at the time of their formation (Blumberg & Curl 1974).

The next two sections are concerned with extremely different patterns that are formed by the growth of solids from liquids and melts—spherulites and dendrites. Spherulites, as the name indicates, are essentially spherical forms, whereas dendrites are branched. Highly branched dendrites are formed from solution under diffusion-limited quasi-stationary conditions when the far field solute concentration is low. The solute concentration field can be described, to a good approximation, by the Laplace equation, with appropriate boundary conditions (Langer 1980; Meakin 1998). The growth of spherulites occurs when the solute concentration is high, and it is usually controlled by the diffusion of impurities that are rejected by the growing solid phase and form an impurity boundary layer (Goldenfeld 1987; Granasy *et al*. 2005). Spherulites are often formed by materials that are highly anisotropic, but this anisotropy is not expressed by their essentially spherical shape, because of their polycrystalline internal structure (Goldenfeld 1987). Despite their highly branched shapes, dendrites are often single crystals, and their internal crystal anisotropies are expressed by their almost symmetric forms.

Small euhedral crystals are almost ubiquitous in geological systems, but large crystals are much less common. The last topic reviewed is the growth of spectacular crystals with dimensions of the order of 10 m.

## 2. Terraces and steps

A sequence of steps or terraces is often formed in fluvial systems (figures 1 and 2), and in these systems each terrace consists of a barrier, dam or ridge that impounds a pool of water. In qualitative terms, it appears that these patterns are formed by the trapping of particles and/or preferential precipitation at the terrace edges where the water is shallowest and flows most rapidly (the terraces are formed by the amplification of perturbations due to the decreased fluid depth and increased fluid velocity associated with them; Buhmann & Dreybrodt 1985; Hammer *et al*. 2008).

In many fluvial systems, the flow is very variable, and when the morphological dynamics is controlled by physical processes such as sediment erosion, transport and deposition, the evolution of the system may occur primarily under high-flow conditions (Burkham 1972). For example, stream step pools are formed when the discharge is high, and vegetation arcs and litter dams may be dry for most of the year (Eddy *et al*. 1999). Abrahams *et al.* (1995) proposed that step pools (figure 1*a*) form where the maximum flow is just sufficient to move the largest debris, and that the step pool pattern evolves to the most stable form in which the flow resistance is maximal. This concept predicts that *h*/*l**s*, where *l* is the distance between steps, *h* is the step height and *s* is the stream channel slope, is slightly greater than unity, which appears to be consistent with field data. In other cases, *h*/*l**s* may be near to unity because steps stop growing if they become immersed in the pool behind a downstream dam. A similar phenomenon occurs in the formation of travertine dams (see below), which stop growing if they become inundated by the pool of water that accumulates behind a much larger downstream dam. Despite substantial effort, the formation and stability of stream step pools is not well understood (Church & Zimmermann 2007). There is no compelling reason to believe that the similarities between step pools and other stepped structures in fluvial systems imply that similar mechanisms are at work. For example, the formation of step pools is a physical process, whereas growth of a series of travertine dams in a stream bed may form a similar morphology (Pentecost 2005) due to chemical (precipitation) and/or biologically mediated precipitation and particle capture.

Figure 1*b*–*d* shows similar patterns formed in fluvial systems. The terraces shown in figure 1*b*,*d* are formed by the precipitation of minerals, and the series of dams shown in figure 1*c* is formed from pine needles that cemented together by the precipitation of silicaceous minerals (it appears that they grow by capturing pine needles and other floating debris at their tops). In all these systems, micro-organisms are abundant, and it appears that they play a role in most precipitation processes at and/or near the Earth’s surface. However, it is not clear whether micro-organisms have a controlling influence on the terrace morphology. Blank *et al.* (2002) found that the composition of the micro-organism community had little influence on the geyserite (hydrated silica) morphology at hot springs in Yellowstone, and the formation of terraces can be explained on the basis of abiotic mechanism alone. Macrobiological processes clearly play an important role in the formation of some terrace patterns such as banded vegetation patterns on sloping ground. Several models have been developed to simulate banded vegetation patterns including cellular automata models (Mauchamp *et al.* 1994; Thiery *et al.* 1995) and a continuum model (Lefever & Lejeune 1997), which represents the vegetation density by a continuous field. These models focus on the ecological and biological factors (competition, symbiosis and reproduction) that control plant population dynamics, although hydrological factors are included to varying degrees.

The rate at which mineral terraces on the Earth’s surface grow is controlled by the rate at which dissolved mineral precursors are supplied by the source and by the rate at which the aqueous chemistry and/or temperature changes to bring about supersaturated conditions. Very often, the growth rate is very fast. For example, figure 1*d* shows a set of large terraces formed by the oxidation of iron in highly acidic water (acid mine leachate) which initially contained about 2 g l^{−1} Fe(II) under essentially anoxic conditions. Absorption of oxygen from the atmosphere and micro-organism photosynthesis resulted in a rapid increase in the oxygen content of the water, and a combination of organic (primarily microbial) and inorganic processes oxidized Fe(II) to Fe(III) leading to the precipitation of a variety of hydrated iron oxyhydroxides and hydroxy sulphates, with a growth rate of the order of 1 cm yr^{−1}.

The growth of travertine terraces, which form when water with a high concentration of dissolved calcium emerges from the ground, may be even faster (rates up to 0.5 cm d^{−1} (Goldenfeld *et al.* 2006) and up to 30 cm yr^{−1} (Fouke *et al*. 2000) have been reported), and this provides a particularly spectacular example of geological pattern formation (figure 2). By comparison, the growth rate of stalagmites is typically 10^{−3}–10^{−2} cm yr^{−1} in cool-region caves and 3–10×10^{−2} cm yr^{−1} in subtropical regions (Fairchild *et al*. 2006). In most cases, the high concentration of dissolved calcium is a consequence of the dissolution of carbonate minerals in water with a very large CO_{2} content, and the precipitation of carbonate minerals (travertine) is controlled primarily by the loss of carbon dioxide. The overall precipitation reaction is
2.1
One molecule of CO_{2} is formed for each molecule of CaCO_{3} precipitated, and, if this CO_{2} was not released, precipitation would stop because of the increased H_{2}CO_{3} activity.

The solubility of amorphous silica in water increases with increasing temperature and depends on the pH value (Alexander *et al.* 1954). The formation of geyserite steps due to the precipitation of silicaceous material could be due to cooling, evaporation, pH change or a combination of these factors. The rapid formation of ice terraces under flowing water suggests that, although micro-organisms are abundant at mineral–fluid interfaces under surface water, they are not essential to the formation of terraces. The solubility of silica is quite low, even at high pH, and the growth of silica terraces is much slower than that of travertine (typically 0.1 cm yr^{−1} or less; Usurov *et al*. (2007)).

A variety of terrace patterns can be observed when water in contact with cold air flows over ice in rivers and creeks. In this case, water loses heat to the cold air, and freezing occurs when water at the ice–water interface reaches a temperature of ≈0^{°}C due to thermal conduction and advection of undercooled water. This process can be extremely rapid because of the essentially unlimited supply of solid precursor (water) and the high rate of heat transport relative to molecular diffusion (heat and chemical transport due to eddy diffusion is rapid in the turbulent fluids, and the freezing process is controlled by heat diffusion across the laminar sublayers near the air–water and/or water–ice interfaces). Similarly, Hardcastle (1920) has reported the formation of a large ice dam near a glacial spring and attributed it to supercooled water issuing from the spring.

A number of investigators have noticed that growth must occur preferentially over the lips of the terraces, where the water depth is small, and the fluid flow is fast. However, several underlying mechanisms have been proposed to explain this correlation. These include: (i) the thickness of the laminar sublayer in the flowing fluid is reduced leading to more rapid diffusive transport to the interface (Buhmann & Dreybrodt 1985); (ii) increased CO_{2} degassing due to a larger water–air interface area (relative to the water volume; Chen *et al*. (2004)); (iii) increased capture of particles (Chafetz & Folk 1985), possibly mediated by biofilm (Hammer *et al.* 2007); and (iv) the water stretches as it flows over the terrace lip, and this increases the concentration gradients at the solid–liquid interface, thus leading to more rapid diffusion towards or away from the interface down larger concentration gradients (Hammer *et al*. 2008). The idea that the growth rate increases with increasing flow velocity, rather than with decreasing water depth, is supported by the observation that height perturbations not only increase the local flow velocity but also the water depth decays (a notch in the terrace edge heals). However, Hammer *et al*. (2008) conclude that the locus of maximum growth (where the effect of concentration gradient compression is greatest) is a little upstream from the locus of maximum velocity.

The growth of mineral terraces is extremely complex. Apart from the uncertain role played by micro-organisms, the fluid flow and chemistry is complex and the constantly changing terraces themselves have a complex morphology over a wide range of scales. For example, centimetre-scale ‘mini terraces’ form on steep slopes, and ‘stalactites’ grow from the overhanging lips of large terraces. In addition, the external conditions, and often the flow rate, chemistry and temperature of the water, may change over time. Nevertheless, the formation of terraces is a robust process, and important progress has been made towards the development of theoretical (Wooding 1991) and numerical models for terrace growth (Goldenfeld *et al.* 2006; Veysey 2006; Chan & Goldenfeld 2007; Hammer *et al.*2007, 2008). Because the time scales associated with the evolution of the solid surface are very large compared with those associated with the fluid dynamics, a quasi-stationary approximation can be used (the fluid dynamics can be solved with a stationary boundary, but this must be done many times during a simulation). This type of time-scale separation is important to the success of many numerical models for geological pattern formation.

The model developed by Goldenfeld *et al.* (2006) (illustrated in figure 3*a*) represents the solid surface by a set of heights, {*h*(*i*,*j*)}, above each element of a square grid. The corresponding water depths, {*w*(*i*,*j*)}, and the air–water interface heights, {*H*(*i*,*j*)=*h*(*i*,*j*)+*w*(*i*,*j*)}, are also important quantities in the model. The flow of water out of each element to lower elements conserves the fluid volume and is proportional to *f*(∇*h*), where for ∇*h*≤∇*h*_{c} and for ∇*h*>∇*h*_{c} (to represent a transition from laminar flow to turbulent flow described by the Chezy (1776)) equation. In the model, equation (2.1) is used to represent the water chemistry, the CO_{2} outgassing is proportional to and the rate of vertical growth is , where *C*_{1} is the calcium ion concentration, *C*_{2} is the dissolved CO_{2} concentration, the term proportional to , where *c*_{t} is a constant, represents the effect of turbulence on CO_{2} outgassing, **F** is the water flux, is the unit vector normal to the solid surface, *R*_{1} and *R*_{2} are constants and . Capillary forces are important near the rim of a dam and these are represented by allowing flow only if the water depth, *w*, exceeds a small threshold value. Undersaturated water ((*C*_{1}−*C*_{2})_{0}<0) is injected at a constant rate onto an initial ‘landscape’, {*h*_{0}(*i*,*j*)}, and the water is levelled in the ponds that form behind the growing dams at all times. Based on their simulation results, Goldenfeld *et al*. (Goldenfeld *et al.* 2006; Chan & Goldenfeld 2007) concluded that the distribution of pond areas is a power law using data that spanned a range of about 100 in area (≈10 in length). This is consistent with a very limited amount of field data from a single spring. Acquisition and analysis of more extensive data are needed to provide a good test for the model.

The model of Hammer *et al.* (2007; illustrated in figure 3*b*) is based on the idea that the rate of growth of the travertine, in the direction normal to the travertine–water interface, increases with increasing flow rate, and this is described by the empirical equation
2.2
where *v*=|**v**| is the flow velocity, *f*(*v*) is the ‘precipitation scaling factor’ and *c* is a constant. In equation (2.2), the factor of (1+(∇*h*)^{2})^{1/2} converts normal growth to vertical growth, the term fills in concave depressions, where ∇^{2}*h* is positive, and suppresses convex mounds where ∇^{2}*h* is negative and the particle transport coefficient, , controls the rate of growth on small scales. In the simulations, *f*(*v*) is assumed to be a linear function of the flow velocity, *ν*, and the flow is simulated using Hydro2de, which uses a finite volume method to solve the depth averaged Navier–Stokes equation (Connell *et al*. 1998).

Both modelling approaches (those of Goldenfeld *et al*. (2006) and Hammer *et al*. (2007)) require quite severe approximations of the fluid dynamics, water and water/mineral chemistry, solute transport, etc. to make the simulations practical on a workstation. However, they both generate patterns that capture some of the most striking features observed in mineral terraces. These models are important because they provide a means of testing conceptual models that encapsulate our understanding of the formation of terrace patterns.

Figures 1 and 2 indicate that the shapes of the terraces are quite different in different systems, so it should not be surprising that the models of Goldenfeld *et al*. (2006) and Hammer *et al*. (2007) also generate terraces with different shapes. This might be due to the parameters that were used to illustrate the models, but more likely it reflects the quite different simplifying assumptions that they use. These models were developed to simulate the terrace patterns formed by hot springs. In this context, it is evident from a comparison of figures 2 and 3 that there are important differences between the patterns generated by computer simulations and those observed in nature. The primary reasons for this are: the limited range of length scales (low resolution) of the computer simulations; the models do not include the variability in the conditions under which travertine terraces are formed; complexities in the natural growth processes, including those mentioned above, are not fully taken into account; and the approximations that are used to reduce the computational burden. One of the most important approximations is the use of simple two-dimensional (depth-averaged) flow equations instead of the time-dependent, nonlinear, three-dimensional Navier–Stokes equation to simulate the fluid dynamics.

The reasons why natural terraces vary substantially in detail are not well understood. More complex, and much more computationally expensive, models combined with extensive field investigations and laboratory experiments would be required to address this issue. Other details such as the complex shapes of individual terrace edges and the growth of ‘stalactites’ from the overhanging lips of some edges would also require much higher resolution and much more computationally expensive simulations with more physical content and less empiricism.

Hammer *et al*. (2008) used a two-dimensional numerical model for coupled aqueous chemistry, precipitation, fluid flow and chemical transport in an open channel to investigate calcite precipitation. This model reproduces precipitation rates observed experimentally and also explains why the precipitation rate often correlates with the flow rate, even though the flow rate is not strictly the controlling factor.

The morphology of most fluvial systems is controlled by erosion, sediment transport and sediment deposition, and erosion is fastest where the water flow velocity is highest and the bed slope is large. In systems that are dominated by mineral precipitation, the rate of precipitation is greatest under conditions that would lead to erosion in most sediment-dominated systems. In terraced systems, the coupling between fluid dynamics, solute transport and precipitation results in unstable growth (nascent terrace lips induce flow conditions that cause them to grow more rapidly than their surroundings, and this results in flow conditions that cause them to grow even faster), which is captured in the models. However, under rapid growth conditions, the deposited carbonate may be a soft mud-like material (Richards 1950), and erosion channels can be observed in this carbonate mud. This material may be formed by the growth and deposition of calcite particles, possibly mediated by micro-organisms. However, the material deposited on the terrace rims is much harder and more erosion-resistant.

## 3. Growth from slowly flowing liquid films

Stalactites, stalagmites, flowstones, drapery, needles and many other patterns in carbonate caves have fascinated mankind for millennia and have been the subjects of scientific investigation for centuries (Ravbar 2003). More recently, interest has grown in the investigation of speleothems as records of climate change (Gascoyne 1992; Baker *et al*. 1993; Bar-Matthews *et al.* 1997; McDermott *et al.* 2001; Fairchild *et al.* 2006). Speleothems are formed as a result of the coupling between fluid flow, aqueous chemistry, heterogeneous chemistry and the exchange of carbon dioxide (and sometimes other reactive gasses) between the aqueous phase and air. In this respect, the growth of these distinctive patterns is similar to the growth of terraces, as described in §2. The most important difference is that the water flows as a thin film over the surface of the solid under low Reynolds number, laminar flow conditions, as opposed to the turbulent flow over at least parts of the terraces discussed in §2. This substantially simplifies theoretical analysis and numerical simulations.

The growth of icicles is another familiar process that is very similar to the growth of stalactites (figure 4), and ice can also form a variety of other patterns that resemble speleothems such as stalagmites and flowstones. The ice grows from a thin film of water that is in contact with air, which has a temperature less than 0^{°}C. The argument that micro-organisms may play an important role in the formation of mineral terraces, but are not essential because similar patterns form rapidly from freezing water, applies equally to stalactites and many other speleothems. The formation of ice patterns is particularly simple because the growth of the ice (from essentially pure water) is controlled solely by the transport of heat across the water film and in the surrounding air. This simplicity, and the rapid growth of ice, makes it a good material for laboratory experiments and a good starting point for the development and evaluation of theoretical analysis and computer simulations that can be applied, with further development, to the formation of speleothems. However, even this very simple system can exhibit complex behaviour (figure 4*c*).

Szilder & Lozowski (2000) developed a three-dimensional model for icicle growth, which is based on a cubic lattice random walk algorithm. Random walks are initiated at the base of the icicle and follow random walk paths on the surface of the growing ‘icicle’, with a downward bias. At each step in the random walk, the particle (representing a small volume of water) moves to one of the six nearest-neighbour sites on the lattice (providing that it is a vacant site on the surface of the growing icicle) or it freezes. In the freezing step, all of the vacant surface sites within a small distance (the freezing range) are located, and the particle is moved randomly to one of the vacant surface sites with the largest coordinate number, and it becomes an ‘ice’ site (or particle). When the particle reaches the tip of the icicle, it remains for a preselected number of steps before it is removed to represent dripping. This model generates ‘numerical icicles’ that resemble real icicles, and a similar approach could be used to simulate stalactites. However, the effects of the randomness associated with the random walk model are quite pronounced, and large-scale simulations or noise reduction (Tang 1985), both of which would increase the computational burden, would be required to reduce the effects of noise to realistic levels. In addition, the model does not accurately represent the flow of the water film on the icicle surface, ripples do not form and the tips of the simulated icicles (figure 5) are much blunter than the icicles shown in figure 4.

Short *et al*. (2006) have developed a model for the growth of icicles that captures the physics of icicle formation in more detail than the model of Szilder & Lozowski (2000). It is assumed that an aqueous film with a thickness of *h*(*z*), where *z* is the vertical coordinate, covers the stalactite. Under the low Reynolds number laminar flow conditions that are typical of icicle growth conditions, the film thickness is given by
3.1
and the liquid velocity at the liquid–gas interface is given by , where *Q*(*z*) is the volumetric flux, *ν*_{w} is the kinematic viscosity of water, *R*(*z*) is the radius at height *z*, *g* is the acceleration due to the Earth’s gravitational field and *θ*(*z*) is the angle of inclination at height *z*. Equation (3.1) and the equation for *V* (*z*) are derived assuming that *h*≪*R* (the flow velocity profile is the same as that of a thin film of water flowing on a flat inclined surface). For stalactite growth, the rate of evaporation (or condensation) is small, and the volumetric flux can be assumed to be constant. For the icicle problem, on the other hand, the water slowly freezes at the ice–water interface, and the volumetric flux decreases with decreasing height, *z*. The thermal Peclet number, *P**e*=*V* (*z*)*h*(*z*)/D_{T}, where D_{T} is the thermal diffusivity of water, lies in the range 0.1≤*P**e*≤1, indicating that thermal diffusion across the water layer dominates advection of heat down the icicle. Consequently,
3.2
and
3.3
where *κ*_{w} is the thermal conductivity of water (*κ*_{w}=*ρ*_{w}*C*_{p}D_{T}, where *ρ*_{w} is the density of water and *C*_{p} is the specific heat), *T*_{m} is the melting point of ice, *T*_{i} is the temperature at the water–air interface and *L* is the latent heat of melting. In terms of the arc length, *s*, measured from the tip of the icicle along its surface in a vertical plane,
3.4
where *v*_{g} is the growth velocity in the direction normal to the surface of the ice. For icicles, the growth is limited by transport of heat through the laminar air sublayer adjacent to the water film, and the temperature at the water–air interface, *T*_{i}, is determined by the melting point of ice, the air temperature and the air flow conditions through the requirement that the heat flux density must be equal on both sides of the water–air interface (*κ*_{w}(*T*_{m}−*T*_{i})/*h*=*κ*_{a}(*T*_{i}−*T*_{a})/ℓ, where *T*_{a} is the ambient air temperature, *κ*_{a} is the thermal conductivity of air and ℓ is the thickness of the air laminar sublayer). Short *et al*. analysed the case of convection driven by the heat transport from the icicle to the adjacent air in otherwise stationary air, based on the analytical result ℓ=*C**λ*(*z*/*λ*)^{1/4}, for the laminar sublayer thickness adjacent to a flat isothermal vertical plate (Schlichting 1979), where *z* is the vertical distance from the bottom of the plate, *C* is a constant that depends on the Prandtl number, Δ*T* is the difference between the temperature of the plate and the ambient air temperature and *λ* is given by , where *ν*_{a} is the kinematic viscosity of air and *α*_{V} is the volumetric expansion coefficient of air.

Since *κ*_{w}/*h*≫*κ*_{a}/ℓ, the temperature, *T*_{i}, at the water–air interface is essentially constant, and this justifies the use of a constant *T*_{i} to estimate the thickness of the laminar sublayer. By dividing the heat flux density estimated from the thickness of the laminar sublayer (using the isothermal vertical plate model) by the latent heat of freezing, *L*, of water, Short *et al*. derived the equation
3.5
for the growth velocity, *v*_{g}, at the ice–water interface, a vertical distance of *z* above the tip of the icicle. Here is the characteristic velocity,^{1}Δ*T*=*T*_{a}−*T*_{i} and *z* is the vertical distance from the tip of the icicle. These ‘thin film’ equations cannot be used at the tip of the icicle. Instead, Short *et al*. assumed that the icicle has a uniformly translating shape, which requires that
3.6
at every position on the growing icicle, where *v*_{t}, the growth velocity at the tip of the icicle, is a model parameter and *ϑ*=*π*/2−*θ* is the angle between the vertical direction and the surface. Substitution of this into equation (3.3) gives
3.7
which implies that *Q*=*Q*_{t}+*π**R*^{2}*v*_{t}, where *Q*_{t} is the rate at which water drips off the tip of the icicle. To obtain a solution, the expression from equation (3.6) is substituted into equation (3.5) giving
3.8
where . Equation (3.8) with the trigonometric identities and , where , gives ((1+(∂*R*′/∂*z*′)^{2}). This can be rearranged to , which has the solution
3.9
For large *z*′, far from the tip, .

A similar analysis for growth limited by thermal diffusion through the liquid film gave the same asymptotic icicle shape. This theoretical model, which appears to work well (figure 6*a*), is based on the assumption that the growth of the icicle results in a constant shape that translates with the tip of the icicle. This idea has a long history (Ivantsov 1947), and it has been used extensively in theoretical work on dendritic growth (e.g. Brener 1993) and the growth of stalagmites (see below).

Typically smaller scale patterns are superimposed on the overall shapes of icicles (figure 4). These shape perturbations grow because heat can be removed more rapidly from a bulge in the surface of an icicle that is more exposed to the surrounding cold air, and consequently icicle growth is subject to a linear instability similar to the Mullins & Sekerka (1963*a*,*b*) instability.

However, the Gibbs–Thompson effect is negligible on the scale of the ‘surface ripples’ that can be clearly seen in figure 4*c* and to a lesser extent in figure 6*a*. Instead, the rate at which short-wavelength perturbations grow is suppressed by the advection of heat by the moving water film. Ogawa & Furukawa (2002) performed a linear stability analysis for the growth of surface waves on icicles, based on diffusive heat transport and laminar film flow, with a quasi-stationary assumption for heat transport in air (∇^{2}*T*=0). Unlike the Mullins–Sekerka instability, where the Gibbs–Thompson effect results in the decay of short-wavelength perturbations, perturbations at all wavelengths are unstable for the model used by Ogawa and Furukawa. However, shape perturbations with a wavelength of , given by
3.10
grow most rapidly. Here, is the thermal diffusivity of water (*ρ*_{w} is the density of water and *c*_{w} is its heat capacity), and . According to this analysis, the characteristic ripple wavelength depends on the ratio *Q*/*R*, and it does not depend on the icicle growth rate, which is controlled by the far-field air temperature. For typical values of the flow rate, *Q*, and radius, *R*, the wavelength predicted by equation (3.10) is of the order of 1 cm, and this is consistent with field observations (Ogawa & Furukawa 2002). Ogawa and Furukawa argued that because the flowing water film transports heat downwards, along the surface of the growing icicle, the maximum temperature gradient at the ice–water interface is shifted downward from the position of maximum (convex with respect to air) curvature, and the minimum temperature gradient is also shifted downward from the position of minimum curvature, and their model also predicts that the ripples migrate downward with a velocity that is about half the growth rate, normal to the surface.

If the rate of icicle growth is limited by the transport of heat through the laminar sublayer in the cold air that is convected upwards because of warming by the latent heat of freezing released by the growing icicle, a similar instability would be expected. The diffusive model used by Ogawa & Furukawa (2002) to represent the heat transport in air is quite different from the boundary layer model used by Short *et al.* (2006). In addition, a quite different perspective on the origin of ripples on icicles has been offered by Ueno (2007). The linear stability analysis of Ueno takes into account the effects of the surface tension at the water–air interface, as well as heat transport and fluid flow in the liquid film. This analysis leads to the conclusion that the wavenumber of maximum growth is given by
3.11
where *a*=[*Γ*_{aw}/*ρ*_{w}*g*] is a capillary length (*Γ*_{aw} is the air–water surface tension) and *P**e*=*u*_{0}*h*_{0}/*κ*_{w} (*h*_{0} is the mean water film thickness and *u*_{0} is the water film velocity at the air–water interface) is the Peclet number. According to this analysis, the characteristic ripple wavelength increases with increasing surface tension and with increasing flow rate. The ripples move upwards (in contrast to the downward motion predicted by Ogawa & Furukawa (2002)), which appears to be consistent with results from experimental investigations (Maeno *et al*. 1994). Both linear stability analyses predict ripple wavelengths of about 1 cm. A schematic of a vertical cross-section through an icicle in a theoretical paper published by Ueno (2004) shows ‘V’-shaped bands of tiny air bubbles, formed by trapping of bubbles in the water film as the freezing front advances, that emerge from the top part of each ripple, and the author points out that this downward pointing chevron pattern is consistent with upward motion of the ripples. This pattern appears to be quite well known, but not well documented in the scientific literature.

Like icicles, stalactites are also covered with a thin film of water, which creeps (equation (3.1)) slowly downwards precipitating CaCO_{3} as CO_{2} diffuses from the CaCO_{3}–water interface through the water film and into the cave atmosphere via the water–air interface. Based on this conceptual model, Short *et al*. (2005) analysed the growth of stalactites. The process is substantially more complex than the growth of an icicle, because of the coupling between the aqueous carbonate chemistry and the transport of CO_{2.} In the humid atmosphere found deep in temporate caves, the water film does not evaporate, and only a small fraction of the dissolved calcium is precipitated (except at low drip rates, stalagmites often have a much larger volume than the volume of the stalactite that feeds them (Baldini 2001)). Since aqueous carbonate chemistry is quite complex, and precipitation may be controlled by the kinetics at the mineral–water interface, the interconversion of CO_{2} and H_{2}CO_{3} and transport of dissolved species (Buhmann & Dreybrodt 1985), it is not immediately clear how the coupling between chemistry and transport can be simplified to a growth velocity law that can be used to simulate stalactites. Fortunately, the characteristic time *τ*_{D}=*h*^{2}/D_{w} for diffusion across the thin (tens of microns) water film is small enough to ensure that diffusion through the water film is not rate-limiting. The derivation of the growth velocity law (Short *et al*. 2005) is too lengthy to reproduce here. Based on the well-justified ideas that diffusion in the cave atmosphere has little influence on the growth rate, that the time that the water remains in contact with the stalactite is much longer than the diffusion time and that the growth time is much longer than the contact time, Short *et al*. (2005) derived the growth equation
3.12
where *V*_{n} is the growth velocity normal to the interface, *v*_{m} is the molar volume of CaCO_{3}, *l*_{Q}=[3*ν*_{w}*Q*/(2*π**g*)]^{1/4} is the characteristic length and R_{c} is an effective rate that depends upon the aqueous chemistry. Using an approach very similar to that employed to develop a shape equation for icicles, Short *et al*. derived the approximate equation
3.13
where is the scaled height and is the scaled radius (*v*_{t} is the tip velocity), which quite accurately describes the shape of simple stalactites (figure 6). Figure 6(ii) compares the theoretical stalactite shape with that of stalactites from Kartchener Caverns, Benson, AZ, USA.

Modelling and simulation of stalagmite growth is also facilitated by the simplicity of thin-film hydrodynamics. The growth of a pendant drop of water in the tip of a stalactite and its eventual release and impact on the top of the growing stalagmite complicates the story, but has relatively little effect on the outcome. Theoretical models for the growth of stalagmites are based on the idea that water spreads out uniformly from the centre of the stalagmite, and the rate of precipitation decreases with increasing distance from the centre. It is clear from the shape of the stalagmite that the growth rate is highest near its centre and decreases with distance from the top. Kaufmann (2003) developed a model for stalagmite growth based on the idea that the growth rate is given by
3.14
where *s* is the distance from the centre of the stalagmite measured along its surface, *V*_{0} is the maximum growth rate at the top of the stalagmite and *R*_{0} is the ‘equilibrium radius’ of the stalagmite. The motivation for this exponential form is not clear. In the case of precipitation on the surface of a stalagmite, a model that takes into account the spreading and thinning of the water film with increased distance from the top of the stalagmite should be used, and the resulting precipitation rate is better approximated by a Gaussian than an exponential function of *s* (Romanov *et al.*2008*a*,*b*). In any event, this model predicts that the stalagmite forms a column with vertical sides and a rounded top, provided that the growth conditions do not change. It turns out (Romanov *et al.*2008*a*,*b*) that, as anticipated by Franke (1965), the shape of the surface evolves towards a constant form that translates upwards at a constant rate for any normal growth rate that decreases monotonically with distance, *s*, from the centre measured in the curvilinear surface coordinate system. If the growth rate is determined by the function *f*(*s*), then the vertical growth rate is given by
3.15
where *s*(*R*), the arc length from the centre of the stalagmite to the locus of points at a distance *R*, from the vertical axis through the centre of the stalagmite in the Cartesian coordinate system, is given by
3.16
and *θ*(*r*) is the angle of inclination of the surface (the slope *m* is given by .

To determine the shape of the stalagmite, an expression that relates the slope ∂*z*(*R*,*t*)/∂*R* to the radius *R* is required. If the shape of the stalagmite evolves towards a constant form that moves upwards with a constant velocity ,
3.17
or and it follows from equation (3.16) that
3.18
or
3.19
which is the desired shape equation. Of course, the function *f*(*s*(*R*)) would change as the shape of the stalagmite changes for a real stalagmite. Romanov *et al*. based their argument for a self-preserving asymptotic shape that translates upwards with a constant velocity on the idea that *f*(*s*(*R*)) is independent of time. They argue that the arc length, *s*(*R*), for some arbitrary radius, *R*, small enough for *f*(*s*(*R*)) to be non-zero must decrease with increasing time because the growth is normal to the surface and *f*(*s*(*R*)) decreases with increasing distance, *s*(*R*).

The evolution of the surface of the stalagmite is described by the equations (Ben-Jacob *et al*. 1983; Brower *et al*. 1983)
3.20
and
3.21
where ∂*κ*/∂*t*|_{n} indicates the time derivative of the curvature along the outwardly directed normal. In equations (3.20) and (3.21), the curvature, *κ*, is positive for parts of the surface that are convex viewed from outside the stalagmite and negative for concave parts of the surface. The first term on the right-hand side, −*κ*^{2}*V*_{n}, decreases the curvature due to uniform dilation (in convex parts of the surface) or contraction (in concave parts of the surface) and the second term describes changes in the curvature due to local variations in the normal growth velocity. For constant (position-independent) *V*_{n}, the first term increases the magnitude of the curvature in concave parts of the interface, eventually leading to the formation of grooves, and convex parts of the surface grow at the expense of concave parts. In general, ∂^{2}*V*_{n}/∂*s*^{2} can have a complex shape, depending on the function *V*_{n}(*s*), but if *V*_{n}(*s*) decreases monotonically with increasing *s*, then is negative and, on average, the term −∂^{2}*V*_{n}/∂*s*^{2} increases the curvature. It is the interplay between the first and second terms on the right-hand side of equations (3.20) or (3.21) that determines how the interface will evolve.

Romanov *et al*. (2008*a*,*b*) performed simulations with three models: (i) a growth velocity that decreases exponentially with the arc length, *s*(*R*), (ii) a growth velocity that was a Gaussian function of *s*(*R*), and (iii) a model that assumes that *V*_{n}=*k*(*c*−*c*_{eq}), where *c* is the calcium concentration at the calcite–water interface, and *c*_{eq} is the equilibrium calcium concentration, which depends on the temperature and the partial pressure of CO_{2} in the water. For model 3, a constant flow rate of *Q* equals *v*_{d}/*τ*, where *v*_{d} is the volume of each drop that falls on the top of the stalagmite and *τ* is the interval between drops that was used. Figure 7 shows the results of simulations performed with sinusoidally varying drip rates. In the simulations, the kinetic coefficient *k* is 1.3×10^{−7} m s^{−1}, which gives a vertical growth rate of about 1.1×10^{−2} cm yr^{−1} at the apex, and each drop has a volume of 0.1 ml. For a periodicity of 500 years, the effects of the variable drip rate are not discernable, and the stalagmite has a simple columnar ‘candlestick’ shape. For much longer periods, the morphology of the stalagmite records the variability in the growth conditions. Simulations with precipitation rates that decrease exponentially or are Gaussian functions for the distance, *s*, from the top of the stalagmite, as well as simulations based on a model that explicitly takes into account film flow on the top of the stalagmite (figure 8*b*), generate columnar forms that resemble the shapes of ‘candlestick’ stalagmites. Figure 8 shows stalagmites with this shape that vary in volume by a factor of 10^{7}–10^{8}. Nothing is known about the internal structure of the very large stalagmite shown in figure 8*a* or the conditions under which it grew. It does not seem likely that it grew from a single dripping source at a more or less constant rate, and it is not clear whether the model can be applied to such large stalagmites. However, the fact that the predicted shape is in good agreement with the most simple shape observed in real stalagmites suggests that the more complex forms are due to changing growth conditions rather than complexity in the growth process. These complexities include changes in the cave atmosphere, dripping rate, temperature, water composition, etc. In addition, the position of the water source may change over time, and there may be more than one source.

The top of the stalagmite shown in figure 8*c*, and the tops of the internal laminations, which record the shapes of the stalagmite at earlier stages of growth, are less convex than the model stalagmite shown in figure 8*b*, and they more closely resemble the model results shown in figure 7*a*. The convexity of cylindrical stalagmites varies substantially (many stalagmites have perfectly flat tops), and additional conceptual and computational modelling will be needed to better understand why stalagmite shapes vary so much.

Some stalagmites have an approximately conical shape (the radius of stalagmites often decreases upwards), or the radius may increase with increasing height, and these shapes can be modelled in terms of a steadily decreasing or increasing supply of water (Dreybrodt 1999). However, conical stalagmites are more commonly formed under intermittent high-flow-rate conditions (Miorandi *et al*. in press).

Ice may also form stalagmites that resemble carbonate stalagmites. If the growth of ice stalagmites, like carbonate stalagmites, is concentrated on the more or less hemispherical dome at the top, this would be expected. Although some ice stalagmites have complex shapes, the complexity appears to be related to variable growth conditions rather than growth instability. It is apparent that stalagmites that grow from thin films of water (under the flow conditions assumed in the model of Romanov *et al.* (2008*a*,*b*)), whether they are formed from calcium carbonate or ice, are not susceptible to the rippling instability associated with icicles, because only the upper part of the stalagmite is growing. Ripples are commonly observed under more rapid flow conditions, and film-flow hydrodynamic instabilities (Liu *et al.* 1995) may play a role in their growth.

A variety of ‘cones’, ‘pyramids’ and ‘columns’ may be formed by the deposition of travertine near hot springs. These features are frequently surrounded by travertine terraces, and complex patterns with characteristic scales of the order of a few centimetres and of the order of a few decimetres are often found on their steeply sloped sides. The shapes of these features likely depend on flow conditions. Travertine columns may form when the flow of water is insufficient to cover the sides of the column, and growth is restricted to its top. In other cases, the supply of water is large enough for the water to spurt from the top. Under these conditions, complex mounds form because there is a strong asymmetry and variability in the way in which water laden with CO_{2} and dissolved calcium is deposited on the travertine. The travertine Liberty Cap (figure 9*a*) was most likely formed under intermediate flow conditions. Unfortunately, it was dormant by the time that the Hayden Expedition surveyed the area in 1871, and the flow conditions are not known. Jettestuen *et al*. (2006) have simulated the growth of Liberty Cap using a ‘normal growth model’ quite similar to the exponential and Gaussian models of Romanov *et al.* (2008*a*,*b*). The main difference between the two models is that the model of Romanav *et al*. is based on the idea that the growth rate is an exponentially decreasing or Gaussian function of the distance, *s*, from the top measured along the surface, while, for the model of Jettestuen *et al*. (2006), the normal growth velocity is given by *V*_{n}=1/(1+*c**s*), where *c* is a constant. Figure 9*b* compares the shape of the 11 m high Liberty Cap with the shape obtained from a simulation performed with the model of Jettestuen *et al*. (2006) with *c*=0.07 m^{−1}.

Goldenfeld *et al.* (2006) have simulated the growth of travertine mounds assuming that the normal growth velocity, *V*_{n}, is proportional to the depth-averaged tangential fluid flow velocity, and the velocity profile is obtained from a laminar flow model. Good agreement between the model and the shape of a travertine dome at Mammoth Hot Springs, Yellowstone, was found for the top part of the dome (figure 9*c*). The departure between the model predictions and the actual shape further down the dome coincides with the appearance of a fluting pattern around the dome, which Goldenfeld *et al*. attribute to the effects of surface tension, and the onset of partial wetting, with the appearance of a water–air–travertine contact line. This is a complex problem that cannot be easily addressed analytically. However, simulations that allow flow only if the water depth measured in the direction of the surface normal (the thickness of the water layer) exceeds a threshold to represent the effects of surface tension do generate a shape that is in good agreement with the observed shape of the dome.

## 4. Scallops

A wide variety of patterns are formed by chemical and/or physical erosion of carbonates and other relatively soluble rocks in karst systems (e.g. Ford & Williams 2007). The formation of a densely packed array of ‘scallop’-shaped depressions on the walls of phreatic cave passages (figure 10) is a striking example, which is important because the characteristic length, *l*, of the scallops is often used to estimate the velocity of the water based on the idea that the pattern evolved towards a critical Reynolds number (e.g. Goodchild & Ford 1971; Curl 1974; Thomas 1979; Gale 1984). This illustrates how a partial but quasi-quantitative understanding of complex patterns can sometimes be obtained from a simple conceptual model and scaling (dimensional) analysis. The Reynolds number is defined in terms of the molecular viscosity, *η*, or kinematic viscosity, *ν*=*η*/*ρ*, of water, the scallop length, *l*, and a characteristic velocity. In some investigations, the mean flow velocity, *V* , has been used as the characteristic velocity and the Reynolds number is defined as *R**e*=*ρ**V**l*/*η*=ℓ*V* /*ν*.

Curl (1966) applied classical scaling analysis to flutes and scallops based on the idea that the relevant dimensional variables are the average fluid flow velocity, *V* , the fluid density, *ρ*, the viscosity, *η*, the channel dimension, *L*, the dimension of the flute or scallop, *l*, and the diffusion coefficient of solute ions, *D*. Three independent dimensionless ratios can be constructed from these variables, and one choice is *Π*_{1}=*R**e*=*ρ**V**l*/*η*, *Π*_{2}=*l*/*L* and *Π*_{3}=*S**c*=*η*/*ρ**D*. These dimensionless ratios are related by *F*(*Π*_{1},*Π*_{2},*Π*_{3})=0 or *Π*_{1}=*f*(*Π*_{2},*Π*_{3}), where *f* is a characteristic function that describes the behaviour of the system. Curl argued that if the length ratio, *Π*_{2}, is sufficiently small and the Schmidt number, *Π*_{3}, is sufficiently large that , where *R**e*^{c} is the critical Reynolds number for stable flutes or scallops, and for flutes in water–carbonate and air–ice systems a value of *R**e*^{c}≈22 500 was estimated. The Reynolds number, *R**e*^{L}=*ρ**V**L*/*η*, based on the channel size is much larger than the critical Reynolds number based on the flute size, and this indicates that the channel flow is strongly turbulent.

Goodchild & Ford (1971) reported scallop lengths varying from 2 mm in steeply sloped conduits to 2 m in large horizontal cave passages. They also investigated the formation of scallops in blocks of plaster of Paris that were cast on the bottom of a 2 m wide flume tank. The blocks were 65 cm long and they were wedge-shaped with a height of ≈0 on the upstream edge rising with a constant slope to 6 or 10 cm on the downstream edge. The results were consistent with a critical Reynolds number (based on the scallop length) but with a smaller critical Reynolds number (*R**e*^{c}≈10 000).

Subsequent investigators (Blumberg & Curl 1974; Thomas 1979) have defined the Reynolds number in terms of the friction velocity, , where *ρ* is the fluid density, *l* is the characteristic scallop length (the ratio between successive integer moments of the scallop length distribution (*M*^{(3)}/*M*^{(2)}, for example, where *M*^{(n)} is the *n*th moment, may be used for the characteristic scallop length), *η* is the dynamic viscosity of the fluid and *u*_{*} is the friction velocity defined as , where *σ*_{w} is the wall stress. This definition of the Reynolds number is based on the wall similarity hypothesis (Townsend 1980), and it provides a more consistent way of comparing scallop formation in different systems. Since the friction velocity is smaller than the mean flow velocity, the critical Reynolds number is smaller than *R**e*_{c}, and typical values reported in the literature lie in the range (Blumberg & Curl 1974; Thomas 1979). Curl estimated a critical Reynolds number of for a small range (about 3.15) of scallop lengths based on experiments in which plaster of Paris was dissolved in flowing water, while Thomas (1979) found that data from a variety of sources could be represented reasonably well by *l*=1000 *η*/*ρ**u** for a number of scallop pattern-forming phenomena over about five orders of magnitude in *l* and *ρ**u**. The uncertainties in the critical Reynolds numbers, , obtained by Blumberg and Curl and by Thomas are quite large, and additional experimental work would be required to narrow the range of uncertainty.

Although it is still often assumed that the mean fluid velocity is inversely proportional to the scallop length, the mean velocity can be more accurately estimated from the scallop length using semiempircal relationships such as , where *L* is the channel diameter and *C* is a coefficient related to the roughness (Blumberg & Curl 1974). Most likely, cave scallop patterns are formed by dissolution or corrasion. However, the formation of scallops does not appear to depend on the microscopic mechanisms that control how solid is removed. Similar patterns are observed due to melting in water-filled ice caves (figure 11*a*), on the undersurface of ice in rivers, on meteorites (Lin & Qun 1986; figure 11*b*), re-entry vehicles and impact ejecta surfaces (chemical and thermal ablation; Ernstson (2004); Duffa *et al*. (2005)), on subliming ice (Bintanja 1999; Bintanja *et al*. 2001), on snow on the Earth and Mars, on cohesive bed forms in fluvial systems (Allen 1971) and on the surfaces of corroding metal pipes (Villien *et al.* 2001). The formation of scallops on re-entry vehicles and their effects on heat transfer were investigated as part of the Passive Nosetip Technology (PANT) Program in the 1970s (Derbridge & Wool 1974; Wool 1975) and other large research programmes related to weapon development and space applications. This work is not readily available, and it is not clear whether it led to results or understanding that could be applied to cave scallops or other geological scallop patterns.

Very likely it is no accident that the critical Reynolds number, , associated with scallop patterns is essentially the same as the critical Reynolds number for the onset of turbulence in rough-walled channels. Near to the onset of turbulence, at relatively small Reynolds numbers, the fluid vortex patterns are more regular than they are at higher Reynolds numbers. Conceptual models for the formation of scallops envisage a single vortex within each scallop (Blumberg & Curl 1974). Figure 12 shows a cartoon of the flow behaviour in a scallop. A wide variety of other patterns are formed by dissolution, corrosion, corrasion, erosion and the deformation of weakly cohesive sediments mediated by turbulent flow (Allen 1971). Because of the complexity of turbulent flow, it is still very difficult to simulate these pattern formation phenomena in quantitative detail. In most cases, the pattern formation process takes place under constantly varying conditions, and in some cases the patterns formed under these conditions may differ in important ways from patterns formed in well-controlled laboratory experiments or in numerical simulations, which usually do not include varying conditions for practical reasons.

## 5. Spherulitic growth

The formation of mineral bodies that have a more or less spherical or botryoidal form (resembling a bunch of grapes) with a subradially oriented polycrystalline fibrous internal structure is common in a wide variety of settings. Although first investigated by geologists in the nineteenth century (Judd 1888; Cross 1891), they have been studied most extensively by polymer scientists since essentially all polymers cooled from a molten state form spherulites (Keith & Padden 1963). The internal structure consists of an array of crystalline fibers that branch at small angles, which do not appear to be crystallographically controlled. Spherulites are polycrystalline, indicating that repeated ‘secondary’ nucleation is required for their formation. This fine-scale polycrystallinity leads to a low-macroscopic anisotropy and a more or less spherical large-scale morphology. Impurities are present in almost all systems in which spherulitic growth has been found (impurities in the form of a distribution of molecular weights, differing copolymer sequences and differing side-branch lengths and positions along the main chain are unavoidable in synthetic polymers), and, in many cases, spherulites grow from very viscous fluids in which diffusion coefficients are small. The observation that impurities are almost always present during spherulitic growth has led to the ideas that it is controlled by the rejection of impurities, and that the highly branched fibrous micromorphology has an origin similar to that of the fibrous and lamellar morphologies associated with unidirectional eutectic solidification (Hunt & Jackson 1966) and eutectoid decomposition.

Spherulites grow outwards from a nucleus (figure 13), branching to fill the available space or they first form as fibers that evolve into a sheaf (like a sheaf of corn; figure 13*b*) that continues to grow and spread until it becomes spherical. Figure 14 shows a conceptual model for the impurity-controlled growth of a spherulite in a quasi-infinite homogeneous medium, assuming constant growth conditions. During the early stages of growth, the impurity rejected by the growing crystalline material, which can be considered to be essentially pure, becomes part of an impurity boundary layer, which is larger than the size of the incipient spherulite. During this stage, the radius, *r*, grows as *r*∼*t*^{1/2}, where *t* is the time after growth started. This transient stage gives rise to a crossover in which the radius and the width of the boundary layer are of comparable size followed by an asymptotic regime in which the interface propagates with a constant velocity, and *r*∼*t*. In this asymptotic regime, the radius is much larger than the width of the boundary layer, and the characteristic morphology length (the distance between the centres of adjacent crystalline fibers), and the propagation of the growth front is essentially equivalent to the propagation of a flat front. In this asymptotic regime, the growth is faster than would be expected from the initial, *r*∼*t*^{1/2} behaviour, because most of the impurity rejected by the crystallization process remains in the region between the growing fibers. Long after the growth front has passed, the fibers may continue to broaden, and the concentrated impurity may become incorporated into impure phases or diffuse away under different conditions. At the stage depicted in figure 14, the width of the diffusion boundary layer, *δ*=D/*V* , where D is the impurity diffusion coefficient and *V* is the growth velocity, is much smaller than the radius, *r*, of the growing spherulite. In the early stages of growth, when *r*≪*δ* the distance that the impurity diffuses through in time *t*, (D*t*)^{1/2}, is larger than the radius of the spherulite.

The morphology scale (*λ* in figure 14) is determined by the Mullins–Sekerka instability (Mullins & Sekerka 1963*a*,*b*; Goldenfeld 1987). The Mullins–Sekerka instability occurs because the growth rate is greatest at the most exposed parts of the interface that have already grown the furthest, where the impurity concentration is smallest due to diffusion into the surrounding liquid. Short-wavelength perturbations to an initially flat circular interface grow more rapidly than long-wavelength perturbations because impurities diffuse more rapidly through short distances than large ones. However the effects of surface energy inhibit the growth of short-wavelength perturbations (through the Gibbs–Thompson effect) and this effect becomes more dominant as the wavelength decreases. The ‘competition’ between these two effects results in a characteristic wavelength at which perturbations grow most rapidly, and this sets the morphology scale, *λ*. These ideas are expressed mathematically, in the analysis of Mullins & Sekerka (1963*a*,*b*).

Spherulitic (botryoidal) morphologies are often formed by mineral growth in systems where micro-organisms are abundant. In these cases, biofilm or biopolymers may serve as the impurity, and there is evidence that polysaccharides and amino acids promote spherical growth (Braissant *et al*. 2003). In these systems, the large increase in viscosity and small ‘impurity’ diffusion coefficients may play an important role (Buczynski & Chafetz 1991). Observation of spherulitic growth in quite pure materials (e.g. Bisault *et al*. 1991) has called into question the need for impurities in spherulitic growth. However, significant impurity concentrations are present in essentially all natural systems, and the formation of spherulites in the absence of impurities does not preclude their importance in most systems of geological interest. Very few examples of the growth of spherulites without impurities are known, and the growth mechanisms with and without impurities may be quite different.

Granasy *et al*. (2005) have argued that spherulites can be formed via a variety of mechanisms that lead to ‘secondary’ nucleation at the crystallization front. These mechanisms include growth front nucleation due to the quenching of orientational defects that may arise from impurities or long-lived dynamic heterogeneities in supercooled liquids. In particular, they argue that the ratio between the rotational and translational diffusion coefficients is very small in highly supercooled liquids, and that this promotes the formation of misoriented crystallites if the rate of molecular reorientation is slow relative to the rate at which the solidification front propagates (Granasy *et al*. 2004). Branching with a preferred crystallographic misorientation is an additional mechanism that Granasy *et al*. (2005) propose to explain spherulitic growth in pure materials with a small supercooling. Based on these ideas, Granasy *et al*. (2005) have developed a phase field model that includes all three mechanisms (quenching of impurity defects, quenching of intrinsic orientational defects and branching with a preferred crystallographic misorientation) modelled as branching with a fixed crystallographic misorientation with respect to the orientation of the parent crystallite that provides the nucleation site and implemented via an orientation-dependent grain boundary energy. The free energy functional depends on a phase field that characterizes the liquid/solid transition, the chemical composition, the temperature field and the crystal orientation in a fixed coordinate system. The phase field model is flexible enough to simulate the growth of spherulites that grow radially from a single nucleus and spherulites that grow via a sheaf of fibrils.

In a uniformly cooled melt, many spherulites are nucleated, and they may continue to grow until they fill the available volume and their boundaries form a Voronoi tessellation. If the growth of spherulites depends on sparsely distributed nuclei (on solid surfaces) they may reach a much larger size. The morphologies produced by spherulitic growth depend on the rate of nucleation of new spherulites relative to the rate of growth. Figure 15 shows patterns generated by a model that combines normal growth with random nucleation (Jettestuen *et al*. 2006).

Spherulitic growth and normal growth with random nucleation can lead to the formation of a variety of quite striking patterns, like those shown in figure 16. Although snowflakes (dendrites) are more familiar, ice may also form spherulites (Hey & MacFarlane 1998).

## 6. Dendrites

When the kinetics of the precipitation chemistry at the fluid–solid interface is fast, and the resulting growth process is limited by the diffusive transport of solute to the interface, dendritic structures are formed. The diffusion-limited aggregation (DLA) model (Witten & Sander 1981) illustrates growth in this limit. In off-lattice DLA models, particles (discs in two-dimensional simulations, figure 17*a*, or spheres in three-dimensional simulations) are added, one at a time, to a growing cluster of particles (initially a single particle that represents a ‘seed’ or ‘nucleus’) via random walk paths that start far outside a growing cluster of particles and terminate, with addition of the particle to the growing cluster, when the particles first contact a particle that has already been incorporated into the cluster. In lattice DLA models, the growing cluster is represented by filled lattice sites, and random walk trajectories that terminate when they reach a vacant site adjacent to an already filled site are used to select which site will be filled to represent diffusion-limited growth from a very dilute solution. Simple versions of the model can be used to grow clusters with the order of 10^{4} particles or sites, but quite sophisticated methods (Ball & Brady 1985; Tolman & Meakin 1989) are required to grow clusters with 10^{6} or more particles or sites (figure 17*a*). The successive branching that occurs during the growth of a DLA cluster is a manifestation of the Mullins–Sekerka instability. In most cases, the effect of anisotropy on the growth of dendrites is large, and the randomness inherent in the DLA model, implemented on a regular lattice, is too large to generate dendrites that resemble natural dendrites (figure 18*a*,*b*). Features that have been added to the basic DLA model to make it more realistic include non-zero solute (random walker) concentrations, sticking probabilities smaller than unity, site selective sticking probabilities, surface diffusion and noise reduction (Meakin 1998).

Lattice-based DLA models are sensitive to lattice anisotropy, and some versions (particularly those with noise reduction and site-selective sticking probabilities) generate patterns that resemble single crystal mineral dendrites (compare figure 17*b* with 18*a*). The approximate symmetry of single crystal mineral dendrites is controlled by the orientation dependence of the surface energy and kinetic coefficient(s) relative to the crystallographic axes. Dendritic structures can be generated by a variety of processes in addition to crystal growth. This is illustrated in figure 18*c* (Daccord & Lenormand 1987), which shows dendritic patterns formed by the dissolution of plaster of Paris (hydrated CaSO_{4}). In the experiment illustrated here, water was injected into a cylinder, which was filled with plaster of Paris. After the dendritic dissolution channels had been formed, the system was allowed to dry and Woods metal was injected into the channels. The figure shows the Woods metal casting. Although the process is very different from dendritic precipitation, the underlying physics is similar. The flow of water in the porous plaster of Paris is described by the Darcy equation, **V**=*K*∇*P*, where **V** is the fluid velocity, *K* the permeability and *P* the pressure in the fluid. The conservation equation, ∇⋅**V**=0, can be expressed as ∇*K*∇*P*=0, which has the same form as the diffusion equation in the quasi stationary limit (d*C*/d*t*→0, where *C* is the concentration).^{2} Consequently there is a close relationship between dissolution in the limit in which the water becomes saturated with CaSO_{4} as soon as it enters the porous medium (in the large Damhohler number limit) and diffusion-limited growth. The formation of black ice dendrites in white ice is a similar process (Knight 1987). In this case, relatively warm water flowing through a hole from beneath the ice flows into a layer of slush on top of the ice. The flow of the water through the slush is described by the Darcy equation, and the warm water melts dendritic channels in the slush, which later freeze into black ice. In most cases, dendrites are quite small. However, Haack & Scott (1992) have suggested that dendrites formed during the cooling of asteroids may have grown inwards by hundreds of metres or greater from the base of the mantle. The electric potential is another scalar field that is described by the Laplace equation, and lightning strikes on sandy soil form branched structures composed of fused sand called fulgurites (Newcott 1993). The formation of manganese oxide pseudo-fossils is another geologically relevant example of random dendritic growth.

## 7. Faceted crystals

The formation of small, macroscopically faceted, single crystals due to precipitation from supersaturated solutions or solidification of melts is extremely common. Small facets that appear as crystals are cooled below the roughening temperature (Balibar *et al.* 2005), and, well below the roughening temperature, most of the surface is occupied by facets at equilibrium (at still lower temperatures, these faceted crystals may no longer be at equilibrium). In general, crystals formed by precipitation are not equilibrium structures, and the crystal shape depends on the growth conditions (the nature and concentration of impurities, supersaturation, flow conditions, stress, etc.) (Pimpinelli & Villain 1999). Despite these complications, and the formation of pseudomorphs (Blum 1843) by replacement reactions (Merino *et al.* 1995), the crystal habit is often used in preliminary mineral identification. To form large macroscopically faceted crystals (figure 19), the morphological scales associated with the Mullins–Sekerka instability and similar growth instabilities must be larger than the crystal size, and this requires small growth velocities (a large diffusion length) and/or small supersaturations (large supersaturation capillary lengths). In addition, the nucleation rate must be very small, and this also requires a small supersaturation. Garcia-Ruiz *et al*. (2007) have investigated the giant (up to 11 m long) gypsum (CaSO_{4}⋅2H_{2}0) crystals found in the Naica mine (figure 19*a*). They propose a scenario in which gypsum grows from a calcium sulphate solution that is slightly supersaturated with respect to gypsum and slightly undersaturated with respect to anhydrite (CaSO_{4}), which is the source of the calcium sulphate. The process occurs at a temperature that is a slightly below the transition temperature of about 58^{°}C, above which anhydrite is the most stable form and below which gypsum is the most stable form. The estimated growth rate for these crystals is about 1.5 μm yr^{−1} (http://www.cprm.gov.br/33IGC/1322921.html).

## 8. Discussion and conclusions

Despite the complexity of almost all geological systems, it is often possible to understand at least some aspects of geological pattern formation in terms of quite simple conceptual and mathematical models. Very often the origin of this apparent simplicity is the dominance of a single mechanism or a small number of mechanisms over a limited range of length scales. However, pattern formation processes that can be described and simulated using a common mathematical model are not necessarily controlled by the same mechanisms. This is illustrated by the similarity between patterns formed by the freezing/melting of ice and precipitation/dissolution of minerals. All of the pattern types discussed in this review may be formed by the precipitation/dissolution of minerals or the freezing/melting of ice. In this case, the reasons for the similarity are quite evident—the transport of solute and heat are described by the same equations (under both laminar and turbulent flow conditions) and the relationship between the transport of solute or heat and the growth rate are described by the same equation (the parameters are, of course, different, but key dimensionless ratios may have similar values).

Computer simulation plays an important role in the investigation of geological pattern formation. It may be used to evaluate conceptual models and mathematical models based on them. Similarity between computer-generated and natural patterns increases confidence in the underlying conceptual and mathematical models, but it does not unambiguously validate them. Clearly, quantitative comparison between computer-generated patterns is of more value than qualitative comparisons, and whenever possible the dynamics of the natural pattern formation process should be quantitatively compared with simulation results. Experiments under well-controlled conditions are also important. However, it is often very difficult or impossible to construct properly scaled experimental models of geological systems. On the other hand, experiments can be used to investigate processes that are too complex to simulate numerically. Pattern formation in systems in which turbulent flow plays an important role is a commonly encountered example. The range of time and length scales that can be investigated experimentally is also often substantially greater in experimental systems.

The formation of geological patterns is a consequence of physical and chemical processes that are, at least in general terms, well understood. However, most geological patterns are formed by a combination of processes acting over a very large range of time and length scales. For this reason, and because of the nonlinearity of many physicochemical processes, the origins of relatively few patterns have been explained in detail. It could be argued that a focus on simple patterns diverts attention and resources from more economically important problems such as resource location and recovery or geohazards, or more scientifically important problems such as large-scale geodynamics or the coupling between biological and geological processes. However, the knowledge, understanding and methods developed to investigate simple pattern formation processes will be applied to more complex systems in the future. The scientific investigation of pattern formation develops and tests our capability of expressing underlying geological processes in quantitative mathematical terms. This capability is crucial to prediction of the future behaviour of important and interesting geological systems, which determines the relevance of geosciences for the general public.

## Acknowledgements

This study was supported by a Center of Excellence grant to P.G.P. from the Norwegian Research Council. Comments by Ian J. Fairchild and an anonymous reviewer led to substantial improvements in this review.

## Footnotes

↵1 A characteristic velocity is a convenient and well-defined representative velocity (reference velocity or typical velocity) of a physical system. For example, in fluid dynamics, the maximum velocity, the average velocity or the root mean square velocity may be used. The characteristic values of other dimensional quantities can be defined similarly.

↵2 In both cases the scalar field (concentration or pressure field) obeys the Laplace equation, which can be simulated by random walks with the appropriate boundary conditions (launching and killing positions).

- Received April 9, 2009.
- Accepted October 16, 2009.

- © 2009 The Royal Society