Influences on the development of volcanic and magma-poor morphologies during passive continental rifting Open Access

The end products of continental extension—rifted passive margins—exhibit extreme variations in the distribution of continental crust, patterns of subsidence and/or sedimentation, and magnitude and timing of volcanism accompanying continental breakup. The most widely used passive margin morphology classification scheme defines two end-member morphologies, volcanic and magma poor (e.g., Franke, 2013). Assignment of these definitions is related to the relative timing of magmatic emplacement and continental breakup. Volcanic margins are characterized by voluminous magmatic emplacement prior to the full thinning of continental crust (Coffin and Eldholm, 1993; Eldholm et al., 2000; Holbrook and Kelemen, 1993; Hopper et al., 2003; Mjelde et al., 1998, 2002, 1997; Mutter et al., 1982; White and McKenzie, 1989), while magma-poor passive margins express negligible volcanism prior to final continental breakup (Boillot and Froitzheim, 2001; Dean et al., 2000; Whitmarsh et al., 2001).

It is important to note that the interpretation of either end-member morphology does not necessarily diagnose the driver of continental rifting (i.e., active versus passive; Şengör and Burke, 1978). While active, plume-induced continental rifting likely precludes the development of magma-poor margins, externally driven, passive continental extension appears capable of generating margins of either morphologic character (Holbrook et al., 1994; Hopper et al., 1992; van Wijk et al., 2001). Because the extensive partial melting of hot, upwelling asthenospheric mantle is necessary to generate magma during passive continental extension, the distribution of strain throughout the lithosphere may exert a dominant control on the timing and extent of rift-related magmatism. Accordingly, understanding the relationship between tectonic and magmatic processes is essential for unraveling the factors that control the development of passive continental rift systems and understanding the varying observed morphologies of passive margins (Pérez-Gussinyé et al., 2001; Pérez-Gussinyé and Reston, 2001). To provide insight into these influencing factors, we have performed two-dimensional thermomechanical numerical modeling experiments. These experiments systematically explore potential controlling factors and offer constrained estimates concerning the timing and magnitude of melt produced during passive continental rifting.

Volcanic passive margins are characterized by significant magmatic emplacement preceding and/or synchronous with the start of continental rifting (e.g., Geoffroy, 2005; Holbrook and Kelemen, 1993; Mutter et al., 1982). This pre-breakup volcanism indicates the relatively early development of mature magmatic systems that are capable of generating and emplacing significant volumes of melt. These magmatic systems may manifest as onshore igneous emplacements, subaerial seaward-dipping reflector sequences, and/or magmatic underplating of the continental crust (Coffin and Eldholm, 1994; Eldholm and Grue, 1994; Mutter et al., 1985, 1982; Planke et al., 2000; Talwani and Abreu, 2000; White et al., 1987). Globally, the majority of passive margins exhibit some form of volcanic morphology (Menzies et al., 2002; Skogseid, 2001), including large portions of the northern, central, and southern Atlantic Ocean, the southern Red Sea, and nearly the entirety of the Indian Ocean (Coffin and Eldholm, 1992; Mahoney and Coffin, 1997; Mutter et al., 1985; Planke et al., 2000). This distribution suggests that the processes responsible for the formation of volcanic passive margins are relatively common and contribute significant volumes to the global igneous activity budget (Coffin and Eldholm, 1994).

Magma-poor margins are less common and are characterized by an apparent absence of significant volcanism prior to the full thinning of the continental crust (Boillot and Froitzheim, 2001; Dean et al., 2000; Whitmarsh et al., 2001). Because melt generation is delayed relative to continental breakup, post-breakup extension is accommodated through the exhumation of the continental lithospheric mantle, rather than through the emplacement and formation of oceanic crust (e.g., Lavier and Manatschal, 2006; Manatschal, 2004). Domains of exhumed lithospheric mantle and magma-poor morphologies have been observed or interpreted along the margins of Iberia-Newfoundland (Boillot et al., 1995; Hopper et al., 2007; Manatschal and Bernoulli, 1999; Péron-Pinvidic and Manatschal, 2009; Reston, 2007; Tucholke and Sibuet, 2007; Van Avendonk et al., 2009; Whitmarsh et al., 2001), Brazil-Angola (Aslanian et al., 2009; Contreras et al., 2010; Contrucci et al., 2004; Mohriak et al., 1990, 2008), Southern Australia–East Antarctica (Direen et al., 2007, 2011; Espurt et al., 2012; Gillard et al., 2015, 2016), the bight of East India (Bastia et al., 2010; Nemčok et al., 2012; Radhakrishna et al., 2012), and the South China Sea (Hayes and Nissen, 2005; Lester et al., 2014; McIntosh et al., 2014; Savva et al., 2013; Yan et al., 2006; Zhou et al., 1995; Zhou and Yao, 2009). This global distribution suggests that development of a magma-poor margin is not a local phenomenon, but a relatively common result of passive continental extension.

