## Abstract

A wide range of explanations has been proposed for the origin of repetitive layering in mafic-ultramafic and in (per)alkaline intrusions. Here we propose that the interaction of mineral grains that sink and float in the crystallizing magma is an alternative mechanism that can explain many of the features of layered intrusions, without the need to invoke extrinsic factors. Similar to traffic jams on a motorway, small perturbations in crystal density develop that impede further ascent or descent of buoyant or heavy minerals, respectively. These “traffic jams” separate layers of magma from the rest of the magma chamber. The magma in the individual layers further evolves as a largely independent subsystem, with gravitational sorting organizing the mineral distribution within each layer. Layering can develop in the intermediate range between full mineral separation in low-viscosity or slowly cooling magma chambers and homogeneous crystallization in high-viscosity or fast-cooling chambers. This self-organization mechanism provides a novel explanation for the formation of rhythmic layering in low-viscosity magmas, for example in the Ilímaussaq igneous complex in southwest Greenland.

## INTRODUCTION

Layered intrusions are characterized by layering on different scales, from >100 m down to the centimeter scale (Wager and Brown, 1967; Parsons, 1987). Although the (large-scale) layering can in some cases be explained by multiple magma intrusions (Irvine and Smith, 1967; Huppert and Sparks, 1980), layering may have developed from a single magma batch in many other cases, such as the nepheline syenites of the Ilímaussaq complex, southwest Greenland (Fig. 1A) (Larsen and Sørensen, 1987). How repetitive layering forms remains a controversial issue. Proposed mechanisms include gravity settling (McBirney and Noyes, 1979; Conrad and Naslund, 1989); the formation of double-diffusive layers (Chen and Turner, 1980); effects of compaction, filter pressing, and interstitial crystal growth (Sparks et al., 1985); alternation between density currents related to roof cooling and layers dominated by local nucleation (Morse, 1986); crystal escape from an unstable boundary layer in a convecting magma chamber (Sparks et al., 1993); degassing (Lipin, 1993; Pfaff et al., 2008); stratification of magma chambers with periodic density inversions (Tegner et al., 2006); repeated “self-stratification” in the mush zone (Nielsen and Bernstein, 2009); liquid immiscibility (Jakobsen et al., 2011); and mineral segregation during aging of mineral assemblages due to concentration gradients resulting from surface energy effects (Boudreau, 2011).

Here we propose a new model for layer formation, based on equations that are analogous, in a first-order approximation, to the inviscid Burgers’ equation (Burgers, 1974), which describes interesting pattern formation phenomena such as traffic jams (Helbing, 2001) and the formation of shock fronts (Kardar et al., 1986). Traffic jams form spontaneously on a motorway in heavy traffic, resulting from the interaction between cars, causing them to slow down where traffic is locally denser (Helbing, 2001). An analogous situation can arise in a cooling magma chamber: crystals form, and they will float or sink, depending on their density difference (Δρ) with the remaining magma. With increasing crystal concentration, the crystals start to impede each other’s movement, slowing them down (Marsh, 1981; Tomkins et al., 2005; Schmidt et al., 2012). This hindered settling causes small perturbations in crystal concentration to catch more crystals and to grow into “crystal traffic jams” or “mats” (Lauder, 1964), especially when different minerals have opposite buoyancies. These porous mats are barriers to further gravitational settling, but would allow melt percolation and redistribution during further compaction, as well as mineral re-equilibration (Boudreau and McBirney, 1997; McKenzie, 2011; Tanner et al., 2014). Here we present a simple numerical model to explore layer formation and patterns that can develop by this process, and compare the results with the layered nepheline syenites of the Ilímaussaq complex in Greenland (Fig. 1A) (Larsen and Sørensen, 1987).

### A FIRST-ORDER NUMERICAL MODEL

Assuming a sheet-like magma chamber with dominantly vertical temperature gradients, a one-dimensional model suffices to simulate the evolving crystal fractions in a vertical profile. The model spans 1500 m divided into 1500 elements. The magma chamber occupies the middle 500 m. Cooling and changes in mineral fractions were modeled in small time steps using a finite-difference scheme. Cooling was assumed to be purely conductive with a thermal diffusivity κ = 10^{−6} m^{2}/s for melt, crystals, and wall rock, starting at a normalized temperature *T*_{norm} = 1 (100% melt) inside the chamber and *T*_{norm} = 0 in the wall rock. *T*_{norm} = 0 was fixed at the top and bottom of the model.

The change in crystal fraction (∂*f*_{i}/∂*t*) of crystals of mineral *i* at height *z* in the magma chamber is a function of the gradient in fluxes (*v*_{i}*f*_{i}; *v* is settling velocity) and the crystallization rate, (∂*f*_{i}/∂*t*)_{cryst}:

