Deep open storage and shallow closed transport system for a continental flood basalt sequence revealed with Magma Chamber Simulator

The Magma Chamber Simulator (MCS) quantitatively models the phase equilibria, mineral chemistry, major and trace elements, and radiogenic isotopes in a multicomponent–multiphase magma + wallrock + recharge system by minimization or maximization of the appropriate thermodynamic potential for the given process. In this study, we utilize MCS to decipher the differentiation history of a continental flood basalt sequence from the Antarctic portion of the ~ 180 Ma Karoo large igneous province. Typical of many flood basalts, this suite exhibits geochemical evidence (e.g., negative initial εNd) of interaction with crustal materials. We show that isobaric assimilation-fractional crystallization models fail to produce the observed lava compositions. Instead, we propose two main stages of differentiation: (1) the primitive magmas assimilated Archean crust at depths of ~ 10‒30 km (pressures of 300–700 MPa), while crystallizing olivine and orthopyroxene; (2) subsequent fractional crystallization of olivine, clinopyroxene, and plagioclase took place at lower pressures in upper crustal feeder systems without significant additional assimilation. Such a scenario is corroborated with additional thermophysical considerations of magma transport via a crack network. The proposed two-stage model may be widely applicable to flood basalt plumbing systems: assimilation is more probable in magmas pooled in hotter crust at depth where the formation of wallrock partial melts is more likely compared to rapid passage of magma through shallower fractures next to colder wallrock.


Introduction
Differentiation of magmas in the crust has been one of the most studied topics in petrology since Daly (1914) and Bowen (1928) argued the relative importance of assimilation and fractional crystallization (FC) in producing the compositional diversity of igneous rocks. Models of closed system crystallization-equilibrium and fractional-entered Communicated by Mark Ghiorso.