A number of possible controls have been suggested to explain observed differences in continental rift style and the varying morphologies of passive margins. The manner and relative significance of how these suggested controls actually impact tectonic processes and passive margin formation are still a subject of ongoing debate (e.g., Armitage et al., 2010; Brune, 2016; Huismans and Beaumont, 2014; Karner et al., 2007; Svartman Dias et al., 2015). Utilizing simplified kinematic models of uniform, instantaneous extension, McKenzie and Bickle (1988) and White and McKenzie (1989) suggested that mantle potential temperature acts as the first-order control on the magnitude of syn-rift volcanism during continental breakup. Following the addition of finite extension rates to these kinematic models, Bown and White (1995) suggested that extension rates may instead act as the controlling factor and that diffusive heat loss may temper the influence of mantle potential temperature. Similarly, using dynamic numerical models of continental rifting, with serpentinization and mantle melting, Pérez-Gussinyé et al. (2006) suggested that slow extension rates, mantle depletion, and low mantle potential temperature may favor the amagmatic exposure of lithospheric mantle during continental rifting. Armitage et al. (2010) used dynamic numerical models of viscous, decompressive mantle melting to suggest that in the presence of a thermal anomaly, the extensional history of a basin, particularly the existing lithosphere structure, significantly influences the resulting margin morphology. A similar influence of pre-rift lithosphere configuration on resulting margin morphology was noted in the analog and numerical experiments conducted by Corti et al. (2003). Additional factors, including lithosphere thermal structure, initial crustal thickness, and lithosphere rheology, have been suggested by a variety of authors as significant influences on the style of continental rift development and passive margin morphology (e.g., Brune, 2016; Buck, 1991; England, 1983; Huismans et al., 2005; Kusznir and Park, 1987; Svartman Dias et al., 2015). These potential controls are summarized in Table 1 and are investigated within our numerical modeling experiments to assess their relative importance and manner of influence. Within this work, we do not investigate the potential influence of mantle fertility and depletion on rifted margin morphology (e.g., Foulger et al., 2005; Korenaga, 2004). For simplicity with our melt parameterization, the mantle within our models is assumed to be fertile, and as such, melt estimates provided by our experiments should be viewed as upper estimates that may be reduced under depleted mantle conditions.

We utilize an adapted version of the explicit, Fast Lagrangian Analysis of Continua (FLAC) algorithm (Lavier and Manatschal, 2006; Poliakov et al., 1993; Svartman Dias et al., 2015; Tan et al., 2012) to test the manner in which proposed controls (Table 1) influence the relative timing of the first emplacement of melt and the full thinning of the continental crust. FLAC has proven reliable for a number of investigations into continental rifting processes (e.g., Lavier and Manatschal, 2006; Svartman Dias et al., 2015) and was found to be suitable for the adaptations and new melting parameterizations required for this work. FLAC implements brittle, elastoplastic deformation following a Mohr-Coulomb yield criterion (e.g., Lavier and Buck, 2002) and simulates localized faulting using a strain-weakening rule (e.g., Huismans and Beaumont, 2002; Lavier et al., 2000). To simulate ductile deformation, FLAC employs a nonlinear, Maxwell, viscoelastic, constitutive update with viscosity determined via experimental flow laws (e.g., Bürgmann and Dresen, 2008). FLAC’s ability to implement a range of brittle-ductile deformation processes makes it particularly suitable for simulating the varying styles of deformation that occur as the brittle continental crust is thinned and the warm mantle asthenosphere upwells.

To capture important lithosphere and asthenosphere deformational processes, our numerical model covers a domain that is 640 km in width and 300 km in depth (Fig. 1). Our model mesh utilizes variable element spacing in both the vertical and horizontal directions. In the upper 150 km depth, element spacing is 1 km, while in the lower 150 km, the resolution is reduced to every 2 km. Horizontal element spacing is symmetric about the center of the model, with finer resolution at the center and lower resolution at the edges. The outermost 120 km on each end of the model has an element spacing every 3 km, the intermediate 100 km domain every 2 km, and the inner 100 km every 1 km. The initial model setup is laterally homogenous in both physical and thermal properties and applies symmetric boundary conditions. The model is vertically stratified into four layers: upper continental crust, lower continental crust, lithospheric mantle, and asthenospheric mantle. To investigate the influence of crustal thickness on final passive margin morphology, the thickness of the lower crust is varied between 18 and 23 km while the upper crust is held at a constant value of 12 km. The thicknesses of the lithosphere and asthenosphere vary as a function of our thermal structure as discussed below. To investigate the influence of crustal rheology on margin morphology, we vary the rheological properties of the continental crust. The rheological properties of the crust and mantle phases are summarized in Table 2. The differential stress profiles of each varying model are shown in Figure 2. To localize deformation, an initial weak inhomogeneity with minimal cohesion and friction angle is positioned in the center of the model. The inhomogeneity has the same mineral phase and temperature properties as the surrounding material. This inhomogeneity simulates a preexisting weakness, such as a fault, and has a dip of 45° from the surface to a depth of 45 km.

FLAC utilizes a Winkler formulation for the bottom boundary condition to simulate regional isostasy. Upwelling mantle is replaced with asthenosphere of equivalent potential temperature (McKenzie and Bickle, 1988). The surface topography is free. Temperatures at the surface of the model are fixed to 10 °C while the basal temperatures are determined by input model values of mantle potential temperature. To simulate full spreading rates, we symmetrically extend the upper 100 km (the viscous lithosphere domain) on both sides of the model at constant half-rates as prescribed by our input parameter extension rate (full rates: 0.5 cm/yr; 1.0 cm/yr; 2.0 cm/yr). The lower, less-viscous 200 km of our model is extended at lower, linearly decreasing extension rates, with zero extension applied at the base of the model. Adiabatic heating and cooling are included within our model to appropriately capture thermal changes related to changes in pressure. No heat flow is allowed through the sides of our model.