Crystallization rate was set using (∂*f*/∂*t*)_{cryst} = –*c* · *f*_{melt} ·∂*T*_{norm}/∂*t*, where *c* is a mineral-dependent rate parameter. Two model minerals crystallized (both with *c* = 5), one with a density of ρ = 3000 kg/m^{3} and the other 2400 kg/m^{3}, from a melt with a density of 2700 kg/m^{3}. For a liquidus temperature of 950° (*T*_{norm} = 1) and a wall rock temperature of ∼120 °C (*T*_{norm} = 0) at the ∼3–4 km depth of the Ilímaussaq intrusion (Konnerup-Madsen and Rose-Hansen, 1984), setting *c*_{i} = 5 results in 95% crystallization at *T* = 710 °C.

The mean settling velocity (*v*) as a function of crystal density (*f*) is often described with:

*m*is an exponent that depends on crystal shape (Tomkins et al., 2005) and is usually assumed to be ∼4 (e.g., Schmidt et al., 2012).

*v*

_{0}is the Stokes settling velocity of an isolated single crystal (

*f*≈ 0), given by:

*g*the gravitational acceleration,

*r*the radius of the crystal, and η

_{0}the viscosity of the pure melt. Equation 2, however, does not describe the complex behavior at increasing crystal fractions, where long-range interactions between crystals come into play (Schwindinger, 1999; Suckale et al., 2012) until the crystal mush forms a solid framework at a few tens of percent of solid crystals (Marsh, 1981; Tegner et al., 2009). To approximate this behavior, we replaced Equation 2 with:

By choosing *n* = 2 and *f*_{max} = 0.7, Equation 4 is almost identical to Equation 2 up to *f* ≈ 0.35, after which the settling velocity rapidly drops (Fig. 2A). Because of the highly nonlinear nature of Equation 4 when *f*_{max} is approached, a low-pass filter was applied to suppress fluctuations less than 5 m in wavelength. Grains were assumed to be spherical (α = 2/9) with a constant radius of *r* = 1 mm. In this first-order model, we thus ignore variations in crystal size and shape, as well as variations in settling velocity due to, for example, clustering (Schwindinger, 1999).

## RESULTS

Keeping all other parameters constant, we varied η_{0}, the viscosity of the pure melt (Figs. 2B and 2C). At a very low η_{0} (≤10^{3} Pa·s), the heavy mineral accumulates at the bottom and the light one at the top. Concentrations in the melt in the middle never reach levels where barriers can develop, enabling full separation of the minerals. At high η_{0} (≥10^{5} Pa·s), the minerals are distributed homogeneously along the column: there is no separation of minerals at all. This is because crystals cannot sink or float any significant distance before the system locks up locally. Layering forms in the intermediate regime where crystals can sink or float a significant distance before lock-up, but not the full height of the magma chamber. A sensitivity analysis (Fig. 2D) shows that the results are robust: as long as settling velocities decrease with increasing crystal fraction (even using Equation 2), there is a viscosity range where barriers develop that define developing layers.

As the chamber cools inwards (Figs. 2C, 2E–2H), perturbations or waves in crystal concentration develop, and several of these self-organize into barriers that consist of heavy crystals lying on top of buoyant ones. These barriers define layers that consist of heavy crystals at the base and buoyant minerals at the top, trapped between barriers, and melt in between as long as the layer is not fully solidified. McBirney and Noyes (1979) argued that crystallization and differentiation of crystals in a magma could not involve crystal settling, because buoyant plagioclase is found in rhythmic layers at the floor of the Skaergaard intrusion on the east Greenland coast. However, our model shows that buoyant minerals can be trapped at the base of the liquid zone, as previously proposed by Lauder (1964) and Nielsen and Bernstein (2009).

The model predicts that dense minerals that crystalized early are concentrated at the base, and buoyant minerals at the top of each layer. This sequence is observed in the layered nepheline syenites (kakortokites) of the Ilímaussaq complex (Fig. 1), where each layer is dominated by dense cumulus amphibole (ρ ≈ 3400 kg/m^{3}) at the base, generally grading upwards into a zone rich in less-dense cumulus eudialyte group minerals (EGM) (ρ ≈ 2920 kg/m^{3}) and topped by a horizon with light cumulus alkali feldspar (ρ ≈ 2580 kg/m^{3}). Compositional data for EGM (Marks and Markl, in press) indicate that within each of the layers, the trapped magma chemically evolves as a quasi-closed system, largely independently from the rest of the magma chamber (Fig. 1B).

## DISCUSSION

