Catastrophic growth of totally molten magma chambers in months to years.

The vertical growth rate of basaltic magma chambers remains largely unknown with available estimates being highly uncertain. Here, we propose a novel approach to address this issue using the classical Skaergaard intrusion that started crystallizing from all margins inward only after it had been completely filled with magma. Our numerical simulations indicate that to keep the growing Skaergaard magma chamber completely molten, the vertical growth rate must have been on the order of several hundreds to a few thousands of meters per year, corresponding to volumetric flow rates of tens to hundreds of cubic kilometers per year. These rates are several orders of magnitude higher than current estimates and were likely achieved by rapid subsidence of the floor rocks along faults. We propose that the Skaergaard is a plutonic equivalent of supereruptions or intrusions that grow via catastrophically rapid magma emplacement into the crust, producing totally molten magma chambers in a matter of a few months to dozens of years.

INTRODUCTION

Knowledge of magma emplacement time scales is critical for understanding how volcanic and igneous plumbing systems operate in the Earth’s crust. The time scales of processes that create space for magma chambers are particularly poorly constrained (, ). Estimates of the volumetric emplacement rates of plutons range from 10−4 to 10−1 km3/year (Fig. 1) (, ). Such estimates are mostly based on radiometric ages of minerals in the rocks of fossilized intrusions (–). Other estimates are based on measurements of surface deformation above growing intrusions (–). Both approaches are, however, indirect and intrinsically ambiguous. For example, since surface deformation measures emplacement rates in magmatically active areas, it provides information on short-term processes rather than an integrated view of long-term fluctuations of emplacement rates and repose periods. In layered mafic intrusions, the use of geochronological data is problematic because the high precision of U-Pb thermal ionization mass spectrometry zircon ages often contradicts basic field observations (–). In felsic plutons, large spreads of zircon ages from individual samples have led to the conclusion that the ages of many zircon grains do not date the timing of magma emplacement (). Geochronology must therefore be used with great caution for rigorous estimation of the growth rate of plutonic bodies. There is also a tacit assumption that magma chamber growth rates must be similar to rates of tectonic processes, which are inherently controlled by horizontal plate motion velocities (1 to 5 cm/year) (). However, at such slow rates, small batches of melts arriving into a chamber will totally or partially solidify between injections (), leaving little chance for the formation of large crystal-poor magma chambers (). Another approach is to use thermal models (Fig. 1) to determine magma emplacement rates that best reproduce observations such as the width of contact metamorphic aureoles or the absence of internal contacts (–). Here, we explore a simple and straightforward thermal approach to obtain the first rigorous estimates of the growth rate of the Skaergaard layered intrusion, Eastern Greenland—a classic example of a solidified basaltic magma chamber, whose study has, historically, constrained many of the fundamental principles of igneous petrology (–).

Fig. 1.

Volumetric flow rates for pluton emplacement.

Volumetric flow rates as compiled by Menand et al. () are estimated from radiogenic isotope geochronology, surface deformation, and thermal models. Also shown is our estimated volumetric flow rate for emplacement of the Skaergaard intrusion.

Skaergaard layered intrusion