The thermal structure of our numerical model is composed of two parts: a steady-state lithosphere geotherm (Hasterok and Chapman, 2011) and an asthenosphere adiabat of equivalent potential temperature (Fig. 2). Several of the variables investigated in this work affect the formulation of this thermal structure, including chosen input values of mantle potential temperature, surface heat flow, and crustal thickness. An input value of mantle potential temperature (1300 °C; 1350 °C; 1400 °C) is utilized to calculate an asthenosphere adiabat from the surface to 300 km depth (McKenzie and Bickle, 1988). The lithosphere geotherm is steady state and is calculated from the formulation of Hasterok and Chapman (2011). For all models, we use constant values for the thickness of upper continental crust (12 km), surface temperature (11 °C), radiogenic heat production in the lower crust and mantle (0.4 µW/m³ and 0.02 µW/m³ respectively), and thermal conductivity in the crust and mantle (2.3 W/m/K and 3.3 W/m/K respectively). We vary input values for surface heat flow (40.0 mW/m²; 47.5 mW/m²; 55.0 mW/m²) and determine upper crust radiogenic heat production (0.87 µW/m³; 1.03 µW/m³; 1.19 µW/m³ respectively) using the empirical relationship between surface heat flow and radiogenic heat production outlined by Hasterok and Chapman (2011). The parameter space of surface heat flow was chosen based on preliminary test runs where these values were found to capture end-member cold lithosphere and warm lithosphere deformation behavior. While these surface heat flow values are somewhat lower than the global continental average of all environments (64.7 mW/m²; Davies, 2013), they are fairly representative of regions of stable, pre-breakup continental lithosphere (Artemieva, 2006; Hasterok and Chapman, 2011).

To calculate temperature with depth, we employ a one-dimensional bootstrapping method starting with our surface temperature and surface heat flow and using the layered definitions of thermal conductivity and heat production (Hasterok and Chapman, 2011). Once we calculate a temperature greater than or equal to our asthenosphere adiabat, we end our geotherm formulation and set the corresponding depth as the lithosphere-asthenosphere boundary (LAB). Regions deeper than the LAB have temperatures equivalent to the asthenosphere’s adiabat, while regions shallower than the LAB have temperatures prescribed via the lithospheric geotherm. Within the lithospheric geotherm formulation, the chosen value of surface heat flow largely controls the overall lithosphere thermal structure (Hasterok and Chapman, 2011). Low values of surface heat flow correspond to a depressed lithospheric geotherm and a deep LAB, while higher values of surface heat flow correspond to an elevated geotherm and a shallow LAB. Chosen values of mantle potential temperature have a minor effect on the depth of the LAB as they shift the adiabat toward lower or higher temperatures. Temperatures at the Moho and LAB vary as a function of our chosen input values for surface heat flux, crustal thickness, and mantle potential temperature. To aid in comparison with other works, we summarize the Moho and LAB temperatures produced by our inputs in Supplemental Table S1¹. The thermal structure of all the differing models is shown in Figure 2.

In order to offer constrained estimates of decompressive melt production, we utilize the peridotite melting parameterization of Katz et al. (2003). This parameterization assumes batch (equilibrium) melting and offers an estimate of the weight fraction of melt as a function of the temperature and pressure conditions, as well as water and clinopyroxene content. We implement this parameterization using markers that track the advection of material phase, weight fraction melt, and pressure-temperature history. Temperature and pressure conditions are updated in FLAC at each time step. For all of the experiments presented in this paper, mantle phases are anhydrous and contain 15 wt% clinopyroxene. Pressure-dependent functions for the mantle solidus and liquidus follow from Katz et al. (2003). Potential melting conditions are checked within the lithospheric and asthenospheric mantle at each time step. For thermal conditions above the solidus, melt production is solved using a fourth-order Runge-Kutta scheme and assuming melting at constant entropy with a thermal correction for latent heat (Katz et al., 2003). For simplicity, melt and matrix phases advect together and do not separate.

To assess the end-member margin morphology produced by each model, we compare the relative timing of first magmatic emplacement to continental breakup. Because we do not allow melt and matrix phases to segregate, we require an estimate for the first age of magmatic emplacement. To estimate this age, we compare model conditions against melt extraction criteria outlined by Schmeling (2006). Schmeling (2006) contended that if melt fractions were above 0.02 and vertically distributed over a critical thickness of 3–5 km, then melt extraction and magmatic emplacement would occur. We follow these criteria and for each model output (every 0.05 m.y.) search for regions of >2% melt that are vertically connected over thicknesses >3 km. The youngest model age that meets these criteria is assigned as the start of magmatic emplacement. It is important to emphasize that although we use these criteria to estimate when extraction and emplacement begin, we do not actually simulate melt extraction or emplacement within any of the models presented in this paper. We estimate the timing of continental breakup more directly. Continental breakup is assumed to occur when the continental crust is first thinned to <1 km, the minimum resolution of our model. If our criteria for magmatic emplacement are met prior to continental breakup, we assume that magmatic emplacement will occur within the overlying continental crust and result in a passive margin reflecting volcanic end-member morphology. Conversely, if continental breakup occurs first, we assume that post-breakup extension is accommodated by lithospheric mantle exhumation and the resulting passive margin will reflect a magma-poor morphology. This volcanic or magma-poor classification scheme, while simple and binary, highlights the most obvious manner in which the simulated margins may vary. Additionally, this classification is tied to our model’s ability to resolve the relative timing of major events, which we believe is well constrained.