It is well established that convection is expected to take place in magma chambers, especially at the low viscosities where layering or complete separation can occur (Marsh, 1989; Sparks et al., 1993). Convection is, however, not included in the numerical model that illustrates the formation of layers in the absence of other processes. Thermal convection arises from density instabilities as a result of the thermal expansion of the magma. At low crystal fractions, this effect is expected to dominate. Near the top and bottom of the chamber, where layering forms, variations in crystal load, with the concomitant variations in density and viscosity, will quickly override thermal instabilities and suppress convection here (Jaupart et al., 1984). It is thus expected that the magma chamber may be dominated by convection at the start of solidification and in the middle of the chamber, with layers forming at the top and bottom by the mechanism proposed here. Processes other than Stokes settling certainly play a role as well. Most notably, compaction of the layers with ensuing redistribution of melt, diffusion, and mineral re-equilibration (Boudreau and McBirney, 1997; Boudreau, 2011; McKenzie, 2011; Tanner et al., 2014) or the effect of latent heat (Morse, 1986) were not incorporated in the first-order model. It is, however, to be expected that these processes modify, but do not fully alter, the layering that results from the formation of barriers.

Above, the interaction between floating and sinking crystals was invoked to explain repetitive or rhythmic layering. Layered intrusions are, however, usually characterized by layering on multiple scales. For example, the exposed agpaitic part of the Ilímaussaq complex (Larsen and Sørensen, 1987) includes, from base to top, ∼300 m of rhythmically layered nepheline syenites (kakortokite) (Fig. 1A), 350 m of fine-grained and melanocratic nepheline syenites (lujavrite), and 600 m of sodalite-rich flotation cumulates (naujaite). A second set of simulations shows that a variety of patterns can develop with only three minerals, including large-scale, non-rhythmic stratification, as well as small-scale, rhythmic layering (Fig. 2I). Our simulations are admittedly oversimplified and ignore other processes that play a role in natural magma chambers. However, even in their simplicity they already develop rich self-organization patterns that are comparable to those found in nature.

### Requirements for the Formation of Layering

The main requirement for the formation of layers in a magma chamber is a relatively low viscosity and slow cooling rate to allow crystals to sink or float a significant distance. Because large intrusions cool slower than small ones, layered intrusions are usually large (hundreds of meters to kilometers in height) and (ultra)mafic (e.g., Skaergaard, Muskox in Canada, Rum in Scotland, Great Dyke in Zimbabwe) or peralkaline (e.g., Ilímaussaq, Khibina and Lovozero, both on the Kola Peninsula, Russia), as these melts have the lowest viscosity (Giordano et al., 2004). Cooling time (τ_{cool}) scales with the squared height (*H*) of a magma chamber: τ_{cool} ∝ *H*^{2}/κ (Marsh, 1989). The time to sink or float (τ_{Stokes}) a distance *H* mostly depends on melt viscosity (τ_{Stokes} ∝ *H*η_{0}/*gr*^{2}Δρ), as the other parameters in Equation 1 do not vary significantly. We define a dimensionless “separation number,” *S*, which is the ratio of these two times:

In our simulations of a magma chamber *H* = 500 m tall with two minerals, we found that layering occurred at η_{0} of ∼10^{4} Pa·s, giving *S* of the order of 100. If *S* << 100 (small chamber and/or high viscosity), there is no time for crystals to sink a significant distance before freezing, and a homogeneous distribution of minerals is expected. If *S* >>100 (large chamber and/or low viscosity), minerals can fully segregate to the top and bottom of the chamber. This indicates that layered intrusions on Earth, where the maximum height of magma chambers is on the order of 10 km, require melt viscosities below ∼10^{6} Pa·s. This would be different on other planetary bodies, where the gravitational acceleration and magma chamber sizes may be different. For example, Suckale et al. (2012) argued that the anorthositic lunar crust was formed by separation of buoyant plagioclase. Using their parameters results in an *S* value on the order of 10^{8}, well in the range of full separation of dense and buoyant minerals.

## CONCLUSIONS

Using a numerical model, we propose that the formation of layering in a magma chamber is effectively pattern formation that results from self-organization. When crystals can float or sink a significant distance, depending on melt viscosity and cooling rate, their interaction may result in the formation of barriers that define developing layers. Gravitational sorting determines the initial mineral distribution within each layer, which is repeated in subsequent layers. The complex interplay between crystal and heat transport can thus lead to rhythmic layering, but also to larger-scale layers or zonation within magma chambers.

This research was sponsored in part by the Ministerio de Economía y Competitividad of Spain under project ENE2012-30832 and the German Research Council (DFG). Soesoo acknowledges support by Estonian Science Foundation grant 8963. Walte acknowledges support by the German Federal Ministry of Education and Research project 05K13WC1. We highly appreciate support by and discussions with Thomas Wenzel and Gregor Markl and the critical comments from three reviewers.