The ca. 55–Ma (million year)–old Skaergaard layered intrusion is a fault-controlled igneous body made up of layered gabbro (–) (Fig. 2A) emplaced at a depth of about 2 km within the East Greenland volcanic rifted margin. The intrusion has a surface exposure of ~8 km by 11 km, a total structural height of ~4 km, and a box-like shape () with a total volume of ~300 km3. The parental tholeiitic melt of this intrusion was evolved [Mg/(Mg + Fe) < 0.5] (, ) relative to mantle-derived melts and with weak contamination by crustal rocks (high 87Sr/86Sri ~ 0.704) (, ), indicating likely derivation from deeper magma staging chambers (, , ). Closed-system crystallization of this parental melt took place from all margins inward with the formation of the Layered Series (LS) in the bowl-shaped floor, the Marginal Border Series (MBS) on the walls, and the Upper Border Series (UBS) below the roof. The UBS and LS meet at the Sandwich Horizon (SH). All three series record the progressive fractional crystallization of a parental tholeiitic melt (, ) as indicated by the successive appearance/disappearance of cumulus phases (e.g., augite, olivine, muscovite, apatite, and bustamite; Fig. 2B) and systematic changes in their compositions (e.g., anorthite content of plagioclase; Fig. 2B). Cumulate mineralogy is correlated between the series and is used to subdivide them into zones and subzones (UBS subzones are denoted with apostrophes, e.g., Lower Zone in UBS is LZa’, whereas MBS subzones are denoted with asterisks, e.g., Lower Zone in MBS is LZa*). The Hidden Zone (HZ) represents an unexposed portion of the LZa in the LS () that shows An content of plagioclase () similar to those from the uppermost portion of the LZa’ in the UBS (Fig. 2B) (). A crucial parameter for our modeling below is the amount of crystals delivered into and/or formed in the Skaergaard chamber before the onset of its bulk crystallization from all margins inward. Two points are most revealing in this respect. First, the melts that filled the Skaergaard chamber are equivalent to the evolved lavas of the contemporary Giekie Plateau Formation (, , ). They have the same phases on the liquidus (olivine and plagioclase), Mg/(Mg + Fe), TiO2/FeO* ratios, and geochronological age. These compositional constraints exclude other suggestions for parental melts [e.g., Milne Land Formation lavas (, )]. The evolved Giekie Plateau Formation lavas (, , ) are all notably aphanitic, suggesting that the melt entering the Skaergaard chamber had either a very low amount of or was totally free of phenocrysts. Some textural evidence interpreted earlier to indicate replenishment by plagioclase-rich slurries at the earliest stage of the Skaergaard chamber () can alternatively be attributed to the collapse of poorly consolidated material from the subvertical sidewalls/roof after the chamber had already started to crystallize (). Second, the LS, UBS, and MBS have similar chemical trends everywhere in the chamber, indicating that the Skaergaard intrusion started crystallizing inward from all margins only when it had been completely filled with melt. The lack of cumulates with the most primitive mineral composition at the base of the chamber, as indicated by the identical An content of plagioclase at the roof contact and in the lowermost portion of the HZ (~68 to 70%; Fig. 2B), suggests that no crystal growth and settling had occurred in the chamber before the onset of bulk crystallization at the margins. This implies that the body grew continuously to its current size by accumulation of compositionally near-homogeneous melt on a time scale that was much shorter than that of its solidification (). During this process, new melt effectively mixed with resident melt in the chamber, impeding the onset of crystallization. A key message is that the bulk of the intrusion eventually started to crystallize from one large volume of nearly crystal-free homogeneous melt. This fundamental physical constraint provides a unique opportunity to estimate the minimum rate of magma emplacement that was required to keep the Skaergaard magma body in a largely molten state (<<1% crystals) while growing to its current size.

Fig. 2.

Geology and stratigraphic section of the Skaergaard layered intrusion, East Greenland.

The Skaergaard intrusion is subdivided into three major petrological units: the LS, UBS, and MBS, crystallized concurrently along essentially parallel trends of differentiation. The LS formed on the floor, while the UBS crystallized under the roof and MBS at the wall. (A) East-West section through the Skaergaard intrusion. Modified from Nielsen et al. (). WSW, west southwest; ENE, east northeast. (B) The two fronts of crystallization converged at the SH as indicated by the mirror images of the cumulus phase arrivals and An content of plagioclase in LS and UBS. Modified from Hagen-Peter et al. () with some compositional points for the HZ from Holness et al. (). LZ, Lower Zone; MZ, Main Zone; UZ, Upper Zone; SL, Sea Level; Ol, olivine; Aug, augite; Mt, magnetite; Ap, apatite; Bu, (ferro)bustamite. An content = 100An/(An + Ab).

RESULTS

Constraining magma emplacement rates for the Skaergaard intrusion