To help tie model results to observations from real-world passive margins, we provide rough estimates for the thickness of igneous crust and the width of various morphologic domains produced in each model. The thickness of igneous crust is estimated at each model output by the calculation of melt thickness. Melt thickness is computed through the integration of new melt produced across the model domain (km²) and division by the applied extension (km) over the time between outputs. Values of melt thickness are rough estimates and should be viewed with a degree of skepticism, particularly at higher values. In real-world systems and in preliminary tests of more complex numerical simulations, the generation of volcanic or oceanic crust via the segregation of melt and mantle matrix phases limits asthenosphere upwelling and tends to equilibrate at melt production conditions of ∼6 km thickness (Bown and White, 1994). We estimate the onset of seafloor spreading as when melt thickness first reaches the 6 km threshold. Depending on model conditions, melt thickness may reach 6 km either before or after continental breakup. If, following continental breakup, melt thickness is only between 0 and 6 km, then we assume that a domain of proto–oceanic crust (<6 km in thickness) is formed. To provide a rough estimate of the domain widths preserved on the conjugate sides of each modeled set of passive margins, the temporal difference between continental breakup, start of magmatic emplacement, and establishment of seafloor spreading can be multiplied by half of the applied extension rate. While useful as a tool to assess the spectrum of margins produced, because melt migration is neglected, the width of this domain should be viewed only as an approximation.

We test a range of potential factors that might influence continental rifting dynamics, the relative timing of magmatic emplacement and continental breakup, and the final margin morphology. These variables include: (1) mantle potential temperature, (2) extension rate, (3) lithosphere thermal structure, (4) crustal thickness, and (5) crustal rheology. Given that a full exploration of this parameter space would require 108 individual model runs, we instead choose a base model (in bold in Table 1) and perform eight experimental runs where we systematically alter only one variable per run. This workflow provides us with nine model runs, from which we can analyze and compare each variable’s relative influence. We recognize that the numerical experiments presented in this paper are unable to cover all potential permutations of continental rifting. However, the systematic exploration of our parameter space can help provide critical insight into the manner in which each variable potentially affects margin morphology, and can help guide future, more-detailed numerical modeling efforts.

We present animations and graphics of all model run results. Model animations display deforming lithosphere and asthenosphere phases, the second invariant of strain, and weight fraction melt. The second invariant of strain highlights deformation and provides insight into the faulting style. For display simplicity, the second invariant of strain is masked in the animations starting immediately prior to the onset of melting. In addition to the animations, we present graphics (Fig. 3–11) displaying melt thickness and minimum crustal thickness from each model run. These graphics compare the relative timing of first magmatic emplacement and continental breakup, and allow for the prediction of either a volcanic or magma-poor morphology. The formulation of melt thickness and criteria for the age of first magmatic emplacement and continental breakup are discussed above in the Methodology section.

The base model is run using an asthenosphere with a potential temperature of 1350 °C, an extension rate of 1 cm/yr, a lithosphere geotherm with a surface heat flow of 47.5 mW/m², a crustal thickness of 35 km, and a dry plagioclase crustal rheology (bold in Table 1). In the animation (^{Animation 1}), we observe crustal deformation progress through several phases and styles. Initial extension exploits our deep, weak inhomogeneity, and through the formation of a conjugate fault, establishes an H-block (Lavier and Manatschal, 2006; Van Avendonk et al., 2009). By 2.50 m.y., the development of a shallow antithetic fault localizes brittle deformation in the upper portions of the continental crust. Ductile deformation in the lower crust begins to reduce the angle of our initial deep-seated faults. At 6.20 m.y., some portions of the crust have thinned to beta factors >2. Starting at 7.00 m.y., multiple synthetic and antithetic faults begin to develop. These faults help to connect the deeper low-angle faults to the upper crust. The first melt is generated 10.30 m.y. after the start of extension. This first melting occurs in upwelling asthenosphere at 2.6 GPa. The pressure and temperature conditions of this initial melt agree with the expected conditions predicted for anhydrous melting of mantle with a potential temperature of 1350 °C (Katz et al., 2003). Utilizing the Schmeling (2006) criteria, magmatic emplacement is estimated to begin at 11.00 m.y. (Fig. 3). Continental breakup occurs later, following the final thinning of superficially exposed lower crust at 15.05 m.y. (Fig. 3). Because initial magmatic emplacement is estimated to occur 4.05 m.y. prior to continental breakup, these experiments suggest that base model conditions would favor volcanic intrusion of the continental crust and formation a passive margin demonstrating a volcanic morphology. We estimate that for the base model, the width of the intruded domain on each conjugate margin would be at least 20.25 km and would have emplaced intrusive phases of up to 11.9 km in thickness at the time of continental breakup.