Electronic supplementary material
The online version of this article (https ://doi.org/10.1007/s0041 0-019-1624-0) contains supplementary material, which is available to authorized users. the realm of chemical thermodynamics and phase equilibria in the late twentieth century following the earlier pioneering studies of N.L. Bowen and others and were accompanied by the introduction of various computational tools (e.g., Nathan and Van Kirk 1978;Ghiorso et al. 1983;Ghiorso 1985;Ariskin et al. 1993;Ghiorso and Sack 1995;Gualda et al. 2012). In contrast, concurrent assimilation-fractional crystallization (AFC) models lagged behind and have largely relied on assessment of trace elements and isotopic ratios without consideration of detailed phase equilibria, major elements, and conservation of energy (Taylor et al. 1979;Taylor 1980;DePaolo 1981). The energy-constrained AFC (EC-AFC) models Spera 2001, 2007;Spera and Bohrson 2001), including the Magma Chamber Simulator (MCS; Bohrson et al. 2014), are able to model crystallization, magma recharge (R), and assimilation processes more comprehensively. In particular, by coupling thermodynamically based phase equilibria solutions to constraints based on trace elements and isotope geochemistry, a more detailed and complete picture of magma petrogenesis can be obtained, in part, because dozens of models (hypotheses) can be efficiently generated and carefully compared against existing petrological, geochemical, and field data.
Continental flood basalts (CFBs) are remnants of the largest subaerial lava eruptions on Earth. Their petrogenesis is complex, in part, because they usually exhibit geochemical evidence of considerable interaction with (or derivation from) continental lithospheric materials (e.g., Lightfoot et al. 1990;Hooper and Hawkesworth 1993;Jourdan et al. 2007a).
Many CFBs also show evidence of polybaric crystallization ± crustal assimilation (e.g., Gibson 1990;Fram and Lesher 1997;Hole 2018), which are difficult to decipher in detail without thermodynamic treatment of phase equilibria. Whereas traditional AFC trace element models have suggested unreasonably high amounts of crustal assimilation for basaltic magmas (up to 60%; see Molzahn et al. 1996;Ewart et al. 1998;Larsen and Pedersen 2009), some recent EC-AFC models have shown that basalts can inherit strong crustal signatures by much lower and more thermodynamically feasible amounts of crustal contribution (< 20%; e.g., Fowler et al. 2004;Jourdan et al. 2007a;Heinonen et al. 2016). This is mainly because thermodynamic models such as EC-AFC and MCS take into account the dynamics of wallrock partial melting, which may effectively enrich incompatible elements, thereby profoundly influencing also the radiogenic isotope compositions of contaminated magma. In contrast, older models (without explicit energy conservation) generally mix constant composition contaminants into magma, without regard for the compositional variability that attends partial melting of wallrock. Thermodynamic models clearly have the advantage over models that achieve precise mass balance but not at all or only roughly incorporate energy considerations.
In this study, we focus on a well-studied flood basalt type (low-ε Nd CT1 of Luttinen and Furnes 2000) from the Antarctic portion of the Karoo large igneous province (LIP) (Fig. 1). Preliminary EC-AFC models that approximated energy balance indicated that assimilation of ≤ 15%  (Luttinen and Siivola 1997;Jourdan et al. 2007b) and the Archean/Paleoproterozoic core complex of Kalahari (Jacobs et al. 2008). Division into North and South Karoo CFBs after Luttinen (2018); for a more detailed map of Karoo LIP exposure and its divi-sion North and South subprovinces, see that reference. b Distribution of Jurassic CFBs and the Precambrian basement in western Dronning Maud Land. Plogen nunatak in Vestfjella indicated. Lithospheric boundary reconstructed after Corner (1994). H.U.S. H.U. Sverdrupfjella of Archean crust by a depleted mantle-derived parental magma is sufficient to explain the incompatible crust-like trace element and Nd isotopic characteristics of the basalts (Heinonen et al. 2016). The present study attempts to constrain the magmatic evolution in more detail by simultaneous modeling of major elements and phase equilibria using the fully thermodynamic self-consistent MCS formulation, which supersedes the older EC-AFC class of models. In addition, thermophysical considerations for the magma ascent network are provided. When assimilation is called upon as part of a petrogenetic scenario, the mechanism by which contamination occurred is rarely specified. Here, we show that, in addition to mantle source heterogeneity, magma transport via a two-stage magma fracture network provides a potential explanation for some of the phase equilibria, major element, trace element, and isotope characteristics of the low-ε Nd CT1 lavas. This analysis illustrates the effectiveness of coupling petrologic modeling with magma transport phenomena to simulate the petrogenesis of complex magma systems.

Regional geology
The generation of the ~ 180 Ma Karoo LIP marked the initial rifting stages of the Gondwana supercontinent (e.g., Jourdan et al. 2005). The Karoo CFBs and associated intrusive rocks consist of several geochemically distinct magma types (initial ε Nd from − 17 to + 9; e.g., Riley et al. 2005;Luttinen et al. 2015) (Fig. 2). Most of the geochemical variation is found among the CFBs and associated dikes close to the socalled fossil "triple rift" structure (recently coined as South Karoo subprovince; Luttinen 2018) that is marked by abundant dike swarms largely parallel to Precambrian zones of weakness in the lithosphere (Jourdan et al. 2006). The South Karoo subprovince also includes most of the primitive and Mg-rich magma types described from the Karoo LIP (e.g., Ellam and Cox 1989;Sweeney et al. 1991;Riley et al. 2005;Heinonen et al. 2010). On the contrary, the widespread and voluminous Karoo lava and sill formations that are found on or within the surrounding sedimentary basins (North Karoo subprovince; Luttinen 2018) are geochemically fairly homogeneous and more evolved (e.g., Neumann et al. 2011).
The Vestfjella CFBs ( Fig. 1) belong to the South Karoo subprovince and may be among the oldest Karoo LIP lavas, because the cross-cutting dikes record 40 Ar/ 39 Ar ages as old as ~ 189 Ma (Luttinen et al. 2015). The CFBs overlie a complex basement that consists of an Archean craton (Grunehogna craton that corresponds to the Kaapvaal craton in southern Africa; e.g., Marschall et al. 2010), which is largely overlain by Proterozoic igneous and metasupracrustal rocks (e.g., Moyes et al. 1995) in the north and a Proterozoic metamorphic belt (Maud Belt; e.g., Jacobs et al. 1998Jacobs et al. , 2003 in the south (Fig. 1). In addition, layers of Phanerozoic fossilbearing sedimentary rocks (e.g., Lindström 1995) are found sporadically overlying the Precambrian rocks. For a more detailed overview of the generalized geology of the area, the reader is referred to Wolmarans and Kent (1982) and Groenewald et al. (1995).
The geochemically diverse Vestfjella CFBs can be divided into four different chemical types (low-ε Nd CT1, high-ε Nd CT1, CT2, and CT3; Luttinen and Furnes 2000). The low ε Nd -CT1 and CT3 are the most common and well characterized. In addition to the CFBs, Vestfjella also contains geochemically variable Karoo LIP dikes (Luttinen et al. 2015) that include magma types that are isotopically indistinguishable from modern mid-ocean ridge basalts (MORBs) of the Southwest Indian Ridge (SWIR) and likely sampled depleted upper mantle sources (Figs. 2, 3; Vestfjella depleted ferropicrites and Low-Nb; Heinonen et al. 2010). Preliminary trace element and Nd isotopic modeling using EC-AFC equations of Bohrson and Spera (2007) (Hawkesworth et al. 1984;Ellam and Cox 1989;Harris et al. 1990;Sweeney et al. 1994;Luttinen et al. 1998Luttinen et al. , 2010Luttinen and Furnes 2000;Riley et al. 2005Riley et al. , 2006Jourdan et al. 2007a;Neumann et al. 2011). Compositions of depleted MORB mantle (DMM; Workman and Hart 2005), and SWIR MORB that are most isotopically similar to the low-Nb dikes (Hamelin and Allegre 1985;le Roex et al. 1989;Mahoney et al. 1989Mahoney et al. , 1992Janney et al. 2005; see Heinonen and Kurz 2015) are also indicated. DMM and SWIR MORB have been back-calculated at 180 Ma using Rb/Sr and Sm/Nd of the DMM reservoir (Workman and Hart 2005) degree partial melts from depleted sources similar to those that formed these dikes are promising candidates for the parental magmas of the four different types of Vestfjella CFB lavas ( Fig. 3; Heinonen et al. 2016).

Low-ε Nd CT1 lavas
The low-ε Nd CT1 lavas are mainly found in northern Vestfjella (Luttinen and Furnes 2000). The great majority of the lavas are subalkaline tholeiitic basalts or basaltic andesites: MgO contents range from 5 to 13 wt% in samples that do not show evidence of olivine or pyroxene accumulation. The lavas exhibit characteristically low initial ε Nd (from − 16 to − 11) (Fig. 2). According to the high-Ti vs. low-Ti classification scheme of Erlank et al. (1988) (high-Ti: Ti/Y > 410, Zr/Y > 6; low-Ti: Ti/Y < 410, Zr/Y < 6), the low-ε Nd CT1 lavas show low-Ti character in terms of Ti/Y (240-360), but transitional character in terms of Zr/Y (3.9-7.4). Several dikes of low-ε Nd CT1 character are also known from Vestfjella; they post-date the generation of the lavas and have been interpreted to record relatively lower degrees of mantle partial melting (Luttinen et al. 2015).
Most of the low-ε Nd CT1 lavas are porphyritic, but all samples considered in this study have been collected from the massive interior parts of the flow units, lack cumulate textures, and are thus considered to closely correspond to melt compositions (Luttinen and Furnes 2000). The predominant phenocrysts are in the order of abundance: plagioclase, clinopyroxene, and olivine. Groundmass, when crystals are optically visible, consists of microlites of clinopyroxene, plagioclase, and Fe-Ti oxides. Plagioclase is usually the only phenocryst phase, although in more primitive samples, it is accompanied by clinopyroxene or olivine, or both. Olivine phenocrysts (up to Fo 90 ;Luttinen 2000;Heinonen et al. 2018) have generally been altered to secondary minerals. Overall major element oxide variations have been explained by a fractionation sequence of olivine + augite (7:3; when Mg# > 0.55) and olivine + plagioclase + augite (2:5:3; when Mg# < 0.55) (Luttinen and Furnes 2000). Orthopyroxene phenocrysts (Mg# up to 84;Luttinen 2000), however, are also known from the most primitive (MgO > 10 wt%) lava flows suggesting some role for orthopyroxene in the differentiation as well. These lavas are highlighted in Figs. 4,5,6,7,8,9. In projection from olivine into the plane CS-MS-A in the QMAS system of O'Hara (1968), most of the lavas plot close to the 0.1 MPa (1 bar) olivine-augite-plagioclase cotectic and far from high-pressure cotectics, which suggests last equilibration at low pressure (Fig. 5). Whereas the other CT magma types seem to record assimilation of lithospheric mantle materials and possibly lower crust based on their incompatible trace element and Nd isotopic compositions, the low ε Nd -CT1 magmas most likely assimilated Archean crust ( Fig. 3; Luttinen et al. 1998;Luttinen and Furnes 2000;Heinonen et al. 2016).
The most extensive continuous stratigraphic cross-section of Vestfjella CFBs, about 900 m, is found at the Plogen nunatak (Figs. 1, 4). This particular cross-section includes all four magma types of the Vestfjella CFB lavas and illustrates the Primitive-mantle-normalized (Sun and McDonough 1989) incompatible trace element patterns shown for a the uncontaminated MORB-like Low-Nb dikes (Luttinen and Furnes 2000;Heinonen et al. 2010Heinonen et al. , 2016Luttinen et al. 2015; isotopically, the most primitive sample P27-AVL highlighted in dark blue) and b representative low-ε Nd CT1 lavas (Luttinen and Siivola 1997;Luttinen et al. 1998;Luttinen and Furnes 2000; picritic orthopyroxene-bearing sample P2-AVL highlighted in orange) and best-fit EC-AFC trace element models with an Archean TTG assimilant presented by Heinonen et al. (2016). Average N-MORB (Sun and McDonough 1989) shown for reference in a; residual garnet is responsible for the relatively lower Y, Yb, and Lu concentrations in the low-Nb dikes relative to average MORB. Parental melt constructed on the basis of the low-Nb dikes for the EC-AFC trace element models presented by Heinonen et al. (2016) is shown in both subfigures variability of the CFB magma transport network and long travel distances of the low-viscosity basaltic lavas from their currently unexposed eruptive centers (Luttinen and Furnes 2000). The cross-section is dominated by low-ε Nd CT1, especially in the lower and upper parts of the section (Fig. 4). We use this cross-section as a reference for temporal evolution of the low-ε Nd CT1 magma type during Vestfjella CFB magmatism. Whereas some trace element ratios (e.g., Ti/Y) show sublinear variation in relation to stratigraphic height, others (e.g., Zr/Y) and major elements (e.g., MgO) do not and are poorly correlated with it ( Fig. 4). In summary, low ε Nd -CT1 is an ideal suite for studying the differentiation history of an early CFB magma type, because there is an abundance of temporally constrained geochemical data available (Figs. 2,3,4) and clear indications for crustal assimilation having a role in its petrogenesis (

Magma Chamber Simulator
The Magma Chamber Simulator (Bohrson et al. 2014) combines EC-AFC equations introduced by Spera and Bohrson (2001) with multicomponent-multiphase models for silicate solid-liquid systems (rhyolite-MELTS software; Gualda et al. 2012;Ghiorso and Gualda 2015;v.1.2.0 used in this study). The input includes major element composition, mass, temperature, and pressure specification for the parental melt (magma), wallrock, and up to five recharge magmas. In addition, a percolation threshold value for mobilization of wallrock partial melts is specified. Solidus and liquidus temperatures and initial melt Fe 2 O 3 /FeO of the different subsystems at given pressures can be constrained with the help of MELTS software. Ti/Y

(a)
Low-CT1 ε Nd with opx High-CT1 ε Nd CT2 CT3 Variations in a MgO, b Ti/Y, and c Zr/Y in the Plogen section of the Vestfjella CFBs. Linear trends with R 2 values are given for the low-ε Nd CT1 lavas above and below 644 m in b and c. Arrows F and A in b and c approximate relative influences of increase in degree of melting in the source and assimilation with TTG crust, respectively; F is more important for variation in Ti/Y, whereas assimilation is more important for variation in Zr/Y. Data sources: Luttinen and Siivola (1997), Luttinen et al. (1998), and Luttinen and Furnes (2000); the full dataset of low-ε Nd CT1 lavas used in this study is available in Online Resource 3 After priming the behavior of wallrock at near-solidus temperatures, MCS proceeds as follows in an EC-AFC run: (1) Fractional crystallization of the parental melt (that composes the resident magma body) occurs until the wallrock has been heated above its solidus and the percolation threshold of wallrock partial melt has been exceeded. More specifically, FC is modeled as individual equilibrium crystallization steps within a temperature range defined by the user. In the end of each temperature step, all equilibrium solids are fractionated from the melt to form cumulates. (2) After the solidus temperature and melt percolation threshold of the wallrock have been exceeded, wallrock partial melts are introduced into the resident magma body. Wallrock changes its modal and chemical compositions as the melting proceeds-a fraction of melt defined by the percolation threshold always remains within the wallrock. The fraction of wallrock melt above the percolation threshold is mixed into and equilibrated with the resident magma such that the major element composition and equilibrium mineral phases are updated accordingly. (3) Subsequent fractionation and assimilation steps follow until the run ends based on userdefined mass or temperature constraints or when thermal  Trace element and isotope models can be implemented after evolution of major elements and mineral phases have been constrained. The trace element calculations utilize the solid-melt-fluid equations of Spera et al. (2007) in each crystallization/melting step. MCS output records the major and trace element and isotope evolution of melt and crystals (and associated fluid phase, if present) in resident magma, cumulates, and wallrock.

Constraining the thermodynamic and geochemical parameters for MCS modeling
All low-ε Nd CT1 lavas and Archean wallrock data used for MCS modeling have been published in previous studies (Luttinen and Siivola 1997;Luttinen et al. 1998;Kreissig et al. 2000;Luttinen and Furnes 2000;Heinonen et al. 2010). The full dataset of low-ε Nd CT1 lavas used in this study is available in Online Resource 3.

Parental melt composition
In a previous study that utilized the EC-AFC model of Bohrson and Spera (2007), we showed that the low-ε Nd CT1 trace element and Nd isotopic compositions can be readily modeled by incorporation of Archean crust into a depleted mantle-derived parental melt (Heinonen et al. 2016). The parental melt was constrained on the basis of the MORB-like Low-Nb dikes found in the area (Fig. 3). This study focuses on the compositional variability and magmatic plumbing system of the low-ε Nd CT1 magma type and, therefore, we use a parental melt composition that has been constrained on the basis of the most depleted and primitive low-ε Nd CT1 lavas that have already experienced some degree of crustal contamination.
The low-ε Nd CT1 parental melt (PM; Table 1) is an average of the most primitive non-cumulate lava samples (P1-AVL, P2-AVL, and P63-AVL; Luttinen and Furnes 2000), except that the lowest values are used for Nd and Y (that show anomalous enrichment in some samples) and the highest values recorded from the low-ε Nd CT1 suite are used for 143 Nd/ 144 Nd (180 Ma) (olivine cumulate sample AL/307;   Heinonen et al. 2010). The Fe oxidation state is not known, but we used the commonly assigned Fe 2+ /Fe tot = 0.85 (e.g., Luttinen and Furnes 2000) as the initial oxidation state for almost all runs (the runs themselves were not buffered for oxygen). To test the effect of reasonable variation in the PM oxidation state, some models with respective values 0.8 and 0.9 were also ran (Online Resource 1). Another unknown but potentially significant factor is the water content of the parental melt. Because the ultimate source of the Vestfjella lavas likely resided in moderately water-enriched depleted mantle (Heinonen et al. 2016(Heinonen et al. , 2018 Phase-specific and temperature-independent partition coefficients used for trace element modeling were estimated based on listings given for mafic melt compositions in EarthRef database (https ://earth ref.org/KDD/). When partition coefficient data were lacking, similar phases and/ or elements were used to estimate the values. Partition coefficients are listed in the MCS output files given in Online Resource 2.

Wallrock composition
Based on the very low ε Nd of the lavas and the previous EC-AFC trace element and isotopic modeling (Heinonen et al. 2016), the most plausible assimilant of the low-ε Nd CT1 lavas is Archean crust. Because of very limited exposure and lack of comprehensive compositional data on the Archean rocks of western Dronning Maud Land (Marschall et al. 2010), we selected the assimilant among the rocks related to the Kaapvaal craton in southern Africa. A tonalite-trondhjemite-granodiorite (TTG) sample 96/203 (Table 1) reported by Kreissig et al. (2000) provides the best fit in trace element EC-AFC models (Heinonen et al. 2016) and we also adopted its Nd isotopic composition (calculated at 180 Ma). As in the case of PM, Fe oxidation state and H 2 O contents for the wallrock are not known. Because the TTGs reported by Kreissig et al. (2000) contain both magnetite and ilmenite as accessory phases, we used Fe 2 O 3 /FeO of Fig. 8 Results of MCS AFC + FC (500 MPa, MgO > 9 wt%) + FC (0.1 MPa, MgO < 9 wt%) models for PM and TTG assimilant (Table 1) shown for major elements together with the low-ε Nd CT1 whole-rock data (for sources, see Fig. 4). The yellow filled circles mark the transitional compositions for the high-P AFC/FC and low-P FC models. MCS FC model at 0.1 MPa (see Fig. 6) also shown for reference. All data and models normalized to 100 wt% volatile-free. Lavas containing orthopyroxene phenocrysts and MgO range for model melts in equilibrium with orthopyroxene indicated. Note that the two lavas within Stage 1 MgO range but without orthopyroxene phenocrysts may have crystallized at lower pressures earlier than suggested by the two-stage model. The yellow field covers the compositions bracketed by the different models. See Sect. "Polybaric two-stage modeling" for more details. 0.5, which has been suggested as the divisive value between ilmenite-and magnetite-series granitoids (Ishihara 1981), both of which are known from the Kaapvaal craton (Ishihara et al. 2006). An average bound H 2 O content of 0.6 wt% from the latter study is also applied here. Preliminary testing indicated that varying Fe oxidation state and H 2 O content of the wallrock within reasonable limits have negligible effects on the results of the MCS AFC models (the only significant effect of adding H 2 O is to lower the solidus temperature of the wallrock). This is because all such models are dominated in weight by the resident magma. Partition coefficients were estimated based on listings in the EarthRef database (but this time data for felsic melt compositions were used) and are listed in the MCS output files given in Online Resource 2.

Pressure, temperature, and mass constraints
The pressures used in the models were 0.1 MPa (atmospheric pressure, only used in pure FC models), 100 MPa (depth of ~ 4 km), 300 MPa (~ 10 km), 500 MPa (~ 20 km), and 700 MPa (~ 30 km). The depths have been calculated to one significant digit to allow for uncertainty in crustal density (2800 kg/m 3 assumed in the calculation). We did not run models at higher pressures, because at and above 700 MPa, olivine is not a stable phase in the fractionating parental magma (primitive Fo 90 olivine phenocrysts have been documented from the lavas; Heinonen et al. 2018). Note that the same parental melt and wallrock compositions are used in every model. The initial temperature of PM was the liquidus temperature as defined by MCS in the beginning of the run. For the wallrock, we assumed that it had already been heated to its solidus temperature by earlier magma pulses and used an initial temperature few degrees below to its solidus at a given pressure. A default melt percolation value of 10% was used for the wallrock (see discussion of melt percolation thresholds in Bohrson et al. 2014). For the purposes of modeling, the mass of the parental melt and wallrock was always 100 units. In the case of AFC in a hypothetical spherical magma body that has a diameter of 1 km, this would mean a wallrock aureole with a thickness of ~ 100 meters.

Modeling the differentiation of the low-ε Nd CT1 flood basalt sequence using MCS
The results of the MCS FC and AFC models are discussed below, illustrated in Figs. 6, 7, 8 and 9, and all the input and output data are included in Online Resource 2 as tabulated Fig. 9 Results of MCS AFC/FC (500 MPa, MgO > 9 wt%) + FC (0.1 MPa, MgO < 9 wt%) models for PM and TTG assimilant (Table 1) shown for trace elements Zr and Y and ε Nd together with the low-ε Nd CT1 wholerock data (for sources, see Fig. 4). The yellow circles mark the transitional compositions for the high-P AFC/FC and low-P FC models. MCS FC model at 0.1 MPa (see Fig. 6) also shown for reference. All data and models normalized to 100 wt% volatile-free. Lavas containing orthopyroxene phenocrysts and MgO range for model melts in equilibrium with orthopyroxene indicated. The yellow field covers the compositions bracketed by the different models. See Sect. "Polybaric two-stage modeling" for more details. For full input and output, see Online Resource 2 worksheets. Fractional crystallization models with varying Fe oxidation state and H 2 O content of PM and an AFC model with five recharge events (R 5 AFC) are illustrated in Online Resource 1. In the following discussion, we concentrate on the relationships of MgO, SiO 2 , Al 2 O 3 , TiO 2 and FeO tot since these are the most abundant major element oxides and have not been notably affected by secondary alteration (CaO, Na 2 O, and K 2 O have been affected to some degree; see Luttinen and Furnes 2000). Trace element (Zr and Y) and Nd isotope models are only presented for the best-fit major element models; Nd is the only isotope system for which there is representative data available for both the lavas and the wallrock.
Any presented model should be consistent with two key observations from the major element data and modal composition of the low-ε Nd CT1 lavas: (1) stability of orthopyroxene in the most primitive lavas (MgO > 12 wt%) and (2) the inflection point at MgO ≈ 7 wt% in major oxide diagrams marking the onset of plagioclase fractionation (Fig. 6).

Fractional crystallization (FC) modeling
It is evident from the FC model results that pressure has a significant effect on the phase equilibria and major element evolution of the parental melt (Fig. 6). The results illustrate a dilemma: the "plagioclase inflection point" in the Al 2 O 3 vs. MgO trend is best reproduced by low-pressure models (0.1 and 100 MPa), but orthopyroxene is only a liquidus phase in picritic melts in models with pressures above 300 MPa.
The effects of varying Fe oxidation state and H 2 O contents in PM are illustrated in Online Resource 1. Summarizing, the effect of Fe oxidation state on the phase equilibria and major element evolution of the models is negligible. Lower MORB-like H 2 O contents stabilize early orthopyroxene at 300 MPa or higher pressure, but higher H 2 O contents akin to purported average values for Karoo CFBs stabilize early orthopyroxene only at 500 MPa or higher pressure. Higher H 2 O contents in PM also delay the onset of plagioclase fractionation.
These observations effectively preclude isobaric FC to explain the evolution of the low-ε Nd CT1 magmatic system.

Assimilation: fractional crystallization (AFC) modeling
On the basis of considerable incompatible trace element enrichment in the low-ε Nd CT1 lavas (Fig. 3), crustal contamination of their parental melts was dominated by interactions with wallrock partial melts instead of bulk assimilation (Heinonen et al. 2016). The MCS models involving assimilation of wallrock partial melts are illustrated in Fig. 7. Such simulations result in profound differences in the chemistry and phase equilibria of the magma relative to the FC-only model at an equal pressure.
In terms of major oxides, the AFC models result in higher SiO 2 , lower or similar Al 2 O 3 , and lower TiO 2 and FeO tot compared to the FC model at the same pressure at a given MgO (Fig. 7). This is because the main contributing phases to the wallrock partial melt are quartz and Si-rich albitic plagioclase. The stability of orthopyroxene in the resident magma increases because of Si added from the wallrock, and orthopyroxene is stable even in low-Mg magmas at high pressures. In contrast, the stability of plagioclase in the magma is significantly decreased, because Ca and Al activities in the resident magma are lower for the AFC models than for the FC models at the same MgO. The models do not thus replicate the inflection point seen in the data. Table 1 Parental melt (PM) and wallrock (WR) compositions used in most of the models Note that the oxides do not sum up to 100 wt%. This is because for the models with varying Fe oxidation state and H 2 O (see text), the listed values for the other oxides were always used before normalization. In MCS, the compositions are always normalized to 100 wt% before the run. The above presentation guarantees that the results are reproducible a The standard Fe oxidation state used in the models for PM (Fe 2+ / Fe tot = 0.85; see Sect. "Parental melt composition"). For testing the effects of varying Fe oxidation state, Fe 2+ /Fe tot values of 0.8 (Fe 2 O 3 = 2.17 wt%, FeO = 7.81 wt%) and 0.9 (Fe 2 O 3 = 1.08 wt%, FeO = 8.78 wt%) were also used in some FC models (Online Resources 1-2) b The standard H 2 O content used in the models for PM (see Sect. "Parental melt composition"). For testing the effects of varying H 2 O, values of 0.4 wt% and 2 wt% were also used in some FC models (Online Resources 1-2). Models ran at 0. Like in the case of the FC-only models, none of the isobaric AFC models explain the key characteristics and the evolution of the low-ε Nd CT1 magmatic system. Adding five recharge events to the AFC model using five aliquots of PM (each 10 mass units) causes negligible geochemical variation relative to AFC-only trend (Online Resource 1). For drastic changes to occur, the recharge magma would have to be of considerably distinct composition relative to PM and other compositions on the AFC trend within a relevant range.

Polybaric two-stage modeling
It is evident from the discussion above and from the model results illustrated in Figs. 6 and 7 that neither FC nor AFC at constant pressure can explain the geochemistry and modal composition of the low-ε Nd CT1 lavas. To weld together the most promising traits of the different models, we performed modeling in two stages: (1) AFC/FC at high pressures (500 MPa; ~ 20 km depth) followed by (2) FC of the variably assimilated magmas at low pressure (0.1 MPa) after the resident magma reached MgO of ~ 9 wt%. In case of the AFC model, this meant assimilation of 6% relative to the mass of the PM during Stage 1. It is important to remind here that PM already represents a contaminated composition (see Sect. "Parental melt composition"), so this number does not correspond to the total amount of crustal contamination experienced by the low-ε Nd CT1 magmas. Based on trace element modeling of a depleted mantle-derived parental melt composition, the total amount of crustal contamination experienced by the most contaminated low-ε Nd CT1 lavas is ≤ 15% (Heinonen et al. 2016).
The results of these two-stage models accompanied by those of a low-pressure FC-only model (Sect. "Fractional crystallization (FC) modeling") are illustrated in Figs. 8 and 9. These two-stage models plausibly replicate the key characteristics of the lavas: (1) AFC/FC at high-pressure stabilizes orthopyroxene in high-Mg melts, and largely explains the variation in major element compositions; (2) Subsequent FC at low-pressure stabilizes plagioclase at ~ 7 wt% of MgO, which explains the inflection point observed in the data (Fig. 8). The total amount of formed cumulates relative to the mass of PM to reach the most evolved lava compositions (MgO ≈ 5 wt%) in all the models presented in Figs. 8 and 9 are in the order of 40-60%. Detailed evaluation of the models on the basis of mineral chemistry is hampered by the lack of representative mineral chemical data on the low-ε Nd CT1 lavas, but it is noteworthy that the earliest olivine and orthopyroxene in the 300-700 MPa models (olivine Fo = 88 mol.%, orthopyroxene Mg# = 89) correspond quite closely to the most primitive compositions measured from the lavas (olivine Fo = 90 mol.%, orthopyroxene Mg# = 84; Luttinen 2000; Heinonen et al. 2018).
Trace element (Zr and Y) and Nd isotopic models for the two-stage scenario are presented in Fig. 9. Zirconium, which is relatively enriched in the wallrock, shows higher variation in the models than Y, which is relatively depleted in the wallrock. Notable enrichment in Zr/Y is only produced by a model that includes assimilation-in FC-only models Zr/Y remains nearly constant throughout the runs; MgO-independent variation in Zr/Y can thus be reasonably explained by variable AFC/FC control during Stage 1 (Fig. 9). Alternatively, the variation in Zr/Y in the lavas may also indicate variable compatibility of Zr during wallrock partial melting (see, e.g., Watson and Harrison 1983) due to wallrock containing zircon (Kreissig et al. 2000). ZrO 2 and zircon are not included in MELTS components and phases, respectively, and thus the D sm (Zr) used for the wallrock in the model is likely to represent a minimum value at least in the early melting steps. Notably, Nd isotopes also exhibit MgO-independent variation that is best explained by variable AFC/FC-evolution of the parental magma, although possibility of variable compatibility of Nd due to allanite in the wallrock (see Kreissig et al. 2000) cannot be completely ruled out either.
In summary, a polybaric two-stage AFC/FC + FC model plausibly explains the mineral identity and major element, trace element, and Nd isotopic composition of the low-ε Nd CT1 lavas. The pressure range for Stage 1 was likely within 300-700 MPa (depth of ~ 10-30 km) that is constrained by early stability of olivine + orthopyroxene as discussed for the previous models. Implications of the best-fit models together with thermophysical considerations are discussed in the following sections.

Thermophysical considerations
To further estimate the plausibility of the two-stage model presented in the previous section, it is possible to approximate the thermophysics of assimilation based on a model of melt transport through cracks. The transport of magma via a crack network has been analyzed by many authors (e.g., Weertman 1971;Spera 1980;Shaw 1980;Emerman et al. 1986;Bruce and Huppert 1990;Spera 1992;Takada 1994;Rubin 1995). Flood basalts have generally been extruded during relatively brief episodes throughout geologic history and are a first-order petrotectonic association. Such igneous events characterized by high volumetric discharge and total volume imply rapid transport of mantle-generated magma in vertical or sub-vertical crack networks created by slightly over-pressurized magma. Typical magma-fracture dike widths are from one to several tens of meters and can extend for hundreds of meters to hundreds of kilometers along strike. Because volumetric discharge rates are high and sustained, vertical (or sub-vertical) heat advection with associated lateral heat conduction into host wallrock can trigger partial melting there and hence, afford the opportunity for basalt contamination by (generally) more silicic partial melts. Here, we use an approximate model that illustrates the dominant effect of local pre-intrusion crustal temperature on the proclivity of magma contamination by wallrock partial melting and concomitant assimilation of partial melts by ascending basaltic magma. Although one should not minimize the uncertainties associated with lack of information regarding the lithologic structure and geothermal gradient of the crust through which the magmas rose, it is useful to frame the thermodynamic evidence for assimilative contamination of low-ε Nd CT1 magmas by Archean crust in the context of magma transport.
According to the presented MCS models, deep AFC contrasts with FC at shallower depths. The essence of the problem for application to low-ε Nd CT1 magmatism is to demark in parameter space the conditions where the wallrock-magma contact migrates either into wallrock (i.e., wallrock meltback, whereby magma is contaminated) or into the magma (frozen magma near wallrock-magma contact chemically insulates ascending magma from contamination). Here, we examine flow in a dike from a thermal perspective once a fracture filled with magma has been generated, perhaps by the magma fracture mechanism discussed elsewhere (Shaw 1980;Spera 1992;Rubin 1995;Detournay 2016).
The key aspect for contamination of mantle-derived magmas by Archean crust is the crustal depth where assimilation occurs. Based purely on the geochemistry, we have argued that primitive mantle-derived magmas assimilated crust at depths of ~ 10-30 km (~ 300-700 MPa) accompanied by olivine and orthopyroxene crystallization. A plausible mechanism for polybaric AFC is by conduit wall partial melting and incorporation of such melts into ascending basalt accompanied by mafic phase crystallization; this is explored semi-quantitatively below.
The model presented here assumes that at time t = 0, a planar crack initially of uniform width (w i ) and length (l) instantaneously forms, and magma flow is initiated along the crack. Thermal boundary layers develop along the walls of the crack, in both marginal wallrock and within the magma conduit per se. Three temperatures are important: T ∞ is the far-field temperature of ambient crustal temperature, T w is the temperature at which rheological blocking occurs in either wallrock or magma (assumed equal for simplicity), and T m is the liquidus (or near liquidus) magma temperature. The rheological blocking temperature is the temperature at which a material makes the solid to fluid transition. In twophase silicate melt + crystal systems, the blocking temperature is defined roughly when the mixture is 50% solid by volume (Lesher and Spera 2015). Heat is transported across the crack by conduction and by advection along the crack parallel to its walls. The migration velocity of the interface between solid and mush at T < T w and melt at T > T w is proportional to the conductive heat flux across the wall and depends on the latent heat and isobaric heat capacities of magma and wallrock. The magma flow rate is driven by a fixed pressure excess (Δp, above the local mean lithostatic pressure) and is associated with a volumetric flow rate (discharge, m 3 /s) Q t related by Poiseuille's relation to the crack width (w t ), l, and magma dynamic viscosity (η).
Solution to the non-linear coupled energy, continuity, and momentum equations including relevant transient terms is formidable as noted by Bruce (1989) and requires numerical solutions that are difficult to generalize for different parameters. However, Bruce and Huppert (1990) present approximate solutions based on a reasonable set of assumptions for the limiting cases when conductive heat transport dominates over advection and vice versa. The latter case, where vertical advection of heat cannot be neglected is most relevant to CFB petrogenesis at Vestfjella and is the case pursued here. In particular, the transition between crack closing by freezing (lateral heat loss dominates over vertical heat advection) and wallrock meltback leading to magma contamination (advection balances conduction in steady state) enables one to define an approximate critical temperature difference between magma and far-field country rock delineating the transition between the processes. There are three dimensionless parameters that govern the critical transition between crack closing and crack meltback. These are the Stefan numbers for wallrock and magma and a dimensionless shear flow term and are defined accordingly: St m = Δh f,m /c p (T m − T w ), St wr = Δh f,wr /c p (T w − T ∞ ), and Π = (ηκl 2 /Δpw i 4 ) 1/3 . The Stefan numbers gage the role of latent heat relative to sensible heat in wallrock heating and magma cooling, where Δh is the enthalpy of fusion of magma (Δh m ) or wallrock (Δh wr ) and c p is the specific isobaric heat capacity (assumed identical for wallrock and magma here). The flow term is a ratio of factors impeding flow to those enhancing flow where κ is the thermal diffusivity. Other variables have been defined earlier.
Although we lack detailed knowledge of the parameters that pertain to any given episode of intrusion or volcanism pertinent to the Vestfjella CFBs, it is possible to test the hypothesis that assimilation during Stage 1 of the low-ε Nd CT1 magmas was governed by larger St wr in the hotter deep crust than during magma transport through the upper cooler crust. If values for the parameters that are roughly constant are set, one can determine an approximate value for the initial crustal temperature that delimits the crack-meltback and crack-closing regimes. For example, scaling the solutions in Bruce and Huppert (1990) and adopting the reasonable constant values Δh f,m = 600 kJ/kg, Δh f,wr = 320 kJ/kg, c p = 950 J/kg K, η = 1000 Pa s, κ = 5 × 10 −7 m 2 /s, l = 10 km, Δp = 10 MPa, w i = 1.5 m, T m = 1250 °C, and T w = 1075 °C (Stebbins et al. 1984;Lange and Carmichael 1990;Rubin 1995;Lange 1997;Petcovic and Dufek 2005;Giordano et al. 2008;Pertermann et al. 2008;Hofmeister et al. 2009;Tikunoff and Spera 2014;Lesher and Spera 2015;Karlstrom et al. 2019; see Online Resource 1 for more details), the transition between crack meltback and crack closing occurs roughly at the initial crustal temperature T ∞ = 280 °C.
In summary, for the parameters listed above, wallrock meltback would occur provided the local country rock temperature was above ~ 300 °C. For a typical pre-intrusion geotherm, this corresponds to depths greater than about 10 km. At shallower levels, by shielding, a selvage of uncontaminated frozen magma would prevent chemical pollution of the upward flowing magma in the interior of the crack. Although this is a rough estimate, it is consistent with the geochemical observations based on the MCS thermodynamic model presented above. Better knowledge of the crustal structure and geotherm as well as the size distribution of the basaltic dikes would provide the necessary information to warrant a more detailed calculation and to test the wallrock meltback hypothesis more thoroughly.

Discussion
The presented modeling provides important constraints on the magma transport system of the low-ε Nd CT1 lavas, and possibly CFBs in general. A two-stage model suggested based on MCS modeling and corroborated by thermophysical considerations achieves the closest match to observed low-ε Nd CT1 lava compositions (Figs. 8,9): (1) AFC/FC at pressures of ~ 300-700 MPa (depths of ~ 10-30 km) followed by (2) FC at shallow pressure (< 100 MPa; see also Fig. 5). In nature, these would translate to (1) AFC/FC affecting rather primitive magmas either in magma bodies or during ascent via a crack network, followed by (2) further differentiation dominated by FC in shallow magma feeder systems (dikes and sills) without significant assimilation (Fig. 10). Near-surface (depth of < 3 km) positive density anomalies revealed by geophysical studies at Vestfjella (Ruotoistenmäki and Lehtimäki 1997) are compatible with Stage 2 shallow feeder systems as suggested by the MCS model.
The rather uniform major element compositions (e.g., MgO = 6.8 ± 0.8 (1σ) wt%, not including the picrite flow; Fig. 4a) of the low-ε Nd CT1 lavas indicate that the crustal magma transport system was replenished by fresh batches of magma and filtered for a rather constant density of the magma (Fig. 4a). The picritic lava flow that contains orthopyroxene phenocrysts obviously bypassed this density filter and ascended towards the surface without forming significant shallow magma chambers (Fig. 10). It is possible that its eruption was initiated by an injection of fresh recharge magma, release of volatiles from the wallrock, or some other process that caused overpressure in the system (see Hartley and Maclennan 2018). Notably, this picritic lava flow already shows signs of crustal assimilation in its trace element signature (Fig. 3), which means that assimilation was taking place very early in the system (see Heinonen et al. 2016). Multiple recharge pulses of primitive magmas may well have heated the wallrock close to its solidus relatively quickly (see previous sections) and before significant Fig. 10 A schematic two-stage model for the crustal architecture of the magmatic plumbing system of the low-ε Nd CT1 basalts of Vestfjella. Stage 1: assimilation and fractional crystallization of early olivine and orthopyroxene take place within the Archean craton at depths of ~ 10-30 km (pressures of 300-700 MPa). amounts of crystallization took place. In addition, the apparent steady-state compositions of the lavas may suggest the involvement of an RTF-type (periodically replenished and tapped, continuously fractionating; O'Hara and Mathews 1981) process that was able to buffer the compositions of magmas, even if they are generated by mixing of very different components (see Luttinen and Furnes 2000).
The gradational stratigraphic variations in contaminationinsensitive incompatible trace element ratios (e.g., Ti/Y: slight increase until ~ 640 m, then slight decrease; Fig. 4b) indicate that over the course of the generation of the low-ε Nd CT1 lava pile, there were slight changes occurring in the parental melt composition (Fig. 4b). These could have been caused, e.g., by variations in the degree of melting or source composition or both. Such processes have been well documented for presently active volcanic systems such as Hawaii (Marske et al. 2008;Greene et al. 2013). In contrast, contamination-sensitive trace element ratios (e.g., Zr/Y, Fig. 4c) exhibit more chaotic behavior in relation to stratigraphy and are more likely be controlled by variable degree of AFC/ FC processes taking place during Stage 1 differentiation or variation in the composition of the assimilant.
Notwithstanding the likely variations in the parental melt and assimilant composition, we suggest that the Al 2 O 3 vs. MgO inflection points in CFB magmas may generally not form under continuous AFC processes, but are governed by fractional crystallization (likely at shallow depths, see Figs. 6,7,8 and 9). On the other hand, existence of orthopyroxene in a CFB sample may hint that at least some of the differentiation took place at high pressures and/or included assimilation of Si-rich crustal materials. We emphasize that these generalized observations always have to be scaled to the specific compositions and magmatic environments of the rocks that are studied. In addition to depth, timing may also be important. Early in the formation of a flood basalt province (such as in the case of the low-ε Nd CT1 lavas that are some of the oldest lavas of Karoo LIP; Luttinen et al. 2015), deeper assimilation is more likely, whereas after up to millions of km 3 of magma have been pumped into the system and the upper crust has been significantly heated, notable shallow contamination can also take place. Such a scenario has been suggested to have taken place during the evolution of the Steens magmatic system of the Columbia River Flood Basalt province (Moore et al. 2018). On the other hand, prolonged melting of the crust after assimilation has begun depletes fusible materials and, in some cases, can inhibit further assimilation (Meade et al. 2014).
It is also relevant to ask if the two-stage differentiation process that took place in the evolution of low-ε Nd CT1 magmas was common in other CFB magmas in Karoo LIP and elsewhere. Neumann et al. (2011) indeed suggested a very similar two-stage crustal differentiation model for the widespread and voluminous North Karoo CFBs. Based on MELTS and EC-AFC modeling, they inferred the first crustal stage to record up to 10% assimilation of lower crustal granulites and the second crustal stage to be dominated by fractional crystallization within the Karoo basin sedimentary rocks. Their and our findings imply that such plumbing systems may be common in other CFB magma types as well. In many studies on other CFB provinces, polybaric AFC/ FC histories are inferred, but without detailed constraints on whether the two processes overlapped or were confined to different depths (e.g., Farnetani et al. 1996;Fram and Lesher 1997;Larsen and Pedersen 2009;Hole 2018). Future studies combining thermodynamical and geochemical approaches hopefully will give more insight into this issue.

Concluding remarks
Modeling with the Magma Chamber Simulator provides a pathway for deciphering the differentiation of a lava sequence (low-ε Nd CT1 magma type) from the Antarctic portion of the Karoo LIP. Fractional crystallization of mantle-derived magmas alone cannot explain the geochemical characteristics of the lavas; assimilation of Archean crust is required. This assimilation must have taken place deep in the crust, where Archean wallrock was available and olivine and orthopyroxene were able to crystallize from the melt (300-700 MPa, depth of ~ 10-30 km). Continuous assimilation and/or crystallization models at these pressures do not stabilize plagioclase at MgO ≈ 7 wt%, however, contrary to what is observed in the lavas. Therefore, a second stage of crystallization at low pressures (≤ 100 MPa) without notable additional assimilation is required. This likely took place in shallow feeder zones, where colder wallrock and/or chilling of earlier near-wall magma shielded the upward flowing magma from wallrock assimilation. Such a two-stage scenario may be feasible for other CFB magmatic systems as well.