We carried out a series of numerical simulations to determine the conditions required for the growth of a crystal-free Skaergaard magma chamber. For magma to remain at high temperatures, the conductive heat lost from the system into the colder, surrounding crust must be balanced by the heat added to the system by emplacement of new magma. Different shapes and modes of emplacement have been proposed for the intrusion (), including a laccolith that lifted the roof rocks (, ). Subsequent field-based studies have shown, however, that the box-like geometry of the intrusion, together with its vertical, fault-controlled walls, is more compatible with magma chamber growth by floor subsidence (, –). The Skaergaard intrusion has a high aspect ratio with its horizontal dimension being at least twice its vertical dimension (, ). Most of the magma heat content will have been lost through the roof and the floor, rather than the side walls, such that its cooling rate and crystal content will be controlled by the vertical emplacement rate or velocity (meters per year) (). Our numerical simulations are constrained by the characteristics of the Skaergaard intrusion, and they compute temperatures and crystal contents of the magma body for different vertical emplacement rates (Fig. 3). The modeled heat transfer through the country rock is conductive, and possible cooling caused by hydrothermal fluid circulation in the intrusion’s roof () is not taken into account. Our results are thus conservative; depending on the efficiency of cooling by hydrothermal circulation, Skaergaard’s emplacement rates might be even higher than those presented below.

Fig. 3.

The crystal fraction (weight %) versus magma emplacement rate (meters per year) in the Skaergaard magma chamber, East Greenland.

(A) The results of numerical simulation show that the minimum vertical emplacement rate that allows for a large, crystal-free chamber to grow without onset of crystallization along its margins is of the order of several hundreds to a few thousands of meters per year. The black stars indicate the positions of movies S1 to S3. (B) Setup of the numerical simulation. The inner rectangle (light yellow) shows the final extent of the intrusion. The numerical domain is indicated by the dashed blue line. qB, qR, and qL are the heat flux at the bottom, right, and left side of the domain. (C) Snapshot showing melt fractions and temperatures at the end of the intrusion’s growth. In this example, the emplacement rate is 500 m/year. The crystal content inside the intrusion (light yellow) is 0.4 wt %.

The relationship between crystal fraction and temperature is based on plagioclase saturation as reported by Thy et al. (). Choosing this relationship also produces conservative results as it is at the lower end of liquidus temperature values obtained from melting experiments (fig. S3). Higher solidus and liquidus temperature values result in higher crystal fractions in the magma chamber (fig. S4). The conductive heat flow out of the magma chamber depends on the temperature difference between the magma and the surrounding crust. Cold (30°C/km) and hot (100°C/km) initial geothermal gradients were tested, and we also ran a simulation to evaluate the possible thermal impact of lava flow accumulation before emplacement of the Skaergaard intrusion (fig. S1). Simulations were also run for different magma injection temperatures and possible magma superheating (fig. S2). The temperatures were calculated by solving the heat balance equation in cylindrical coordinates using a two-dimensional (2D) finite difference computational approach. Both conductive heat transfer and latent heat released by crystallization were included (table S1). The relationship between crystal fraction and temperature is constrained from petrological experiments (see fig. S3) (, ). To account for the effect of internal magma convection and mixing at low crystal fractions, temperatures were homogenized over the extent of the magma body where crystal fractions are <10 weight % (wt %). In the model, the intrusion grows vertically until it reaches a thickness of 4 km.

We find that for magma emplaced at its liquidus temperature (1160°C) (, ), the minimum vertical emplacement rate that permits a large, crystal-free chamber to grow without freezing at its roof and floor is extremely high, on the order of several hundreds to a few thousands of meters per year (Fig. 3C and movie S1). This corresponds to volumetric flow rates of tens to hundreds of cubic kilometers per year. Slower growth rates result in a crystal-rich magma chamber. For even slower rates (less than a few tens of centimeters per year), large parts of the intrusion are already solidified when the intrusion reaches its final size. The thermal modeling thus indicates that the volumetric flow rate of magma required for the Skaergaard magma chamber to remain crystal-free is three to five orders of magnitude larger than those typically inferred from the geochronology of felsic and mafic intrusions (, , ). This corresponds to an emplacement time for the entire mafic intrusion of a few years to a few decades. The results are similar for both cold and hot geotherms. The accumulation of lava flows before intrusion results in a temperature of 70°C at the estimated 2-km roof’s depth of the Skaergaard intrusion (fig. S1). This represents a small deviation from the initial 30°C geotherm and therefore does not affect the results. If the magma is superheated and emplaced at a temperature above its liquidus, the minimum emplacement rate reduces to a few tens of meters per year (for 15°C of superheating), resulting in emplacement durations on the order of centuries (fig. S1).