We present two model runs that help demonstrate the influence of mantle potential temperature (T_p) on margin morphology during passive continental extension. These include a low T_p model corresponding to a potential temperature of 1300 °C (lower than the base model by 50 °C) and a high T_p model corresponding to a potential temperature of 1400 °C (higher than the base model by 50 °C). The low T_p model undergoes crustal deformation that is initially similar to that of the base model (^{Animation 2}). However, by 11 m.y., after the development of more complex fault systems, the rates of crustal thinning between the base and low T_p models begin to diverge. The first melt in the low T_p model is generated at 10.85 m.y. and at pressures of 2.1 GPa (lower pressure due to lower potential temperature). Magmatic emplacement is estimated to begin at 11.55 m.y., with continental breakup following at 13.05 m.y. (Fig. 4). Because initial magmatic emplacement occurs 1.5 m.y. prior to continental breakup, the low T_p model demonstrates a volcanic morphology. Interestingly, even after continental breakup and formation of a volcanic passive margin, the low T_p model has not reached melt production rates capable of sustaining seafloor spreading (Fig. 4). We suggest that during the 1.6 m.y. between continental breakup and establishment of seafloor spreading melt production rates, a domain of proto–oceanic crust would form outboard of the seaward limit of continental crust. We estimate the intruded continental domain to be 7.5 km in width and the proto-oceanic domain to be 8 km in width for each conjugate margin of the low T_p model. Maximum thickness of intrusive crust is estimated to reach 3.2 km.

The high T_p model undergoes crustal deformation in a manner that appears to diverge from that of the base model by as early as 3 m.y. (^{Animation 3}). Similar to the base model, the high T_p model develops shallow antithetic faults that localize brittle deformation in the upper crust, however the location of these faults differs from that in the base model. The first melt in the high T_p model is generated at 8.45 m.y. at pressures of 3.25 GPa (higher pressure due to higher potential temperature). Magmatic emplacement is estimated to begin at 9.10 m.y., with continental breakup following at 12.20 m.y. (Fig. 5). Because initial magmatic emplacement occurs 3.10 m.y. prior to continental breakup, these experiments suggest that the high T_p model would form a passive margin demonstrating a volcanic morphology. The width of the intruded continental domain is estimated to be 15.5 km for each conjugate margin and up to 14.5 km in thickness. Based on our binary classification scheme, both the low T_p and high T_p models appear to favor formation of volcanic margins, and therefore mantle potential temperature does not appear to be a significant and independent control on determining end-member passive margin variability. However, model results demonstrate clear differences in the width and thickness of the volcanic domain and suggest that mantle potential temperature provides an important influence on the magmatic characteristics of passive margins.

Two experimental runs were conducted to investigate the influence of fast (2.0 cm/yr) and slow (0.5 cm/yr) extension rates on passive continental rifting. Please note when comparing the crustal thickness graphics for the extension rate experiments (Figs. 6 and 7) that the scales of the x-axes have been altered compared to other model graphics. In the fast model (^{Animation 4}), we observe crustal faulting that resembles a more rapid form of the deformation seen in the base model. By as early as 3 m.y., the crust has thinned to beta factors >2. The first melt is generated at 5.15 m.y., at pressures (2.6 GPa) identical to those of the base model. Magmatic emplacement is estimated to begin at 5.55 m.y., with continental breakup closely following at 5.65 m.y. (Fig. 6). Melt production does not reach seafloor-spreading thicknesses until 7.00 m.y. As suggested by Bown and White (1995), rapid extension reduces time-dependent diffusive heat loss. In only the fast model, we observe a significant portion of continental lithosphere being induced to melt (^{Animation 4}), suggesting that very little lithosphere heat content was lost due to lateral heat diffusion during upwelling. Magmatic emplacement occurs slightly earlier than continental breakup in the fast model, suggesting that these conditions might favor a very limited volcanic end-member margin morphology. The intruded crustal zone for each conjugate margin is estimated at 1.0 km in width for the fast model with intrusive phase thicknesses of only 0.3 km. The 1.85 m.y. between first melt generation and establishment of seafloor spreading is the fastest increase in melt production observed in all models. However, despite this rapid increase in melt productivity, this experiment predicts the formation a zone of proto–oceanic crust outboard of the seaward limit of continental crust. We estimate the conjugate width of the proto-ocean domain to be 13.5 km.

In the slow model (^{Animation 5}), we observe crustal deformation processes that resemble a sluggish version of the base model. The first melt is generated at 19.35 m.y. and at pressures identical to those of the base model (2.6 GPa). As expected for asthenosphere with higher diffusive heat lost, this initial melt is produced within asthenosphere that is furthest from the laterally adjacent cooler lithosphere. Asthenosphere proximal to the lithosphere is not induced to melt until lower pressures are reached. Magmatic emplacement is estimated to begin at 20.70 m.y. Continental breakup occurs much later, at 28.0 m.y. (Fig. 7). The 7.3 m.y. between first magmatic emplacement and continental breakup is the longest time delay of all model runs. Results from this slow model indicate that lower rates of continental extension tend to favor formation of volcanic margins. We estimate the conjugate width of this volcanic domain to be 18.25 km with intruded thicknesses of up to 10.0 km. Both the fast and slow models highlight the strain-rate dependence of the lower crust and its influence on the timing of continental breakup. Both models predict a formation of a volcanic passive margin, suggesting that extension rate alone does independently or fully control margin morphology. However, similar to mantle potential temperature, extension rates clearly demonstrate an important influence on the overall margin magmatic character, with higher rates of extension favoring a narrower intruded volcanic domain and lower values of intrusive thicknesses.

