Chemical oscillating patterns are ubiquitous in geochemical systems. Although many such patterns result from systematic variations in the external environmental conditions, it is recognized that some patterns are due to intrinsic self-organized processes in a non-equilibrium nonlinear system with positive feedback. In rocks and minerals, periodic precipitation (Liesegang bands) and oscillatory zoning constitute good examples of patterns that can be explained using concepts from nonlinear dynamics. Generally, as the system parameters exceed some threshold values, the steady (time-independent) state characterizing the system loses its stability. The system then evolves towards other time-dependent solutions (‘attractors’) that may have an oscillatory behaviour or a complex chaotic one. In this review, we describe many of these pattern types taken from a variety of geological environments: eruptive, sedimentary, hydrothermal or metamorphic. One particular example (periodic precipitation of pyrite bands in an evolving sapropel sediment) is presented here for the first time. This will help in convincing the reader that the tools of nonlinear dynamics may be useful to understand the history of our planet.
Rocks and minerals often exhibit rhythmic patterns, whereby their chemical composition and/or physical properties vary repetitively across a spatial transect. The pattern length scales could vary from the micrometre scale to the kilometres. Many such patterns are the result of systematic variations in the external environment that generated the deposit. For instance, varves in lake sediments reflect climate annual variations. However, it is now being recognized that many rhythmic patterns result from self-organized processes in a non-equilibrium nonlinear geochemical system [1,2], i.e. from the intrinsic nonlinear dynamics of the system (with positive feedback between some of the variables) without requiring the presence of a fluctuating external template. The concepts behind nonlinear dynamical systems can thus be used to understand the genesis conditions of many petrological or mineralogical rhythmic patterns.
In this review, we present examples of rhythmic patterns in rocks and minerals that are purely geochemical in origin, emphasizing the variety of geological environments where they are found and their mechanism of formation, both qualitatively and quantitatively (when available). Space limitations will not allow us to discuss patterns that involve hydrodynamic flows or mechanical effects, but a recent review can be found in Jamtveit & Hammer . Section 2 presents the general characteristics of one type of rhythmic patterns (Liesegang bands), emphasizing their mechanisms of formation and modelling strategies. Section 3 illustrates examples of Liesegang patterns taken from igneous, sedimentary, metamorphic and hydrothermal environments. The sapropel example presented in §3b is new and preliminary results are shown here for the first time. We briefly present in §4 another type of geochemical rhythmic pattern: oscillatory zoning (OZ). Concluding remarks are given in §5.
2. Liesegang patterns: general considerations
‘Liesegang bands’ are a type of rhythmic pattern characterized by a macroscopic (millimetre-to centimetre-scale) spatial arrangement of precipitate aggregates, each made up of a large number of crystallites of various sizes. These chemical oscillations result from the coupling between the diffusion of a dissolved reactant and the nucleation, growth and/or ripening of the precipitate product. At the end of the nineteenth century, Liesegang  investigated the striking spatial patterns obtained when two counter-diffusing solutions in a gel substrate produce a precipitate in the wake of a moving reaction front. The resulting pattern consisted of regularly spaced precipitate bands (in one dimension) or rings (in two dimensions) . Three-dimensional patterns may also occur in the form of helicoidal or chaotic structures. These types of ‘periodic precipitation patterns’ were first mentioned by F. F. Runge in 1855 and Ostwald proposed in 1897  to name them ‘Liesegang bands’ (or ‘rings’). These patterns are easily produced in the laboratory and also occur in engineering applications and in biological or geochemical systems. It is the latter category that concerns us here. Sadek & Sultan  review a range of examples of Liesegang bands emphasizing their ubiquity and the richness of their dynamical features.
The formation of Liesegang bands can be understood in terms of the coupling between the diffusion of dissolved reactants and the kinetics of precipitate formation (nucleation, growth and ripening). To fix the ideas, assume that a dissolved species A (outer species) moves into a reservoir containing species B (inner species) to produce a colloidal precursor C to the precipitate . In simple systems, empirical laws were found to rule the spatial arrangement of the Liesegang bands and their formation timing: (i) the Jablczynski law establishes that the relationship between the position of the nth band and that of the previous one is linear (for n sufficiently large), (ii) the Matalon–Packter law indicates that the slope of this linear relation is itself a linear function of the inverse of the outer species concentration, and (iii) the time scaling law says that the formation time of the nth band is proportional to the square of its position (again for n not too small), suggesting that diffusion plays an important role. However, in a natural system exhibiting Liesegang bands, the spatial pattern does not necessarily exhibit these simple scaling laws, as the system may be too complex or its spatial extent too short.
To understand how self-organized precipitation bands are generated, Ostwald  proposed a mechanism that would be called in modern parlance a ‘pre-nucleation’ model. In this category of models, nucleation and growth of precipitate occur when the reactant activity product exceeds some threshold (critical supersaturation for nucleation or solubility product for growth). Assume that a precipitate band forms at some distance from the reactor edge. Then, the concentration of reactant locally decreases, the activity product decreases and, if diffusion is not too fast, growth ceases. Meanwhile, the reaction front moves further away until the activity product reaches again its threshold and another band is created. The cycle proceeds in this way but with modification from diffusion, thus producing a sequence of bands with progressively larger widths. This sequence of events can be mathematically described by using a simple macroscopic phenomenological approach [5,8], cellular automata , critical dynamics models [10,11] or models based on classical nucleation theory and growth kinetics [12–14]. These models are consistent with the Jablczynski law, the Matalon–Packter law and the time scaling law. In particular, Dee’s model  can be reformulated as a set of seven coupled nonlinear partial differential equations (one each for the concentration of species A, B, C and P plus the number density of crystallites, their average size and area) and has been found to be quite useful [14–16].
Post-nucleation models assume that the nucleation phase is terminated and couple diffusional transport, growth and Ostwald ripening. The latter phenomenon is owing to the fact that the concentration of a component in equilibrium with a crystallite depends on the size of the crystallite, so that smaller crystals are more easily soluble than larger ones. As such, post-nucleation models provide an explanation for the formation of banded patterns in an otherwise homogeneous system. It also gives a more complete description of the time evolution of the pattern, including the creation of satellite doublets within each band. Homogeneous concentration profiles are inherently unstable in systems described by post-nucleation mechanisms. Indeed, imagine that the size of a crystallite at a certain site increases slightly over its steady-state value, owing to environmental fluctuations. Then, the local equilibrium reactant concentration decreases slightly, and the supersaturation increases, leading to an increased growth and local depletion of the reactant. A diffusional flux is then generated from the neighbouring sites, which cause the reactant concentration to decrease and the precipitate to dissolve at these sites. Eventually, the reactant concentration increases as a result of dissolution and reactant is transported away from the original perturbation by diffusion, enhancing precipitate growth further down. The initial perturbation is thus seen to cause a sequence of bands with enhanced precipitate concentration, separated by zones with little (or no) precipitate. A simple mathematical model that encapsulates these phenomena is provided by the competitive growth model . This model is compatible with the Jablczynski law, the Matalon–Packter law and the time scaling law. It also leads to the generation of satellite doublets: as time progresses, each band evolves to a pattern characterized by a high precipitate concentration on the wings and a minimum in the middle. This characteristic pattern is due to the high flux of dissolved inter-band material that favours growth on both margins of the band, followed by intra-band ripening. Models that combine both pre-nucleation and post-nucleation mechanisms are also available [15,16,18].
3. Examples of Liesegang patterns in geochemical systems
Without being an exhaustive list, we describe here many rhythmic growth patterns seen in a wide range of geological settings (igneous, sedimentary, hydrothermal and metamorphic), in which a Liesegang periodic precipitation mechanism operates. This review emphasizes the richness and variety of Liesegang patterns in geological settings.
(a) Liesegang band patterns in igneous systems
(i) Layered intrusions
The first example is seen in many layered intrusions. These systems are generated by the slow crystallization of a magmatic melt intruding a system of older rocks. The term ‘layered’ is associated with large-scale variations in the textural and mineralogical characteristics of the resulting rock. These systems often constitute mining sources of platinum, copper, chromium and nickel. Examples are found in the Stillwater complex (Montana), the Bushveld complex (South Africa) and the Skaergaard intrusion (Greenland). These layered intrusions often exhibit  complex large-scale banded series characterized by alternations in rock type. For example, the Stillwater intrusion shows layers of norite, gabbro and anorthosite, presumably formed by external variations in the characteristics of the parent magma. However, a fine-scale layering pattern is also observed in the norite layer of the Stillwater complex . It is characterized by a series of fine dark pyroxene doublets segregated from the plagioclase matrix (separated by a few millimetres). In the so-called Ultramafic Zone of the Stillwater complex, another fine-scale layering is also observed [21,22]. It is characterized by thin dark layers of chromite (FeCr2O4) and olivine (Mg2SiO4) grains alternating with thicker light-coloured layers composed of olivine alone over a scale of a few millimetres. Similar layering is also seen in other layered intrusions, such as the Skaergaard [23,24] where spinel (MgAl2O4) bands alternate over a scale of 10 cm with bands of silicate minerals.
A post-nucleation mathematical model has been successfully proposed [20,22,24] to explain these band patterns. In the simplest version , the dynamics was described by the one-dimensional reaction–diffusion equation for the reactant concentration in the melt. The growth rate was assumed to be first order in the supersaturation, measured with respect to an equilibrium reactant concentration that follows the Gibbs–Thomson relation. The initial condition corresponded to a uniform sol subjected to a small local perturbation. These equations suffice to generate a banded precipitate pattern consistent with the layering observed in the Stillwater complex, including the formation of ‘doublets’ (sharply defined associated pairs of precipitate bands).
The more complex banded pattern seen in the Ultramafic Zone has also been successfully modellized  by a version of the post-nucleation mechanism in which two minerals grow from a two-component melt. Heat conduction and latent heat release were also taken into account. Finally, the decimetre-scale layering of the Skaergaard has been modellized  using a simplified version of this latter model, in which temperature was maintained constant, but where strong advective mass transport (resulting from compaction or magma flow) dominated over diffusion. Again, a post-nucleation mechanism appears responsible for the generation of doublets.
Other models [25,26] have also shown that thermal effects can generate centimetre-scale oscillations in igneous systems, such as those observed in the Stillwater complex. There, a contact between a cool host rock and a hot magma is considered. Thermal conduction is coupled with nucleation and growth of a single-component crystal in the magmatic melt and includes latent heat release, but not ripening. In some conditions, damped thermal oscillations and oscillations in the crystal size are obtained as a function of the distance from the contact point. This type of model has been recently generalized  to the crystallization of a binary eutectic melt. Complex undamped oscillations in the crystal characteristics are also obtained in this model which can be applied to layered intrusions.
(ii) Orbicular rocks
Orbicular granitic rocks provide another spectacular, albeit rare, occurrence of banded patterns in igneous systems . The orbicules consist of 10 cm spheroidal structures in which a central core is surrounded by one or many concentric shells of light and dark granitic material. The orbicules are embedded in a plutonic rock matrix, which is petrologically distinct. The core is typically a xenolith of plutonic or metamorphic origin, whereas the shells consist of radially or tangentially oriented feldspar crystals (plagioclase and/or K-feldspar) mixed with quartz and mafic minerals (biotite, amphibole) in the dark shells. The feldspar crystals often exhibit a feather-like (‘plumose’) habit, indicative of a fast crystal growth in the radial direction, from the core outwards. Another indication of their diffusion-controlled moderately fast growth rate is the reverse zoning exhibited by the plagioclase crystals : in the radial direction of the plagioclase crystals, the anorthite (CaAl2Si2O8) content tends to increase systematically. The orbicule origin is still poorly understood, but studies of specific outcrops (Fisher Lake, California ; Ploumanac’h, France ; or Sierra de Velasco, Argentina ) point to a complex nucleation/growth process that may involve Liesegang band formation in the generation of multi-shell orbicules. The basic mechanism is as follows. Consider the intrusion of a pocket of ‘orbicular’ magma in a preexisting slightly cooled granitic mush. If the intrusive magma is hotter, the contact zone will become superheated. Similarly, if the intrusion has higher water content, the liquidus temperature is lowered  and an effective superheating ensues. In both cases, most of the small crystals in the preexisting mush melt and only the largest ones survive. The next step involves the formation of the feldspar shell from these cores. After a quick cooling phase, or if water gets lost as a result of fracturing or volcanic eruptions, the mush could become moderately strongly undercooled. The homogeneous nucleation rate is negligibly small. On the other hand, the growth rate of feldspar crystals that have heterogeneously nucleated on the core crystals is large. Plumose habit and reverse zoning (for plagioclase) result. An orbicule could thus form in a matter of a few weeks! However, the presence of many light and dark shells in an orbicule is still a challenge to explain. A first step in a mathematical description of the generation of dark and light shells is provided in Wang & Merino , where plagioclase and biotite are seen to crystallize from a silicate mush in an alternating fashion when the model parameter values fulfil some conditions. But more work is needed.
(b) Liesegang band patterns in sedimentary systems
(i) Ferruginous sandstones and claystones
Classic examples of Liesegang bands in sedimentary systems are seen in ferruginous sandstone deposits. For instance, embedded in sandstones from the Middle Jurassic Hojedk formation (Iran), there are centimetre-scale structures that exhibit a variety of ring patterns , whereby dark concentric bands alternate with light-coloured material. Mineralogically, the light bands are grains of quartz, feldspar and ferromagnesian minerals, mica and various rock fragments, with some iron oxide and hydroxide minerals precipitated between the grains. The dark bands are enriched in iron oxide and hydroxide minerals that act as cementing material as well as substituting for the grain minerals. The formation mechanism proposed involves weathering of the sandstones with iron-rich meteoric waters infiltrated through fractures. This configuration generates two counter-diffusing solutions at the weathering front: one inward solution enriched in iron and iron oxides and one outward flow enriched in aluminium and silica. Although no detailed quantitative model is presented, this configuration has the potential to generate classic Liesegang bands. A similar mechanism has been proposed  to explain the formation of haematite (Fe2O3)-enriched dark bands in dolomite (CaMg(CO3)2) from the Ordovician Upper Arbuckle Group (Oklahoma) as a result of the infiltration of fluids enriched in iron and calcium through fractures in the rock.
Liesegang bands have also been reported  in an Oligocene–Miocene claystone formation from the North Alpine Foreland basin (Germany). Here, a 10 cm thick layer shows symmetric bands of coarse quartz crystals alternating with zones of finer quartz grains. The authors proposed that the pattern is the result of homogeneous nucleation and growth of quartz from a silica-rich gel-like organic mud deposited in an ancient anoxic lake. After partial burial (a few metres), the sediment was then presumably exposed to oxygen-rich groundwaters coming from the top and bottom directions. The counter-diffusion of oxygen in the sediment then led to the formation of a symmetric quartz precipitation pattern.
(ii) Ruin marble
The so-called ‘ruin marble’ is an instance of banded porous calcareous rocks used for decorative purposes. The characteristic texture of ruin marble consists of fractures separating the rock into independent Liesegang systems in which periodic precipitation of (typically) ferric oxy-hydroxides occurs. Marko et al.  have studied Cretaceous–Palaeogene clay-limestone from the western Carpathians (Slovakia) exhibiting a ruin marble texture. They reconstruct the history of the deposit as follows. After deposition of a calcareous mud, burial diagenesis incorporated a ferrous component in the rock. As a result of fluid overpressure build-up, a network of fissures was created, which were then filled in by calcite. Uplift brought the formation closer to the surface, where it disintegrated into boulders and was subjected to the infiltration of oxidizing meteoric waters. Capillary forces encouraged infiltration of these waters along the healed fissures creating periodic precipitation bands of iron oxy-hydroxides throughout the boulders. The scenario is consistent with the observations but no quantitative model is presented. McBride  also described Liesegang patterns in ruin marble from the Late Paleozoic Dimple Formation (Texas) and from Pakistan. These are characterized by millimetre to centimetre-scale bands of iron oxides (mostly goethite FeOOH) in a calci-siltite. The formation mechanism appears similar to the marble studied by Marko et al., although the details of where the iron originated are not known.
(iii) Pyrite bands in sapropels
We present now an early diagenetic system as a final example of a sedimentary system exhibiting Liesegang bands: Holocene sapropel sediments from the eastern Mediterranean. Sapropels are defined as sediments characterized by an alternation of decimetre scale organic matter-rich layers (up to 30 wt%) and organic matter-poor sediment . The alternation is believed to be linked to changes in the palaeoclimatic environment [40–43]. In anoxic conditions, bacteriological decomposition of organic matter consumes iron hydroxide and/or sulfate ions, thus enriching the sediment in hydrogen sulfide and dissolved ferrous iron, which in turn can act as reactants for pyrite (FeS2) precipitation. The fact that organic matter-rich layers are spatially separated suggests that sulfide and iron counter-diffusing flows are possible between successive sapropel layers, and therefore that Liesegang pyrite bands may exist. Indeed, such bands have been detected by Passier et al.  in a sapropel system at more than 3 km depth. The pattern is characterized by five pyrite bands occurring just below the most recent sapropel layer, between 25 and 33 cm below the sediment–water surface.
We have recently proposed  a one-dimensional early diagenetic model that considers the main biogeochemical reactions occurring in an anoxic sediment and incorporates a pre-nucleation mechanism for the precipitation of pyrite:
(a) Microbial decomposition of organic matter (CH2O) by reduction of sulfate and iron hydroxides (ferrihydrite Fe(OH)3) and by methanogenetic fermentation: The rates are assumed to follow Monod kinetics [46,47] with respect to the electron donor and first order with respect to organic matter.
(b) Secondary redox reactions: The appropriate rates are assumed to be proportional to the concentration of the electron donor and to the concentration of the electron acceptor.
(c) Equilibrium reactions for adsorption of iron and sulfide dissociation: Here, (S−H), (S−Fe+) symbolize free and adsorbed iron surface sites, respectively. A simple constant adsorption isotherm was chosen.
(d) Precipitation/dissolution of pyrite and its colloidal precursor (FeS2)col via iron monosulfide formation:
Note that the hydrogen-forming pyrite pathway dominates and was the only one considered in the early version of our model . We are now adding the zero-valent sulfur pathway to our reaction network. The precipitation rate of iron monosulfide is first order in the supersaturation (where [ ] denotes concentration, KFeS is the solubility constant and aH the activity of hydrogen ions, assumed constant), whereas its dissolution rate is assumed proportional to the undersaturation −ΩFeS and to the concentration of FeS. The kinetics of pyrite precipitation from the hydrogen-forming pathway assumes a steady-state elimination of the intermediate (FeS)col according to the method described in the study of L’Heureux . The result reads  as follows: 3.1where sFeS2=[FeS][H2S]/KFeS2−1 is the supersaturation of pyrite with a solubility constant KFeS2, G is a kinetic coefficient, N is the number density of pyrite crystals per dry sediment mass, v is the molar volume of pyrite, J is the nucleation rate (number of new crystals created per unit time per pore water volume) and r* the critical nucleation radius. Also, F=ρsφs/φ is a concentration unit conversion factor that involves the mean sediment density ρs, the porosity of the sediment φ and the solid volume fraction φs=1−φ. The model adopted here applies a pre-nucleation mechanism. The first factor in equation (3.1) describes nucleation of new pyrite crystals, whereas the second term corresponds to crystal growth proportional to the surface area of pyrite. Finally, the dissolution of pyrite adopts an expression similar to the growth term.
Our reactive-transport model is based on the continuity equation for the mole numbers of each chemical species. For dissolved and solid species, respectively, this is 3.2and 3.3Here, x is depth measured downward from the sediment–water interface, Ci is the concentration (mole per pore water volume) of dissolved species and Cs the concentration (mole per mass of dry sediment) of solid species. The net reaction rates are Ri and Rs for dissolved (mole per pore water volume per unit time) and solid species (mole per sediment mass per unit time), respectively. Di is the porosity-dependent diffusion coefficient and Db is the bioturbation coefficient. Finally, U and V are the advection velocity of the pore water and solid sediment matrix, respectively, measured from a reference fixed to the sediment–water interface. Here, compaction is ignored, so that U=V corresponds to a constant sedimentation rate. The crystal number density obeys a similar continuity equation 3.4In our model , the porewater sulfate concentration and the dissolved iron concentrations were fixed at the boundaries of the system. Pyrite Liesegang bands were obtained for intermediate values of these concentrations. Other parameters, such as the pH, organic matter concentration at the top of the sapropel, the organic matter decomposition rate and the size of the system, also play a role in determining the number of pyrite bands and their width and position.
We are currently generalizing the situation to a more realistic situation where the boundary conditions are time-dependent in order to mimic a plausible depositional history of the top two sapropel layers over the last 12 000 years. First, a first sapropel layer with enhanced organic matter flux is deposited. This is followed by the deposition of an organic matter-poor layer enriched in detrital iron Fe(OH)3. A layer of organic matter-poor material (‘normal’ sediment) is then deposited, followed by the top sapropel layer. Finally, a cap of normal sediment closes the system at the top. The details will be published elsewhere. Nevertheless, figure 1 illustrates a numerical solution using reasonable parameter values. It is seen that Liesegang bands similar to the observed pattern are effectively generated. This suggests that the presence and modelling of cyclic precipitation patterns in an early diagenetic system can be used as a tool to reconstruct the depositional history of that sediment system and give clues as to what the palaeo-environmental conditions could have been.
(c) Liesegang band patterns in metamorphic systems
Liesegang precipitation patterns in metamorphic systems are rather uncommon. We discuss here two examples involving contact metamorphism.
(i) Liesegang bands in skarns
Skarns are formed when a granitic magma intrudes a carbonate host rock. Magmatic hot fluids may bring about chemical alterations of the host rock in the contact area (metasomatism). In an example of a skarn system from the Ocna de Fier-Dognecea ore deposit (Romania, Upper Cretacian), a very wide range of macroscopic banded textures were reported : millimetre to centimetre-scale bands of magnetite (Fe3O4) alternating with calcite (CaCO3) or diopside (MgCaSi2O6) or garnet (some of which exhibiting branching of bands and magnetite dendrites), nodules (garnet core surrounded by rings of garnet and magnetite), rings of galena and sphalerite (ZnS) in orbicules, mottled patterns (magnetite needle-like crystals aligned in a calcite matrix) and mossy patterns showing skeletal magnetite crystallites in marble. Ciobanu argues that these patterns were formed in the initial (prograde) phase of the deposit formation when a hydrothermal fluid issued from a magmatic source and enriched in copper, iron, lead and zinc infiltrated a limestone bed. The banded and mossy patterns would result from a pre-nucleation mechanism (with some influence from the flow and from Ostwald ripening), whereas the orbicular and nodular patterns would mostly be manifestations of Ostwald ripening (post-nucleation mechanism). However, no quantitative modelling is offered.
(ii) Metamorphosed chert nodules
Another instance of contact metamorphism has generated Liesegang-like bands: Holness  observed millimetre-scale banded patterns of forsterite olivine (Mg2SiO4) alternating with calcite (CaCO3) at the rim of chert nodules in the Kilchrist intrusive complex (Scotland, Palaeocene). Cherts are fine-grained silica-rich sedimentary rocks and the nodules were formed when a magmatic intrusion chemically altered a diopside (MgCaSi2O6) seed in a dolomite (CaMg(CO3)2)-rich country rock. The relevant chemical reaction is: 3 dolomite+diopside→4 calcite+2 forsterite+2 CO2. Metasomatism thus produced a calcite–olivine rim around the diopside nodule. It appears that the size distribution of olivine grains is consistent with a post-nucleation mechanism. Macroscopically, mineralogical segregation owing to Ostwald ripening is enhanced at the point of contact between the host rock and the magmatic intrusion, where the fluid is hotter. Thus, the position of the banding pattern in the complex can give useful information as to where the intrusive contact was located.
3.4 Liesegang patterns in ore deposits: Mississippi Valley-type sphalerite
Mississippi Valley-type (MVT) deposits are lead and zinc sulfide ore minerals (galena (PbS), sphalerite (ZnS)) in association with sedimentary carbonate host rocks. Their name comes from the Mississippi Valley area in central USA (Missouri, Iowa, Illinois, Wisconsin) but they are actually found all over the world [50,51]. They constitute important economic sources of lead and zinc and other metals. MVT sphalerite often exhibits a striking pattern of colourful (white, yellow, reddish brown and black) sub-millimetre clusters of bands radiating from a centre (figure 2). These bands are made up of a large number of ZnS acicular crystallites oriented normally to their boundary. Band doublets are also sometimes observed . Also, galena crystals often exhibit a dendritic habit, suggesting far-from-equilibrium growth conditions. Using scanning electron microscopy and electron microprobe analysis, we have shown  that the band colour correlates with iron substituting for zinc in the sphalerite: in a sample from the Pine Point (Northwest Territories, Canada) deposit, light bands have an iron content of 3 mol%, whereas dark bands have 10 mol% iron.
It is believed [50,54,55] that MVT deposits result from the mixing of two hot brines within a porous carbonate rock: one brine containing dissolved ferrous iron and zinc chloride complexes, and another enriched in sulfides. The corresponding reactions are , ZnCl2+H2S→k2ZnSaq+2H++2Cl− and (1−p)ZnSaq+pFeSaq→Zn1−pFepS, where p is the iron molar composition. Using this reaction scheme, we have modelled [52,56,57] the band formation in MVT sphalerite with Liesegang post-nucleation mechanisms. The concentration (mole per pore water volume) of the relevant species Fe2+, ZnCl2, FeS and ZnS evolved according to a reaction–diffusion equation combined with the radial growth velocity 3.5where βi are kinetic coefficients and [FeS]eq, [ZnS]eq are the concentrations of FeS and ZnS in brines that would be in equilibrium with a crystallite of radius R. They depend on R according to the Gibbs–Thomson relation. The Liesegang pattern is defined by the spatial variation of the composition p. In our model [52,57], p was an explicit dynamical variable ruled by 3.6where L≪R is a length scale that parametrizes the roughness of the crystallite surface . In the steady state, p is then defined by the relative growth rate of FeS and ZnS. The model consists then of six differential equations for the concentration of Fe2+, ZnCl2, FeSaq and ZnSaq, the crystallite size R and its iron molar composition p. In Katsev & L’Heureux , we have solved such a model in a two-dimensional system with cylindrical symmetry, in which the metal brine is injected (with fixed concentration) at the core of the cylinder. Such MVT pipes are actually observed at the Pine Point deposit. Homogeneous initial conditions were chosen for the crystallite size. The results indeed show the existence of bands of iron molar composition varying between 1 and 10 mol% as a function of distance from the core. Formation of doublets is also possible. However, the inter-band regions are crystallite-free in the model. If the initial crystallite radius is chosen randomly (over a uniform distribution) instead, then a continuous band pattern is obtained. A more realistic model should also consider the nucleation of new crystallites and competition for space with the growth of other minerals, such as galena.
4. Oscillatory zoning patterns
Another important class of rhythmic geochemical patterns is found in the so-called OZ. This type of pattern is characterized by more or less periodic fluctuations in the concentration of a chemical species within a single crystal, along a core-to-rim profile. The fluctuations typically occur over a length scale of the order of sub-micrometres to tens of micrometres and constitute therefore a microscopic pattern, in contrast to macroscopic Liesegang patterns. Perhaps, the most well-known mineral exhibiting OZ is plagioclase feldspar (a solid solution series between a calcium-rich end member (anorthite) and a sodium-rich one (albite)). For some time, it was thought that OZ occurred rarely in other minerals, but the review by Shore & Fowler  makes it clear that OZ is actually ubiquitous. It can be found in all mineral classes in a variety of geological environments. Typically, it occurs in solid solutions (e.g. plagioclase, grossular-andradite garnets [59,60]) or in trace components adsorbed on the mineral surface as it grows (e.g. Mn or Sr in calcite [61,62]). For lack of space, we concentrate on the case of solid solutions here.
Two broad classes of mechanisms are proposed in the study of Shore & Fowler  to explain the formation of OZ. In ‘extrinsic’ mechanisms, the chemical composition fluctuations result from variations in the external conditions defining the growth environment: pressure, temperature, pH, fugacities, water content in silicate melts, etc. In ‘intrinsic’ mechanisms, OZ is the result of self-organization in a nonlinear growing environment. Of course, both mechanisms could occur simultaneously, generating fluctuating patterns with a rich spectral content.
(a) Oscillatory zoning in plagioclase
OZ in plagioclase is commonly found in volcanic rocks, such as andesite and diorite (with anorthite content smaller than 50 mol%). It is characterized [63,64] by variations in the anorthite content of typically 5–15% over a scale of tens to hundreds of micrometres. By plotting the thickness of one zone as a function of the previous one (first-return map), Higman & Pearce  suggested that an intrinsic dynamics is at least partly responsible for the formation of OZ, without precluding the influence of systematic or random changes in the external growth conditions.
In Brandeis et al. , a model including nucleation from a magma melt, crystal growth, heat transport by conduction and latent heat generation was proposed in which damped oscillations are possible, but on a centimetre-scale zoning. In Haase et al.  and Ortoleva , an isothermal diffusion model of plagioclase growth is proposed with an ad hoc autocatalytic reaction scheme, which is difficult to interpret physically. In Allègre et al. , another isothermal diffusion model is proposed in which the growth rate is delayed with respect to the value expected from kinetic theory. There, transient OZ is possible, but the delay is also difficult to interpret. Lasaga  proposed a diffusion model with a value of the growth rate that is expected from kinetic theory. No OZ was found in the model. Wang & Wu  proposed a model for which the concentration of a component in the crystal is proportional to its growth rate. The possibility of OZ arose as the system has the potential to switch between two growth regimes, but no detailed dynamics is presented and diffusion is ignored. In Tsune & Toramaru , OZ is obtained by assuming that the system is made to switch between two growth regimes (onto a soft surface and a rough surface). Notwithstanding the variety of OZ formation models that are available, much work still needs to be done in order to obtain a clear picture of OZ formation in plagioclase.
OZ in plagioclase occurs often in far-from-equilibrium conditions (large undercooling) . L’Heureux and co-workers [71–73] have proposed an isothermal nonlinear model of OZ in plagioclase for such a situation. The model is based on the diffusive transport of anorthite component in the silicate melt with a nonlinear boundary condition at the growing crystal interface (since the growth velocity is a strongly increasing function of the melt concentration). A nonlinear partition coefficient K(co)=cs/co, which relates the anorthite molar composition in the crystal cs to its value co in the melt in the vicinity of the crystal, is also used. K is a function of co and of a single partition parameter KD. The diffusion equation has only one steady state (fixed point) and its linear stability analysis can be performed. In near-equilibrium conditions, the partition function K (and KD) is larger than one, and it can be shown that no oscillations are found (the fixed point being stable), in accordance with Lasaga’s result . However, in far-from-equilibrium conditions, it is possible for KD to be smaller than one (see the discussion below). Then, K is also smaller than one and a more interesting dynamical behaviour is possible: the fixed point can undergo a Hopf bifurcation, whereby it makes a transition from a stable focus (where the system approaches it with underlying damped oscillations) to an unstable one. Numerical solutions confirm that an autonomous periodic solution appears just beyond the bifurcation point. As KD gets smaller, the system goes through a rapid sequence of period doublings. It is interesting to note that two frequencies were detected in the power spectrum of observed OZ sequences in plagioclase . The fact that one frequency is twice the other one suggests that period doubling could be detected in such a system, in agreement with the model prediction. For still smaller KD, the system becomes chaotic. The first-return map in the chaotic regime shows features that are similar to the one obtained from observations .
But is it physically possible for KD to be smaller than unity? In fact, in far-from-equilibrium conditions, one can use a molar composition cs of anorthite in the crystal that is determined by kinetics rather than thermodynamics. Indeed, Wang & Wu  assumed that cs is proportional to the growth rate of the anorthite component. They obtained an expression relating cs to co that involves the growth kinetic coefficient and the anorthite–albite interchange energy. This formalism provides a microscopic expression for KD. It turns out that, for cs<0.5, KD can be made smaller than one when the absolute value of the exchange energy is sufficiently large or when the rate coefficient for the growth of albite components is much larger than that for anorthite. In fact, if the growth rate of albite is much larger than that of anorthite, albite is preferentially incorporated in the plagioclase crystal (small cs), leaving a high anorthite composition in the residual melt (large co), both factors contributing to make K and KD smaller than one.
The physical mechanism for the generation of OZ in the model is then as follows. For KD>1, growth consumes anorthite from the melt, so that the anorthite composition co in the melt in the vicinity of the crystal front is smaller than its bulk value. Assume that a perturbation increases co. The growth velocity will therefore increase, which results in the consumption of more anorthite and a decrease in the perturbation (stabilizing effect). However, for KD<1, growth results in a larger anorthite composition at the crystal front (owing to the faster growth rate of albite components, for instance) and the perturbation gets amplified in a positive feedback effect. Of course, diffusion will eventually reduce the concentration gradient, so that co will evolve towards its steady-state value. The idea is that the fixed point loses its stability, so that the potential to generate cyclical variations exists when KD becomes smaller than unity. This is basically the same mechanism as the one discussed by Wang & Merino .
(b) Oscillatory zoning in an experimental system
In spite of its ubiquity in nature, OZ is very difficult to obtain in the laboratory under controlled conditions. One notable exception is the growth of barite (BaSO4)–celestite (SrSO4) crystals from aqueous solutions at room temperature by Putnis’s group [75–78]. In their experimental set-up, a 28 cm long gel column linked one reservoir of Na2SO4 solution to another reservoir containing a mixture of BaCl2 and SrCl2 solutions. As a result of the counter-diffusion of these mother solutions in the column, precipitation of barite–celestite crystals occurred. After a time of about one month, the crystals were analysed by back-scattered electron imaging and their chemistry determined by scanning electron microscopy. Depending on the concentrations of the mother solutions, various zoning patterns were obtained within the single crystals. In particular, for larger concentrations of the mother solutions, OZ was obtained, where fluctuations in barium content occurred along a core-to-rim profile over a length scale of 1–10 μm. This important experiment then shows that OZ is possible when external conditions are maintained fixed (apart from possible small random environmental fluctuations), thus suggesting that self-organization (intrinsic mechanism) plays an important role in explaining the formation of this pattern.
L’Heureux & Jamtveit  proposed a nonlinear model of OZ in the barite–celestite solid solution system that exhibits many features consistent with the experimental facts. The basic formation mechanism can be described as follows. Consider two end-members A and B of a binary solid solution series in solution in proximity to a growing AB crystal. When a dissolved A unit, say, gets close to a kink site on the growing surface, we assume that it has a higher probability of being integrated in the crystal if that kink site is locally enriched in A units (this is expected to be true on energetic grounds). Thus, an A-rich crystal surface will favour the growth of more A units. However, diffusion limits the transport of A units in the solution. As the solution in close proximity to the crystal gets deleted in A units, B units are more likely to grow and the system switches from an A-mode growth to a B-mode growth. Eventually, B units become deleted in turn and freshly imported A units grow preferentially on the surface, and the cycle repeats itself.
The nonlinear model developed in the study of L’Heureux & Jamtveit  reduces the dynamics of the system to a set of four coupled nonlinear ordinary differential equations for the SO42−, Ba2+ and Sr2+ concentrations in the gel and for the molar composition of BaSO4 in the crystal. Besides diffusion, the essential ingredient in the model is the implementation of an autocatalytic (nonlinear) factor in the crystal growth velocity, based on the idea that the probabilities of finding a kink site on a two-dimensional surface are proportional to the square of the inverse average distance between kinks , which is itself a linear function of the molar fraction. There always exists at least one steady state (fixed point) for the dynamical system. The linear stability analysis of this system  suggests a rich self-organized behaviour: for some values of the parameters, the steady state undergoes a Hopf bifurcation. Numerical solutions confirm that a periodic attractor (limit cycle) appears as the fixed point loses its stability, which is compatible with the occurrence of self-organized OZ. Moreover, there are parameter values for which three fixed points coexist, one of them unstable and the two others stable: one Ba-rich and the other Sr-rich. This is an example of bistability and is a feature of many nonlinear dynamical systems: according to the initial conditions, the system will evolve towards one or the other stable fixed points. These findings were confirmed by a direct numerical solution of the original partial differential equations . We have also shown  that the system is particularly sensitive to concentration random fluctuations in the neighbourhood of the bistable region: in this region of parameter space, small random fluctuations are sufficient to induce large oscillatory excursions in the crystal composition profile.
In the model described previously, the growth mechanism includes only incorporation of material at kink sites. A more realistic growth law has also been proposed by Kalischewski et al. , where a similar approach is used except that the growth law includes species adsorption on the crystal surface (with appropriate binding energy parameters), surface diffusion and incorporation or desorption. The boundary conditions for the concentration fields far from the crystal were also different (Neuman rather than Dirichlet). In that model, the steady state can also undergo a Hopf bifurcation, thus generating OZ. However, the steady state is unique and no bistability is found.
In this review, we have presented a variety of geochemical rhythmic patterns and some of the modelling strategies that are used to understand their formation. We have concentrated on two broad classes of patterns in a variety of geological settings: a macroscopic pattern defined by an aggregate of crystallites (Liesegang bands) and a microscopic pattern characterized by chemical composition fluctuations within a single crystal (OZ).
There are many other examples of self-organized structures in geochemical systems, but space limitations do not allow us to review them. Let us just mention two interesting examples of geochemical oscillators: agates and banded iron formation (BIF). For a long time, bandings in agates were believed to be due to a classic case of Liesegang periodic precipitation of silica resulting from the filling up of basaltic cavities with hydrothermal fluids. However, it was shown  that patterns in agates are rather due to a self-organized geochemical oscillator in a closed system: an autocatalytic reaction that transforms in situa chunk of hydrous silica gel into amorphous silica. Another recent example of a geochemical oscillator can be found in the famous Archean/Early Proterozoic BIFs, which are deposits of chert alternating with iron-rich minerals formed in the early ocean. It was shown  that BIF may have been formed by the autocatalytic oxidation of ferrous iron in the presence of the complex leached out from low-Al oceanic crust rocks (such as komatiites). When reasonable parameter values are used, periodic generation of iron hydroxides and silica (chert) results from such a model.
We hope that we have convinced the reader that, in order to understand the origin of geochemical rhythmic patterns exhibited by a variety of rocks and minerals, the tools used in the study of nonlinear dynamical systems and the concepts of self-organization are essential. Thus, not only do these examples provide beautiful illustrations of the manifestation of nonlinear dynamics in nature, but they allow Earth scientists to better understand the genesis conditions under which these patterns are formed. In this manner, Earth scientists can obtain more information on the past history of the Earth and, perhaps, better understand its future evolution.
6. Funding Statements
I wish to thank the Natural Sciences and Engineering Research Council of Canada for financial support.
I wish to acknowledge the help from many collaborators: Rimma Bektursunova, Sergei Katsev, Anthony D. Fowler and Bjørn Jamtveit.
One contribution of 13 to a Theme Issue ‘Pattern formation in the geosciences’.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.