The validity of the model and the robustness of the results were tested by calculating the maximum temperatures reached in the country rock and by comparing them with temperature estimates from contact metamorphic mineral assemblages (Fig. 4). At high emplacement rates (≥5 m/year), the modeled temperature overlaps with data obtained from titanium-in-quatz thermometry, and the model closely reproduces the temperatures and distances to contact for the onset of country rock melting and the microcline out isotherm. At lower emplacement rates (<5 m/year), the modeled host-rock temperatures are too low. At even lower emplacement rates (<0.5 m/year), no large molten magma chamber forms (Fig. 3A), and using thermal aureoles to compare the model with observations loses its relevance. At low emplacement rates, aureole temperatures vary with depth and strongly depend on where the successive batches of magma are emplaced ().

Fig. 4.

Thermal aureole in models of the Skaergaard intrusion.

(A) Maximum temperatures reached for an emplacement rate of 500 m/year. The horizontal red line shows the location of the temperature profiles in (B). (B) Maximum temperatures in the country rock. The curves show the simulation results for different emplacement rates. The symbols show temperature estimates obtained by Bufe et al. () from titanium-in-quartz thermometry, the onset of country rock melting, and the microcline out isotherm. The temperatures estimated by Bufe et al. () are on samples collected on two profiles perpendicular to the vertical intrusion wall and located at an estimated depth of 2.3 ± 0.7 km below the intrusion’s roof.

Limits on magma emplacement rates from feeders and crustal deformation

Field mapping and forward modeling of gravity anomaly data indicate that the space for the vertical growth of the intrusion was accommodated by a down drop of its floor relative to the bounding normal faults (, , ). Lack of solid-state deformation fabrics (except for later fractures) with the Border Series gabbros adjacent to the faults () suggests that floor subsidence occurred while the intrusion was in a magmatic state and crystallized along undeformed margins. Gravity modeling has suggested the presence of a subvertical feeder zone below the central part of the intrusion (), but these may reflect anomalies caused by later dikes (). Possible feeder dikes include the contemporaneous swarm of Skaergaard-like dikes with lengths and widths on the order of 1 to 20 km and 10 to 100 m, respectively (, ). This provides further constraints to evaluate the feasibility of the high growth rates predicted by the thermal modeling.

Growth of the Skaergaard chamber can be approximately modeled by either subsidence of the floor or lifting the roof of a rectangular tabular body with horizontal dimensions ~8 km by 11 km. On the basis of field evidence (, , ), floor subsidence is the most likely mechanism for growth of the intrusion, but it should be noted that the calculations below are equally valid for both roof lifting and hybrid roof lifting/floor depression emplacement scenarios (, ). Piston- or cauldron-like subsidence is accommodated by displacement on subvertical faults, which could also act as conduits for magma flow from the source to the intrusion (–). We assume that the base of the underlying crustal block subsided into a depleting lower magma reservoir or stack of underlying reservoirs, possibly aided by ongoing regional extension (, ). Such reservoirs below the Skaergaard intrusion have not yet been detected but are expected to occur because deep staging chambers are inferred to underlie almost all shallow-level magma reservoirs (, ). In addition, the emplacement rates of the Skaergaard intrusion estimated above are simply too high to be attributed to the ascent of magma directly from the mantle. Such rates and the evolved nature of Skaergaard parental melt (, ) are more consistent with a rapid upward transfer of magma from mid- to upper-crustal staging chambers (). Taking minimum vertical growth rates dh/dt = 200 to 2000 m/year required to maintain a crystal-poor chamber (Fig. 3), the time to reach the ~4-km thickness of the Skaergaard intrusion ranges between 20 and 2 years, respectively. The corresponding volumetric magma inflow rates, Q, required to sustain these growth times vary between 17.6 and 176 km3/year (Fig. 5 and table S2).

Fig. 5.

Vertical growth rates (dh/dt) versus volumetric flow rates (Q) for growing a tabular, box-shaped magma chamber with dimensions 11 km (length) by 8 km (width) by 4 km (thickness) equivalent to the Skaergaard intrusion (thick grey line).