To evaluate the manner in which lithosphere thermal structure influences margin morphology, we conduct and compare two experimental runs using a warmer elevated (surface heat flow of 55 mW/m²) and a colder depressed (surface heat flow 40 mW/m²) geotherm. In the elevated model (^{Animation 6}), initial crustal deformation proceeds in a manner that strongly mimics the base model. Melt production in the elevated model begins at 10.10 m.y. (Fig. 8), slightly earlier than in the base model but at identical pressures (2.6 GPa). Magmatic emplacement is estimated to first occur at 11.20 m.y. Continental breakup occurs later at 16.20 m.y. (Fig. 8). The 5.00 m.y. between first magmatic emplacement and continental breakup is the third longest of all models, and suggests that an elevated lithosphere geotherm would result in a volcanic passive margin. The protracted time between first magmatic emplacement and continental breakup is likely a result of warmer, ductile lower crust accommodating the final stages of crustal thinning over a lengthier period. For the elevated geotherm model, we estimate a volcanic domain of 25.0 km in width, with intrusive thicknesses up 10.9 km, for each conjugate margin.

The colder, depressed geotherm model (^{Animation 7}) demonstrates a style of crustal faulting distinct from that of the base and elevated models. Extension is accommodated on deeply penetrating faults that bound the H-block and couple the crust and mantle. At 5.00 m.y., following significant H-block subsidence, a pair of initially shallow faults begin to form and localize deformation within the block. By 8.00 m.y., these faults appear to have coupled to the upper mantle, and a new, antithetic, upper-crust fault is formed. Crustal thinning then progresses rapidly and reaches continental breakup by 10.70 m.y. (Fig. 9), the third fastest breakup of all model runs. Melt does not start being produced until 12.95 m.y. with magmatic emplacement estimated to first occur at 13.50 m.y. (Fig. 9). Seafloor-spreading levels of melt production are delayed until 15.25 m.y. The timing of the all magmatic processes, including the first melting, magmatic emplacement, and establishment of seafloor spreading, are delayed in the depressed model relative to the base model. This delay in melt production likely result from the deeper LAB of the cold model, which forces the asthenosphere to upwell over a larger distance to intersect the solidus. Numerical model results indicate that during the 2.8 m.y. between continental breakup and first magmatic emplacement, extension is accommodated through exhumation of lithospheric mantle. Therefore, depressed lithosphere geotherm conditions appear to favor the formation of passive margins with magma-poor morphology. We estimate that the exhumed mantle domain would be 14.0 km in width for each of the conjugate margins produced by the depressed geotherm model. Melt thickness estimates (Fig. 9) suggest the development of a 22.75-km-wide proto-oceanic domain outboard of each of these exhumed mantle domains. Given that elevated and depressed geotherms demonstrate a clear ability to influence and generate end-member morphologies, we interpret the lithospheric geotherm to be a significant control on passive continental extension.

To investigate the influence of initial crustal thickness on margin morphology, we compare results from our base model against an experiment with a thin crust of only 30 km (12 km upper; 18 km lower). In the thin crust model (^{Animation 8}), crustal deformation is dissimilar to that of the base model and instead closely resembles that of the depressed geotherm model. With a thinner lower crust, brittle faults penetrate deeply. Similarly to the depressed geotherm model, this acts to couple the crust and mantle and rapidly thin the crust. Continental breakup occurs rapidly (8.85 m.y.; Fig. 10), the second fastest time to breakup. Although first melting (8.70 m.y.; 2.6 GPa) precedes continental breakup, first magmatic emplacement (9.35 m.y.) postdates breakup (Fig. 10). The timing of continental breakup, first melt, magmatic emplacement, and establishment of seafloor spreading are all earlier in the thin crust model than in the base model. This is likely related to the thinner domain of both crust and lithosphere in the model. Model results indicate that during the 0.5 m.y. between continental breakup and first magmatic emplacement, extension is accommodated via lithospheric mantle exhumation, and therefore the thin crust passive margin would demonstrate a magma-poor morphology. We estimate that this exhumed mantle would occupy a domain 2.5 km in width outboard of the limit of continental crust for each conjugate margin. Similar to the depressed geotherm model, the thin crust model predicts a zone of proto–oceanic crust outboard of the exhumed lithospheric mantle. This proto-oceanic domain is estimated to be 13.75 km in width for each conjugate margin. Comparison of the varying morphologies of the thin crust and base models indicates that crustal thickness acts as a significant control on the development of end-member passive margin morphologies.