Dark red box (indicated as “minimum”) shows conditions for the lowest range of emplacement rates required to grow a crystal-free Skaergaard magma chamber constrained by thermal modeling (upper and lower dh/dt values indicated by dashed horizontal lines). The higher emplacement rates (indicated as “maximum”) are, however, possible but cannot be constrained by our thermal modeling. Vertical arrows indicate volumetric flow rates in a 1-km-long vertical feeder dike with widths of 1, 2, 5, and 10 m. Time-averaged volcanic eruption and felsic intrusion rates () are added for reference as are estimates of magma feeder and flow rates in flood basalt provinces and LIPs (, , ). Horizontal green box indicates the range of fault slip rates (indicated by dashed horizontal lines) estimated to accommodate caldera supereruption rates (indicated by dashed vertical lines) ().

We can now evaluate the viability of these time scales and rates by estimating the corresponding flow rates in the magma feeder zone(s) below the Skaergaard intrusion, which we approximate here as a single vertical dike. Volumetric flow rates in a dike-like feeder are driven by the density difference between the magma and the host rocks and are governed by the magma viscosity and the length and width of the conduit (). For a 1-km-long dike (i.e., the minimum length of dikes in the Skaergaard swarm) (), a conservative density difference of 50 kg/m3, and a typical, crystal-free basaltic magma viscosity of 1000 Pa s volumetric flow rates of 1.3, 10.3, 161, and 1290 km3/year are predicted for dike widths of 1, 2, 5, and 10 m, respectively (Fig. 5 and table S3). These widths are below the estimated feeder thicknesses for the Skaergaard intrusion, which, for the same parameters, could have sustained much higher volumetric flow rates. Reducing the density difference to 10 kg/m3 gives Q values of 0.3, 2, 32, and 260 km3/year for 1-, 2-, 5-, and 10-m-wide feeder dikes, respectively (table S3).

DISCUSSION

Feeder dikes in LIPs

Our calculations show that even a single, relatively short and thin feeder dike would have been capable of delivering magma into the growing Skaergaard intrusion at flow rates that fall well within those required for the melt to remain crystal free (tens to hundreds of cubic kilometers per year). The eruptive flow rates estimated for the basaltic lava flows associated with some Large Igneous Provinces (LIPs) (–) are well above our conservative estimates of Q, probably because the feeder dikes are longer and thicker, while other LIPs have similar values. Taking the Columbia River Basalts as an example, such flow rates range from 1 to 14 km3/day (365 to 5110 km3/year) for the Maxwell feeder dikes (), 20 to 40 km3/day (7300 to 14,600 km3/year) for the Wapshilla Ridge Member basalts (), 72 km3/day (26,280 km3/year) for Umatilla Member basalts (), and up to 100 km3/day (36,500 km3/year) for the Rose Member basalts (). It is also noteworthy that although the time-averaged rate of magma emplacement of basaltic sills and dikes in the Karoo LIP is estimated to have been ~0.78 km3/year, intrusion of individual shallow sills in this system may have occurred at rates of 60 to 300 km3/year (). Such pulsating events within a narrow time frame (, ) are in line with our calculations for the Skaergaard intrusion.

Fault slip rates and magma accommodation of rapidly growing magma chambers

The minimum vertical growth rates (Figs. 3 and 5) also require the minimum slip rates on the intrusion bounding faults to be on the order of hundreds to thousands of meters per year. Such slip rates are well above typical time-averaged velocities on active faults (approximately millimeters to centimeters per year) but well below those of single-slip earthquake-generating events, which have coseismic rupture velocities >0.1 m/s (). Both the fault slip and volumetric flow rates required for the rapid growth of the Skaergaard intrusion are comparable to estimates for catastrophic, ash-flow erupting caldera collapse events (). For example, calderas with similar surface areas to the Skaergaard intrusion (=88 km2; e.g., Atitlan, Guatemala = 113 km2; Crater Lake, USA = 50 km2) record >0.1 to >3 km of continuous vertical displacement in hours to days during certain stages of caldera collapse (–). These observations translate to slip rates of 0.1 to 10 cm/s (≡3.1 × 104 to 3.1 × 106 m/year) on the ring faults that permit crater floor down drop and corresponding eruption rates of 105 to 106 m3/s (≡3154 to 31,536 km3/year) (Fig. 5) (). Slip rates of these magnitudes are facilitated or lubricated by frictional melts that are generated on so-called “superfaults” formed during landslide, hypervelocity impact, and impact crater collapse events, in addition to ash-flow erupting calderas (a.k.a. “superuptions”) (). Hence, although the fault slip and volumetric flow rates required to maintain a crystal-free Skaergaard magma chamber seem implausibly high compared to time-averaged tectonic fault displacement and intrusion and volcanic eruption rates, they are one to two orders of magnitude less than those associated with catastrophic subsidence rates that occur during some ash-flow caldera eruptions. In the case of Skaergaard, very high vertical slip rates may have been facilitated by melt lubrication of the bounding faults ().