In order to investigate the potential influence of crustal rheology, we compare results from our base model against those of an experiment in which a weaker dry quartz rheology is used for the continental crust. Crustal deformation observed in the dry quartz model (^{Animation 9}) differs from that observed in the base model. Following initial H-block formation, two sets of antithetic faults develop in the upper crust at 2.5 m.y. While these faults confine brittle deformation in the upper crust, ductile deformation in the lower crust begins to reduce the angle of the H-block bounding faults. Around 5.5 m.y., a new fault forms in the upper crust outside of the H-block and begins to accommodate delocalized deformation. At 8.5 m.y., a new fault forms directly over the upwelling asthenosphere, couples the crust and upper mantle, and helps re-localize crustal deformation. The first melting occurs at 9.15 m.y., and magmatic emplacement follows closely at 9.85 m.y. (Fig. 11). Seafloor-spreading rates of melt production are reached at 12.05 m.y. Continental breakup occurs at 15.25 m.y. (Fig. 11) and is likely assisted by the shallow, rapidly upwelling asthenosphere. Results from the dry quartz run demonstrate the second longest period (5.4 m.y.) of magmatic emplacement prior to continental breakup. The intruded volcanic domain is estimated to be 27.0 km in width with intrusive thicknesses up to 15.8 km for each conjugate margin. As both the stronger dry plagioclase rheology of the base model and the weaker rheology of the dry quartz model demonstrate volcanic morphologies, it does not appear that crustal rheology is a significant influence on formation of end-member passive margin morphologies. However, a stronger rheology does appear to favor a more limited volcanic domain width and thinner intrusive crust.

We estimate passive margin morphology (volcanic versus magma poor) by comparing the relative timing of continental breakup and first magmatic emplacement. Volcanic passive margins are interpreted to form in numerical simulations where magmatic emplacement begins before continental breakup. Magma-poor margins are interpreted to form in models where magmatic emplacement does not begin until after continental breakup. The timing of continental breakup, initial magmatic emplacement, and establishment of seafloor spreading, as well as estimates for the width of the various domains and estimates of intrusive thickness, are summarized for all models in Table 3. In Figure 12, melt thickness and first magmatic emplacement are graphically displayed relative to continental breakup for all model runs. Each investigated variable demonstrated some ability to influence the timing of magmatic emplacement relative to continental breakup and the magmatic character of the preserved passive margins. A spectrum of margin morphologies was produced with varying degrees of magmatism. Cooler lithosphere geotherms, thinner crust, faster extension, lower mantle potential temperature, and stronger crustal rheologies appear to favor formation of passive margins with more magma-poor affinities. Warmer lithosphere geotherms, thicker crust, slower extension, higher mantle potential temperatures, and weaker crust will tend to favor formation of more magmatic, volcanic passive margins.

Similar to real-world observations (Menzies et al., 2002; Skogseid, 2001), the majority (7 out of 9) of our numerical experiments predict the formation of passive margins demonstrating some volcanic characteristics. Two of these volcanic models, the fast extension and low T_p cases, failed to reach melt production rates capable of sustaining seafloor spreading prior to continental breakup, suggesting formation of proto–oceanic crust seaward of the limit of continental crust. Only in two of our model runs, the depressed geotherm and thin crust experiments, did the tectonic and magmatic processes develop in a manner conducive to the formation of a magma-poor morphology. Both of these magma-poor models have volcanic counterparts (i.e., elevated geotherm; base model–thick crust). Because multiple variables appear capable of reproducing both morphologies, we can conclude that no single variable is alone responsible for determining passive margin end-member morphology. However, even though multiple variables may influence passive margin morphology, our experiments indicate that some factors are more influential than others. It is important to emphasize that our experiments only reflect the relative influence of each variable at the conditions modeled. At different conditions for the base model, the relative effectiveness of each investigated variable would likely be modified.

In our experiments, the lithosphere geotherm appears to be the most significant variable in determining passive margin morphology. A warmer, elevated lithosphere geotherm induces earlier magmatic emplacement relative to continental breakup, and tends to favor a volcanic morphology. Conversely, a cooler, depressed lithosphere geotherm delays magmatic emplacement until after breakup and results in a well-defined magma-poor morphology (Fig. 12). The lithosphere geotherm largely controls passive margin morphology through two complementary influences on: (1) the strength of the continental crust and (2) the depth of the LAB. Cooler geotherms would tend to have both a brittle lower crust that favors rapid continent breakup, and an initially deep LAB that would delay the onset of melting by distancing the asthenosphere from the solidus. Conversely, warmer geotherms would have ductile lower crust that could stymie rapid crustal thinning, as well as a shallow LAB that would favor relatively earlier melting.

Crustal thickness appears to be a significant influence as well, with thinner crust favoring a magma-poor morphology and thicker crust favoring a volcanic morphology. Crustal thickness influences margin morphology by significantly altering both the style and rates of crustal thinning. When the crust is thin, brittle faulting is able to penetrate more deeply and induce early continental breakup and mantle exhumation. The variation in margin morphology resulting crustal thickness experiments is not as extensive as the variability observed for the lithosphere geotherm experiments (Fig. 12). This leads us to suggest that, although a significant influence, crustal thickness is likely not as important as the lithosphere geotherm in determining margin morphology during continental extension.

Despite not being able to produce both end-member morphologies, extension rates and mantle potential temperature exert an influence on the relative timing of magmatic emplacement and continental breakup, an important quality in determining the character of passive margins. Rapid extension appears to favor less pre-breakup volcanism, while slower extension appears to enhance the magnitude of pre-breakup volcanism. This result stems from the strain-rate dependence of the ductile lower crust. In our rapidly extended model, the lower crust thins rapidly, leading to continental breakup while the asthenosphere is just beginning to upwell and is still at considerable depth. This deep asthenosphere is less productive and is ultimately responsible forming the proto-oceanic domain following breakup. Conversely, in the slowly extended model, the asthenosphere upwells more quickly than the crust can thin, and ultimately forms one of the most well-defined volcanic margins of all experiments. It seems likely that under slightly altered initial conditions, variation of the extension rates would be able to demonstrate both types of end-member margin morphology. The modeling results from our experiments disagree with the experimental results found by Pérez-Gussinyé et al. (2006). This discrepancy likely stems from different model resolutions and treatments of lithosphere thermal structure, crustal deformation, and mantle melting. The most important difference between these works appears to be the increased thermal conductivity that Pérez-Gussinyé et al. (2006) utilized to simulate hydrothermal circulation. This increase in conductivity leads to mantle melting that is highly time and extension-rate dependent. The high conductivity and low extension rates of their experiments facilitate heat loss by diffusion, which reduces the temperature of the upwelling asthenosphere, and leads to significant reductions in mantle melt production. It seems likely that if we chose to utilize a similar enhanced thermal conductivity within our model, the results of these studies might have been more similar.

Mantle potential temperatures offer a somewhat confounding picture of their influence on margin morphology. As expected, the lower potential temperatures demonstrate shallower melting, delayed magmatic emplacement, and lower rates of melt production. Higher potential temperature conditions demonstrate deeper, earlier melting and elevated melt production rates. Given the shallow depths the asthenosphere needs to reach to melt in the lower potential temperature model, we expect it capable of developing a magma-poor morphology under marginally different model conditions. Although the high potential temperature model demonstrates a more volcanic morphology than the low potential temperature model, the high potential temperature model does not appear to produce a volcanic domain as wide as the base model. This is due to earlier continental breakup in the high potential model (2.85 m.y. earlier). The highly sensitive nature of strain localization is responsible for differences in the location of faulting and timing of continental breakup between the high T_p, low T_p, and base models. Differences in the thermal and density properties of the asthenosphere affect isostatic forces and strain localization within the crust. It is worth noting that the high potential temperature model produces a thicker intrusive package than the base model due to enhanced melt productivity. Variations in mantle potential temperatures demonstrate a clear influence on the depth, timing, and production rates of decompressive mantle melting and are therefore considered an important factor in determining passive margin morphology.

Differences in crustal rheology were able to produce observable changes in the relative timing of magmatic emplacement and continental breakup. Both the stronger dry plagioclase and weaker dry quartz experiments demonstrated timings of magmatic emplacement and continental breakup that indicated a volcanic morphology. Utilization of a dry quartz rheology slightly extended the timing of continental breakup, causing larger thicknesses of volcanic intrusives relative to the base model. This result suggests that weaker rheologies might favor development of more volcanic passive margins. Conversely, strong rheologies, such as dry plagioclase, might assist in the development of magma-poor morphologies. It is likely that under altered model conditions, crustal rheology might be a more significant factor in determining passive margin morphology. However, for the experiments investigated in this work, crustal rheology did not demonstrate a clear ability to promote a diverse range of potential outcomes, and is therefore judged to be the east important variable in these experiments.

We have conducted nine numerical experiments that investigate the influence of (1) lithosphere geotherm, (2) crustal thickness, (3) extension rate, (4) mantle potential temperature, and (5) crustal rheology on passive continental rifting and development of end-member margin morphologies. Experiment results indicate that all of these variables are capable of influencing continental rift development, and a spectrum of magmatic morphologies is produced. However, only two variables were found to independently and significantly alter the predicted end-member passive margin morphology. Out of all of the investigated variables, changes in the lithosphere geotherm demonstrated the most profound influence on passive margin morphology. Warmer geotherms encouraged development of volcanic margins, and cooler geotherms clearly formed magma-poor margins. Crustal thickness was identified as the second-most influential variable. Changes in lower crustal thickness, as small as 5 km, were able to alter the final passive margin morphology. Thin crust appears to promote development of magma-poor margins, while thick crust encourages formation of volcanic margins. Although changes to extension rate and mantle potential temperature did not alter the binary classification of the resulting margins, fast extension and low mantle potential temperatures appear to favor development of magma-poor characteristics. Conversely, slow extension rates and high mantle potential temperatures might favor formation of volcanic morphologies. Crustal rheology was found to be a more minor influence on margin morphology at the conditions investigated in these experiments, but might be more important at other initial and boundary conditions. Our results lead us to conclude that lithospheric thermal structure and crustal thickness are the most likely first-order controls on passive margin morphological development. Extension rate and mantle potential temperatures likely influence this margin development, but alone cannot account for the observed morphologic variations at passive continental margins.

We thank Sascha Brune and Jolante van Wijk for their reviews which greatly improved the quality and clarity of the manuscript. We thank Jean-Arthur Olive for his helpful discussions. This paper is UTIG Contribution Number 3127. A portion of this work was financially supported by the PLATES consortium.