Growth rate of LIP magma chambers

The maximum eruptive rate reported for volcanic islands and mid-oceanic ridges () is ~0.1 km3/year. Assuming a conservative 10:1 intrusive:extrusive ratio, this corresponds to a maximum intrusive rate of ~1 km3/year for the underlying magma chamber. The volumetric flow rates estimated for the emplacement of Skaergaard are one to two orders of magnitude higher than these values, which suggests that the growth of Skaergaard-type magma chambers was an exceptionally fast geological event. Mafic layered intrusions, such as Skaergaard, are closely associated with LIPs that are characterized by the generation and emplacement of enormous volumes of mainly basaltic magmas (>1 million km3) over short periods of time (often <2 Ma) (, ). Layered mafic intrusions are an integral part of LIP plumbing systems (, ) through which magmas are either delivered to the Earth’s surface [open-system chambers, e.g., Bushveld Complex ()] or in which magmas stall and differentiate [closed-system chambers, e.g., Skaergaard intrusion ()]. Details of the processes leading to the rapid growth of the Skaergaard intrusion remain to be established. However, we can speculate that they are linked to the same massive magma input to the crust during LIP events that leads to large and rapid eruptions of basaltic lava at the surface. We can also surmise that vertical stacking of magma reservoirs that ultimately feed surface eruptions and the presence of faults enable the rapid and continuous ascent of magma. Once an open conduit (e.g., a dike) connects two vertically separated reservoirs, gravity-driven subsidence of the floor rocks of the overlying chamber into a depleting lower chamber may be facilitated by subvertical and possibly melt-lubricated faults (, , ). Such an exchange of magma between the reservoirs may facilitate the vertical growth of the overlying chamber (and simultaneous vertical compaction of the lower chamber) at rates consistent with those estimated here for the Skaergaard intrusion.

Similar to Skaergaard, many layered mafic intrusions tend to have cumulate stratigraphy that displays systematic crystallization sequences and mineral chemical trends indicating their crystallization in large and mostly molten magma chambers in Earth’s crust (, , , –). This suggests that the growth rate of the magma chambers of these intrusions is likely comparable to that of Skaergaard and that the filling times of such bodies might be unexpectedly fast. For instance, taking a volumetric flow rate of 100 km3/year, a body of the size of Skaergaard (300 km3) can be filled in a matter of only 3 years. Since our numerical simulations constrain the minimum rate of magma emplacement, the actual filling time for such a body may be even faster, perhaps several months or even a few weeks. Following this logic, for bodies the size of Kiglapait (3500 km3) () and Windimurra (26,000 km3) (), the filling times could have been <35 and 260 years, respectively. Again, the filling times for these intrusions may be on the order of a few months if greater than minimum chamber growth rates are used. Therefore, leaving aside the matchless giant Bushveld Complex (600,000 km3) (), the actual filling times for many layered intrusions may be in the range of only a few months to a maximum of dozens of years. This conclusion is valid for the cases of magma chamber growth from both crystal-free [e.g., Skaergaard ()] and crystal-bearing magmas [e.g., Stillwater ()]. These estimates are consistent with the durations of volcanic eruptions, which most commonly last from about a month to a year, with the overall spectrum of times spanning about 100 years (, ). We propose that the Skaergaard and possibly some other layered intrusions can be viewed as products of catastrophic events associated with a rapid and voluminous emplacement of basaltic magmas into the Earth’s crust. In other words, some layered mafic intrusions may represent the plutonic analogs of LIP-related volcanoes that are responsible for eruptions of enormous volumes of flood basalts on the Earth’s surface (, , ).

MATERIALS AND METHODS

Thermal modeling

Temperatures and melt fractions are computed on a 2D numerical grid. The dimension of the cells is 50 m by 50 m. The size of the numerical domain is a function of the simulation duration and hence of the emplacement rate. It is calculated using the heat diffusion length to avoid boundary effects. The coordinates are cylindrical ().

The emplacement of magma is simulated by setting the corresponding cells to liquidus temperature or higher temperatures in case of overheating. Convection is simulated by averaging the heat content over the cells with less than 10 wt % crystal. The heat transfer between cells is calculated with the finite difference expression of the heat balance equationwith

H is the heat content per unit volume, t is the time, k is the thermal conductivity, T is the temperature, ρ is the density, c is the specific heat capacity, X is the melt fraction, and L is the latent heat of crystallization. The values of the parameters used in the simulations are reported in table S1. The conductivity depends on temperature. We used the equation from () for the basaltic roof and for the Skaergaard magma, and the equation from () for the gneiss below the intrusion.

The relationship between temperature T and crystal fraction X is based on liquidus temperatures estimated by Thy et al. () from Skaergaard plagioclase compositions (fig. S3). T-X data points used for the simulation were read from Fig. 4 in (), and the T-X relationship is interpolated linearly between the data points (fig. S3). The T-X relationship is sensitive to the assumed composition of the parental magma and to the oxygen fugacity. Our chosen relationship is conservative in terms of minimal emplacement rates. Figure S4 shows that a T-X relationship that involves higher solidus and liquidus temperatures, based on experiment FG1_AC of Thy et al. (), results in higher crystal fractions in the magma chamber. We used a constant latent heat of crystallization per mass of crystal crystallized without distinguishing between the different phases. The latent heat released is not linear, since it depends on the mass crystallized, which itself depends on the crystallinity-temperature relationship. Considering the uncertainty on the crystallinity-temperature curve and on the crystallizing phases, neglecting the difference in latent heat release of the different phase does not affect the accuracy of the results.

The intrusion grows as a series of twenty 200-m-thick sills. Each new sill is emplaced at the top of the magma chamber. In the absence of magma chamber, sills are added on top of former ones. Different emplacement rates were tested by varying the time between sills. The space for new magma is accommodated by moving the cells below the new sill downward. The numerical domain size remains constant so that cells disappear at the bottom of the domain with new injections.

The boundary condition at the top of the domain has a fixed temperature of 0°C. The initial temperatures increase downward following a geothermal gradient of 30°C/km (cold geotherm) or 100°C/km (hot geotherm). The boundary condition at the bottom of the domain is a heat flux qB that is determined by assuming that temperatures are restored to the initial geothermal gradient 2 km away from the boundary. The left and right boundaries are insulating (no heat flux).

The roof of the Skaergaard intrusion is formed by lava flows. We tested the possible effect on the initial temperatures of lava flow accumulations before the intrusion. This was done in 1D by adding layers at magmatic temperatures at the top of the numerical domain. We found that after the emplacement of 2 km of lava over 200,000 years, the temperature is 70°C at the base of the lava column, which represents only a slight departure from the initial cold geothermal gradient of 30°C/km (fig. S1).

Intrusion growth and magma supply rates

The vertical growth rate dh/dt and volumetric infilling rate Q of an h-m-thick horizontal tabular intrusion bounded on all sides by vertical faults iswhere A is the surface area of the intrusion given by length × width for a rectangular body or πr2 for a circular body of radius r (). The time to grow an intrusion to vertical thickness h is

The volumetric flow rate in a vertical dike feeding a tabular intrusion is given bywhere Δρ is the density difference between the host rocks and the magma, μ is the magma kinematic viscosity, g is the gravitational acceleration, and a and b are the dike width and length, respectively (). To calculate Q in a vertical feeder dike, we assume a typical crustal density of 2800 kg/m3 and a range of basaltic magma densities between 2750 and 2790 kg/m3 (), which gives Δρ = 50 to 10 kg/m3. We also assume a typical basaltic magma kinematic viscosity of 1000 Pa s and a minimum feeder dike length b = 1 km. Q is then calculated for dike widths a = 1, 2, 5, and 10 m (table S3).