Magmatism at Passive Margins: Effects of Depth‐Dependent Wide Rifting and Lithospheric Counterflow

Rifted passive margins exhibit a large variety in the timing, distribution, and amount of magmatism. The factors controlling magmatism during rifted margin formation, remain, however, incompletely understood, partly owing to the complex rifting styles. In this study, we use 2‐Dimensional numerical models to investigate the effects of depth‐dependent wide rifting and lithospheric counterflow on magmatism during rifted margin formation. Results show that a strong crust promotes narrow margins with a sharp transition to normal thickness oceanic crust whereas a weak crust promotes depth‐dependent wide rifting, with preferential removal of mantle lithosphere, leading to formation of wide margins with over‐thickened (>18 km) igneous crust in the distal margin. Counterflow of depleted lithospheric mantle, in contrast, may delay syn‐rift magmatism, and result in exhumed a‐magmatic continental mantle at narrow margins. The combination of wide rifting and lithospheric counterflow results in magma‐poor wide margins, with in some cases a transition to excess magmatic activity at breakup. Our models provide an explanation for the contrasting magmatic productivity at narrow and wide rifted margins, such as observed in the Lofoten‐Vesterålen, Newfoundland, Kwanza Basin, and Orange Basin margins.

2 of 23 et al., 1987), and high-velocity lower crust bodies (LCB) interpreted as magmatic intrusions and underplating at the base of the crust White et al., 2008). In contrast, non-volcanic margins are lacking these characteristics, in some cases exhibiting a broad zone of exhumed mantle with little to no magmatism at the sea floor preceding formation of mature oceanic crust. Examples of non-volcanic margins include the narrow Iberia-Newfoundland margins with exhumed a-magmatic continental mantle preceding the formation of normal oceanic crust (Sutra et al., 2013;Whitmarsh et al., 2001) and ultra-wide central South Atlantic margins (Contrucci et al., 2004;Unternehr et al., 2010).
Magmatism during rifting arises from decompression melting of the sub-lithospheric mantle. Theoretical calculations and numerical models have shown that decompression melting is to first order controlled by the potential temperature of the sub-lithospheric mantle (McKenzie & Bickle, 1988;White & McKenzie, 1989). Higher potential temperature may generate a larger extent of melting and thicker oceanic crust (White, 1997). Given this dependency, it is commonly accepted that the thick igneous crust generated at volcanic margins results from melting of mantle with anomalously high temperatures (McKenzie & Bickle, 1988;Nielsen & Hopper, 2004;White & McKenzie, 1989). Yet, at many margins, such as the Lofoten-Vesterålen (Faleide et al., 2008), Gulf of Aden (Mohriak & Leroy, 2013), Newfoundland (Sutra et al., 2013), Kwanza Basin (Fernandez et al., 2020), and Orange Basin (Blaich et al., 2009) margins, the early oceanic crust has normal thickness in the range of 4-8 km (Figure 1), indicating normal mantle potential temperature. However, the degree of magmatism at these margins can still be highly variable, ranging from the volcanic Orange Basin margins (Figure 1d) to the magma-poor Newfoundland margins ( Figure 1b). It is unlikely that all these margins have abnormally high or low mantle temperature during rifting but suddenly recover to normal temperature at the time of breakup. Therefore, other explanations than variation of mantle temperature are required. Alternative models for high magmatic productivity at volcanic margins include more fertile mantle composition (Brown & Lesher, 2014;, rift history (Armitage et al., 2010), active upwelling (Holbrook et al., 2001;, and/or small-scale convection (Boutilier & Keen, 1999;Kelemen & Holbrook, 1995;Mutter et al., 1988;Simon et al., 2009). For non-volcanic margins, the lack of breakup volcanism has been attributed to low mantle temperature (Reston & Morgan, 2004), slow spreading rate (Bown & White, 1994;Hopper et al., 2004), lithospheric counterflow (Beaumont & Ings, 2012;Huismans & Beaumont, 2011, 2014, or serpentinization (Lavier et al., 2019).
Most previous models of melt generation at rifted margins have focused on variations in composition or temperature of the sub-lithospheric mantle, implicitly assuming simple and uniform break-up of the continental crust and mantle lithosphere. However, observations have shown that many margins exhibit complex, depth-dependent extension (Davis & Kusznir, 2004;Kusznir et al., 2005;Reston et al., 2007;Royden & Keen, 1980). For example, the Vøring margin is characterized by larger mantle stretching factors compared to the crust (Davis & Kusznir, 2004;Kusznir et al., 2005), indicating preferential removal of the lithospheric mantle, while the Galicia margin indicates larger crustal thinning than for the whole lithosphere (Davis & Kusznir, 2004;Reston et al., 2007). Huismans and Beaumont (2011) proposed that such contrasting rifting styles may arise from depth-dependent wide rifting, where the ductile crust extends through distributed deformation and lithospheric mantle breaks earlier than the crust, and from lithospheric counterflow, where the lower part of lithospheric mantle flows in the opposite direction as extension (Beaumont & Ings, 2012;Huismans & Beaumont, 2011, 2014. It is unclear how these contrasting rifting styles with depth-dependent extension interact with melt generation. In our earlier study (Lu & Huismans, 2021), we explored, using thermo-mechanical numerical models, the influence of wide rifting together with mantle potential temperature. Here we build on our previous models to investigate melt production at normal mantle temperature, with a focus on the effects of depth-dependent extension and lithospheric counterflow. We present four end-member models that allow for narrow and wide rifting and for counterflow of depleted mantle lithosphere and show how these influence magmatism during rifted margin formation. Finally we compare our models with observations from natural examples. (1)  (Faleide et al., 2008), (b) Newfoundland (Sutra et al., 2013), (c) South Atlantic central segment (Kwanza Basin) (Fernandez et al., 2020), and (d) South Atlantic southern segment (N. Orange Basin) (Blaich et al., 2009) margins. Thickness of early oceanic crust is annotated with a double arrow. While these margins have normal thickness of mature oceanic crust, they exhibit large variation in magmatic activity.

of 23
where is the spatial coordinate, v velocity, P pressure, η viscosity, ρ density, g gravitational acceleration, T temperature, t time, k thermal conductivity, c p specific heat capacity at constant pressure, A r radioactive heat production per unit volume, ϕ melt fraction, and ΔS is the change of entropy on melting. Repeated indices imply summation (i = x, z). The last term in Equation 3 expresses the absorption of latent heat during melting.
When the stress is below frictional-plastic yield, the rocks deform following temperature dependent power-law viscous flow. The effective viscosity takes the general form: where ̇′ 2 = 1 2̇̇ is the second invariant of the deviatoric strain-rate tensor, n is the power-law exponent, A is the pre-exponential factor, Q is activation energy, V is activation volume, and R is the universal gas constant. The values of A, Q, n and V are derived from laboratory experiments (Table 1). f is viscosity-scaling factor that is used to produce weaker or stronger versions of these materials (Huismans & Beaumont, 2011).
When the stress exceeds the plastic yield, model materials follow a frictional-plastic Drucker-Prager yield criterion: where ′ 2 = 1 2 ′ ′ is the second invariant of the deviatoric stress tensor, C is cohesion, φ eff is the effective internal angle of friction. The presence of pore fluid pressures will decrease the effective internal angle of friction φ such that sin eff = ( − ) sin , where P f is the pore fluid pressure, and is the dry internal angle of friction. For hydrostatic pore pressures, the value of the effective internal angle of friction is eff ∼ 15 • (Huismans & Beaumont, 2014). When rocks are subject to increasing strain, they may lose strength, which is known as strain-softening. Strain softening is applied by linearly decreasing eff( ) from 15° to 2°over a range of accumulated visco-plastic strain of 0.5< ϵ <1.5 (Figure 2d).

Melt Parameterization Model
The approach to predict melt productivity is identical to Simon et al. (2009) and Lu and Huismans (2021), which is further based on (Scott, 1992) and (Nielsen & Hopper, 2004). The solidus temperature is parameterized as a function of depth and depletion: where 0 is the solidus temperature at surface, z is depth, X is the depletion factor which represents the concentration of perfectly compatible element in the solid with an initial value of 1, and ( ∕ ) | and ( ∕ ) | are the gradients of the solidus temperature with depth and depletion, respectively (see Table 2 for parameters). Mass conservation of perfectly compatible element leads to evolution of the depletion factor as: which following Equation 6 implies that the solidus temperature increases with on-going melting. The amount of incremental melt in each time step is calculated as: Our melt parameterization model considers both damp and dry melting assuming the presence of a small amount of water in the upper mantle. Damp melting occurs when the temperature exceeds the wet solidus, which has the same form as the dry solidus but with a reduction of ΔT wet = 200°C, such that sw = − Δ wet . Melting of the host rock will turn to "dry" at some point as water preferentially enters into the melt phase. For an initial amount of ∼810 ppm water in the upper mantle, the transition from damp melting to dry melting ( lim ) is inferred to be ∼1%-2% melt fraction (Hirth & Kohlstedt, 1996). Following (Simon et al., 2009), we assume that the melt fraction in the damp melting regime is linearly parameterized to be 0 on the wet solidus and lim = 0.02 on the dry solidus. Note that melt productivity in the damp melting region is significantly lower than that in the dry melting regime.
It has been suggested that a small amount of melt (<1%) can be retained in the host rock before the melt becomes interconnected and can be extracted to the surface (Forsyth, 1992). When the total melt fraction exceeds the threshold of maximum melt retention of ϕ ret = 1%, the additional volume of melt (ϕ ext ) is "extracted" and converted to the equivalent incremental igneous crustal thickness, which is tracked by a separate set of Lagrangian points that move at surface velocities. The melt parameterization model is coupled to the thermo-mechanical model through feedbacks of temperature, density, and viscosity (Appendix A). A benchmark (Lu & Huismans, 2021) with potential temperatures in the range of 1280-1320°C shows that the predicted igneous crustal thickness with T p = 1300°C is in good agreements  (Gleason & Tullis, 1995). b Dislocation Creep law for Wet Olivine (Karato & Wu, 1993). c Values of pre-exponential Constants have been converted to tensor form.

Table 1
Parameters Used in the Thermo-Mechanical Models with global observations of oceanic crustal thickness (Bown & White, 1994;Grevemeyer et al., 2018), which is used as the reference normal temperature in this study.

Model Design
We build the models upon our previous work (Lu & Huismans, 2021), which demonstrated the effect of margin width on magmatism. Here we explore in addition the influence of lithospheric counterflow. Earlier studies have shown that counterflow of lower lithospheric mantle is controlled by the presence of a thick and depleted mantle lithosphere (Beaumont & Ings, 2012;Huismans & Beaumont, 2011, 2014. Our model setup includes two contrasting types of continental lithosphere ( Figure 2):"A" models with a normal thickness mantle lithosphere and "C" models with a thick and depleted mantle lithosphere. (a) Initial setup for "A" model. The initial composition field is laterally uniform including 25 km upper crust (orange), 10 km lower crust (brown), 90 km lithospheric mantle (green), and 475 km sub-lithospheric mantle (yellow). A small weak seed (pink) is imposed in the center to localize deformation. All boundaries are free slip except the top free-surface boundary. The lithosphere is extended at constant velocity of ±V ext /2 on each side, which is balanced by an inflow (V b ) in the sub-lithospheric mantle. Also shown is the initial temperature profile, which has Moho temperature of T m = 550°C and mantle potential temperature T p = 1300°C. (b) Initial setup for "C" model. The C model is identical to A model except that it has a thicker (165 km) and depleted (Δρ = 30 kg/m 3 ) mantle lithosphere, which consists of 90 km upper (dark green) and 75 km lower (light green) mantle lithosphere. The thermal gradient in the lithospheric mantle is accordingly lowered, while the Moho temperature is configured to be the same (T m = 550°C) as A model. (c) Strength envelop for models A1, A2, C1, and C2 for a given strain rate of 10 −14 s −1 . Type "1" and "2" indicate strong (f c = 30) and weak (f c = 0.02) crustal rheologies, respectively. The dashed lines indicate strength envelopes after strain softening. (d) Strain softening parameterization. The effective frictional angle is linearly decreased from 15° to 2° over accumulated visco-plastic strain range of 0.5-1.5. WQ: Wet Quartz; WO: Wet Olivine.
The initial model is 1,200 × 600 km in size and includes a horizontal layered lithosphere with a 35 km thick crust, 90 km lithospheric mantle, and 475 km sub-lithospheric mantle for the A models. The crust is divided into upper crust (25 km) and lower crust (10 km) for visualization purpose and both layers have the same properties. A weak seed is included in the center of the model at the top of the strong mantle lithosphere to localize initial deformation. A Wet Quartz flow law (Gleason & Tullis, 1995) is used for the crust, which is scaled by a crustal viscosity-scaling factor, f c . Two end-member crustal rheologies are used: a strong (f c = 30, Type 1) and a weak (f c = 0.02, Type 2) crust ( Figure 2c). For mantle materials, a Wet Olivine flow law (Karato & Wu, 1993) is used. The lithospheric mantle is scaled by a factor of five. Velocity boundary conditions are all free slip except for the top boundary that is a free surface. Extension is applied at the side boundaries within the lithosphere with constant horizontal velocities ( ext∕2 = 0.75 cm/yr, half rate), while a return velocity V b is applied in the sub-lithospheric mantle to ensure mass balance.
Thermal boundary conditions are specified temperatures at the top (T 0 = 0°C) and bottom (T bot = 1540°), and zero heat flux at side boundaries. The initial temperature is laterally homogeneous and consists of three segments delimited at the Moho and the base of the lithosphere. The sub-lithospheric mantle has a constant potential temperature = 1300°C with an adiabatic thermal gradient of 0.4°C/km, which results in a base lithosphere temperature of = 1350°C for the A models. The initial geotherm within the mantle lithosphere is linear between the Moho, with a temperature of = 550°C, and the base of the lithosphere. The initial temperature in the crust follows a steady state geotherm for uniform crustal heat production = 0.8776 W/m 3 and Moho temperature of 550°C and a mantle heat flux of 20 mW/m 2 . Adiabatic heating and cooling is taken into account.
C models are designed to investigate the effect of lithospheric counterflow on magmatism. Archean, Proterozoic, and Phanerozoic mantle lithosphere is depleted to varying degrees and less dense than asthenosphere by 20-70 kg/m 3 (Griffin et al., 2009). It remains stable over geological time scales as part of a horizontally stratified lithosphere and is not subject to convective removal owing to its lower density and higher viscosity. Huismans and Beaumont (2011) have suggested that partly refertilized (metasomatized) depleted lower mantle lithosphere may flow toward the rift axis during passive margin formation. The setup of C models is similar to the A models, except that the lithospheric mantle is thick (165 km) and depleted and is divided into 90 km thick upper lithospheric mantle and 75 km lower lithospheric mantle. The upper and lower lithospheric mantle follow a Wet Olivine flow law (Karato & Wu, 1993), and are scaled by viscosity-scaling factors of f ulm = 5 and f llm = 3, respectively. The depleted mantle lithosphere is 30 kg/m 3 lighter at T 0 . The potential temperature and the Moho temperature are the same as in the A models, which results in a base lithosphere temperature of T l = 1380°C and a lower thermal gradient in the mantle lithosphere. This setup ensures that the temperatures in the crust and in Dry-to-wet viscosity ratio Ω -5 Table 2 Melt Parameterization Model Parameters the sub-lithospheric mantle in the C models are the same as in A models. The depleted lithosphere has a higher solidus temperature such that it is not subject to melting in our models.

Results
We investigate melt generation during rifted margin formation by taking into account the effects of (a) different crustal rheologies and (b) the presence of thick depleted lithospheric mantle. We present four end-member cases (see Movies S1-S8 for model animation): Models A1 ( Figure 3) and A2 (Figure 4) with standard mantle lithosphere, and Model C1 ( Figure 5) and C2 ( Figure 6) with thick depleted mantle lithosphere. "1" and "2" denote strong and weak crustal rheology, respectively.

Model A1: Strong Crust, Standard Mantle Lithosphere
Model A1 (Figure 3), which has a strong crust (f c = 30) coupled to the reference thickness (90 km) mantle lithosphere, promotes localized deformation in the center of the model, leading to formation of deep frictional-plastic faults cutting through the crust into the upper mantle lithosphere, tilted crustal blocks, and uplifted rift flanks in the early stage of the model (Figure 3a). Crustal breakup occurs at 6 Ma after 90 km of extension, soon followed by rupture of the mantle lithosphere after which the model evolves to stable seafloor spreading (Figures 3 and 7a). Early rupture of the continental crust leads to the formation of narrow margins, where the width of the margin, defined from the start of crustal thinning to the most distal edge of continental crust (COB; Lu & Huismans, 2021), is less than 50 km wide.
The model produces normal magmatism during both continental rifting and seafloor spreading ( Figure 3). Predicted thickness of igneous crust increases seaward in the narrow continent-ocean transition (COT) and reaches a peak thickness of 7 km. At the center of the model domain, the location of the mid-ocean ridge spreading (MOR) is indicated by the thinnest melt thickness, which is a consequence of our vertical melt "extraction" scheme and the divergence at the spreading center. During seafloor spreading, the stable thickness of oceanic crust is ∼6 km on average, in the range of global normal oceanic crustal thickness (Bown & White, 1994;Christeson et al., 2019). The predicted melt thickness during spreading is smooth as no secondary small-scale convection occurs. This model can serve as a reference case to investigate the effects of more complex depth-dependent rifting styles.

Model A2: Weak Crust, Standard Mantle Lithosphere
Model A2, which has a weak crust ( = 0.02) and standard thickness (90 km) mantle lithosphere, promotes distributed deformation in the crust and the development of wide margins ( Figure 4). In this model, the upper crust and upper mantle lithosphere are decoupled by the weak middle and lower crust. The rifting process can be divided into two stages. Stage 1 is characterized by localized deformation and early rupture of the mantle lithosphere below moderately extending crust (Figure 4a). Stage 2 is characterized by continuous wide and distributed crustal extension, ending with the final crustal rupture (Figures 4b and 4c). Rupture of the crust occurs at 30 Ma, more than 20 Ma later than that of mantle lithosphere (Figure 7b), resulting in formation of ultra-wide (>400 km) margins. No deep faults or uplift of rift flanks develop in this model.
While mantle temperature is the same as in model A1, melt productivity in model A2 is highly contrasting (Figures 7a and 7b). Early rupture of the strong mantle lithosphere results in upwelling of sub-lithospheric mantle leading to decompression melting beneath the extending crust and continuous syn-rift magmatism between 6 Ma and ∼ 30 Ma (Figures 4b and 4e). Ultra-thick igneous crust is accreted to the distal margin, with the peak thickness of igneous crust exceeding 18 km close to the COB (Figure 4d), much thicker than normal oceanic crustal thickness (Bown & White, 1994;Christeson et al., 2019). After crustal breakup, the enhanced melt thickness evolves to a steady state value of 5.5 km (Figures 4d and 7b), consistent with the normal mantle temperature in the model. We note that the model has a normal mantle potential temperature and does not exhibit secondary convection.    crust and upper mantle lithosphere undergo narrow localized deformation in the center, leading to formation of deep faults and uplift of rift flanks (Figure 5a). The crust and upper part of the mantle lithosphere rupture early at 10 Ma (Figure 7c), leading to the formation narrow crustal margins (Figure 5b). Subsequently the model evolves differently from model A1 owing to the presence of the thick and buoyant mantle lithosphere. The depleted lower mantle lithosphere flows in a ductile manner into the rifted zone (Figure 5b), opposite to the extension of the lithosphere, a characteristic process known as "lithospheric counterflow" (Beaumont & Ings, 2012;Huismans & Beaumont, 2011, 2014. In this model, the lower continental mantle is exposed at sea floor and develops a zone (∼70 km) of exhumed continental mantle on both conjugate margins (Figure 5c). Rupture of exhumed depleted mantle occurs at ∼20 Ma, 10 Ma later than crustal breakup (Figure 7c). Despite the same early crustal breakup and narrow margin formation (total width of 96 km) as in model A1, magmatic activity is highly contrasting. Counterflow of the depleted lower lithospheric mantle leads to later magmatism (Figure 5e). After rupture of the exhumed lower lithospheric mantle at about 20 Ma, the model evolves to a mature mid-ocean ridge spreading system with normal thickness of igneous crust (Figure 5d).

Figure 7.
Melt thickness evolution over time for the four end-member models. White arrows show key events during model evolution. Dashed lines represent surface velocity at which predicted melt thickness is moving. Note in all models that melt thickness moves at horizontal surface velocity after the establishment of mid-ocean ridges. In models A2 and C2, the horizontal velocity in the crust, at which the predicted melt thickness moves, is lower than extension velocity at side boundaries before crust breakup, which is a key characteristic of depth-dependent extension.

Model C2: Weak Crust, Thick Depleted Mantle Lithosphere
Model C2, which has a weak crust (f c = 0.02) and thick (165 km) and depleted (Δρ = 30 kg/m 3 ) mantle lithosphere, promotes distributed deformation in the crust and formation of wide margins underplated by lower lithospheric mantle ( Figure 6). Rifting of this model can be divided into three stages ( Figure 7d). As in model A2, stage 1 is characterized by localized deformation and early rupture of upper mantle lithosphere at 7 Ma concomitant with wide distributed extension of the continental crust above. Stage 2 following rupture of the upper lithospheric mantle is characterized by flow of lower lithospheric mantle into the necking area and ends with rupture of the depleted mantle at around 22 Ma. Stage 3, with continuous extension is characterized by juxtaposing of the sub-lithospheric mantle to the highly extended crust and is marked by the final rupture of the crust at 26 Ma.
The magmatic behavior in this model varies with time. Rifting remains a-magmatic until ∼22 Ma at which time the depleted lower lithospheric mantle ruptures, leading to formation of ultra-wide a-magmatic margins. However, thick (∼10 km) igneous crust is accreted to the distal margin. After crustal breakup, the thickness of igneous crust evolves to steady state value of 5.5 km (Figure 6d), consistent with the normal mantle temperature in the model. The model exhibits contributions of both depth-dependent extension and lithospheric counterflow as described in models A2 and C1. As in model C1, the onset of magmatism is significantly delayed owing to the counterflow of depleted lower lithospheric mantle (Figure 6e). After rupture of the depleted lower mantle lithosphere, the very distal margin exhibits igneous magmatic thickness in excess of normal oceanic crust (Figure 6d), a characteristic effect of depth-dependent extension as in model A2. The maximum value of predicted melt thickness is ∼8 km, larger than normal oceanic crust (e.g., model A1), but significantly lower compared the peak thickness in model A2.

The Role of Depth-Dependent Extension in Enhancing Magmatism
Our results demonstrate that wide rifting may lead to increased magmatic outputs. Such an enhancement is explained by depth-dependent extension during rifting, which is characterized by preferential removal of mantle lithosphere and subsequent melting beneath extending crust. As distributed extension in the crust occurs over a much wider horizontal length scale than corner flow mantle upwelling, the horizontal velocity at which the crust moves is significantly lower compared to the rate of mantle upwelling below. Therefore, the crust collects more melt as it stays longer above the area of mantle melting.
The effect of narrow versus wide rifting is illustrated by comparing total melt fraction (depletion) fields between models A1 and A2 at the same time (Figures 8a and 8b). The comparison shows that the total melt fraction fields for the two models share a number of similarities: (a) Materials that experienced significant amount of melting (ϕ > 0.02) are bounded by ruptured lithospheric mantle; (b) these melted materials have similar lower boundaries at 60-70 km in depth; (c) the total melt fraction fields are mostly horizontally homogeneous. These characteristics suggest that the total volume of magmatism in the system is proportional to the distance between the two parts of ruptured mantle lithosphere at the time of crustal breakup, or approximately to margin width defined as the distance between the termination of un-thinned continental crust and the most distal location of continental crust (see Lu and Huismans (2021) for details). This melt comes in excess to what is produced after final lithosphere breakup and adds therefore igneous crustal thickness.
Earlier studies have advocated that small-scale convection may explain enhanced magmatic activity at volcanic margins without a mantle temperature anomaly (Armitage et al., 2013;Boutilier & Keen, 1999;Keen & Boutilier, 2000;Nielsen & Hopper, 2004;Simon et al., 2009). The main issue with this hypothesis is that models with sufficient low mantle viscosity to allow small-scale convection remain convective and unstable long after breakup, resulting in highly non-uniform magmatic productivity (e.g., Boutilier & Keen, 1999;Simon et al., 2009), inconsistent with observed stable oceanic crustal thickness after breakup. Alternatively, other studies (Holbrook et al., 2001; proposed that active mantle upwelling may provide an enhanced flux of mantle through the melting zone. In this model, however, a more buoyant mantle is required to allow active upwelling, which implicitly assumes the involvement of mantle plumes. Our models, with depth-dependent wide rifting, provide an alternative scenario that not only predicts the magmatic pulse at breakup but also provides a mechanism for previously inferred high rates of mantle flow compared to the crust at volcanic rifted margins without requiring the existence of mantle plumes and/or high mantle temperature. While we have emphasized that depth-dependent wide rifting may systematically lead to increased volume of magmatism at margins for given mantle potential temperature, this is not to be understood that the effect of high mantle temperature is excluded. In our earlier work (Lu & Huismans, 2021), we have illustrated that both margin width and mantle potential temperature are important. At narrow margins, high mantle temperature is required to explain excess magmatism.
It is also worth to note that the effect of increased magmatism with margin width does not exclude feedbacks from magmatism, which are not included in this study. For example, melt intrusion may weaken the lithosphere and ultimately limit the width of margins (Bastow et al., 2010;Buck, 2006;Keir et al., 2006;Lang et al., 2020). Such weakening may be particularly efficient at high mantle temperature when the melt production rate is high. The final architecture of rifted margins is expected to be a consequence of both tectonic stretching and magmatic extension. Observed wide volcanic margins in the southern South Atlantic (Chauvet et al., 2021;Clerc et al., 2018) suggest that tectonic stretching might be dominating at these margins.

The Role of Lithospheric Counterflow in Delaying Magmatism
In contract to the effect of margin width that increases magmatic accretion during rifting, lithospheric counterflow may delay melt productivity (e.g., Models C1 and C2). Early melting is inhibited because the depleted mantle that is flowing into the necking area has a significantly higher solidus temperature. The delay time of magmatism, Δt mag , is controlled by the volume of depleted lower lithospheric mantle that flows into the rift area and is proportional to the width of exhumed mantle (Figure 9d). Earlier analysis shows that the amount of lithospheric counterflow scales with the depletion density and inversely with effective viscosity of the mantle lithosphere and with rifting velocity (Beaumont & Ings, 2012). Lower viscosity, higher depletion, or thicker lower mantle lithosphere lead to more efficient counterflow and hence a larger amount of depleted mantle in the necking area. Beaumont (2011, 2014) showed that mantle lithosphere with 15 kg/m 3 depletion and a "Wet Olivine" rheology scaled by a factor of 3 allows for the development of lithospheric counterflow, explaining underplating of wide rifted margins such as in the non-volcanic central South Atlantic with depleted lower mantle lithosphere. In a later study, Beaumont and Ings (2012) show that when using a stronger "Dry Olivine" rheology a higher depletion of 50-80 kg/m 3 is required for significant counterflow. In this study, we use the same "Wet Olivine" rheology scaled by a factor of 3 as in Huismans and Beaumont (2011) and a depletion of 30 kg/m 3 that is higher than that of Huismans and Beaumont (2011). However, as we have included the density feedback of melting (see Appendix A 5.1), which results in a density decrease of 5-15 kg/m 3 for the sub-lithospheric mantle owing to melting, the effective density contrast between the depleted mantle and residual sub-lithospheric mantle after melting is consistent with (Huismans & Beaumont, 2011). The resulting width (∼160 km) of exhumed mantle of Model C1 (Figure 9b) is comparable with observations at the Iberia-Newfoundland conjugate margins (Sutra et al., 2013). Additional models with lower depletion (model C1-a, Δρ = 25 kg/m 3 ; Figure 9a) or higher depletion (model C1-b, Δρ = 35 kg/m 3 ; Figure 9c) of the mantle lithosphere show that the width of exhumed mantle becomes narrower (120 km) or wider (180 km), respectively. The width of exhumed mantle scales with the depletion of lithospheric mantle (Figure 9d), consistent with previous analysis (Beaumont & Ings, 2012).

Interaction Between Wide Rifting and Lithospheric Counterflow
At wide margins lithospheric counterflow interacts with wide rifting. The final magmatic output depends on the competing opposite effects of wide rifting and counterflow on magmatism ( Figure 10). Wide rifting with preferential removal of the lithospheric mantle allows for upwelling of sub-lithospheric mantle beneath the extending crust and generates enhanced magmatic accretion during rifting. Underplating these wide margins by depleted mantle lithosphere suppresses the upwelling of sub-lithospheric mantle and thus inhibits magmatism. If rupture of the depleted mantle below wide margins occurs before crustal rupture, excess igneous crust may by accreted to the distal part of the margin (model C2). In contrast, if the wide rifted crust ruptures before the underplated depleted mantle lithosphere, it may exhume to the surface and form a-magmatic mantle preceding a mature mid oceanic spreading center (model C2-c). This is illustrated by models with varying margin width and lithospheric mantle depletion (Figure 10).

Summary of Magmatic Behaviors for Different Rifting Styles
In previous sections we have shown how end-member rifting styles control magmatic activity. Crustal rheology is the first-order factor controlling margin width (narrow Type 1 vs. wide Type 2 cases). Strong crust coupled to strong mantle lithosphere promotes formation of narrow margins with early establishment of normal thickness oceanic crust (Figure 11a). A weak crustal rheology promotes formation of wide rifted margins and leads to increased magmatic accretion during rifting owing to preferential removal of mantle lithosphere and subsequent decompression melting (Figure 11b). The presence of thick and depleted mantle lithosphere may allow for lithospheric counterflow, during which the lower part of the lithospheric mantle flows into the rifted zone opposite to the direction of extension. Such replacement of sub-lithospheric mantle by counterflow of depleted mantle lithosphere leads to exhumation of continental mantle with little to no magmatism at narrow margins ( Figure 11c). The combination of wide rifting and lithospheric counterflow may lead to complex magmatic activity. If the depleted mantle is insufficient to fill the whole space created by rupture and advection of the upper mantle lithosphere during wide rifting, the margin may exhibit a-magmatic early rifting followed by magma-rich late breakup (Figure 11d). In contrast, if the wide rifted crust ruptures before the underplated depleted mantle lithosphere, it may exhume to the surface and form a-magmatic mantle preceding mature seafloor spreading at wide margins (Figure 10c).

Comparison to Natural Systems
We next compare our model results with four natural systems (Figure 1): the narrow Lofoten-Vesterålen margin (Faleide et al., 2008), the narrow Newfoundland margin (Sutra et al., 2013), the wide central South Atlantic Kwanza Basin (Clerc et al., 2018;Fernandez et al., 2020), and the wide South Atlantic Orange basin (Blaich et al., 2009;Franke et al., 2010). These margins all have normal thickness oceanic crust in the range of 4-8 km following inception of a mature mid oceanic ridge system and are thus consistent with a normal mantle potential temperature. While rifted margins often exhibit a poly-phase rift evolution (Lovecchio et al., 2018;Sutra et al., 2013), we will focus on the last rifting phase leading to breakup.
The Lofoten-Vesterålen margin is the northern segment of the mid-Norwegian margin, separated from the Vøring segment to the south by the Bivrost Lineament/transfer zone (Faleide et al., 2008). Like other Norwegian margins, the Lofoten margin has a prolonged history of extension, with main rift phases widely held to be Early Triassic, Middle to Late Jurassic, Early Cretaceous, and Late Cretaceous to Paleocene leading to break-up at the beginning of the Eocene Skogseid et al., 2000). The Early Eocene continental breakup was magmarich and formed part of the North Atlantic Igneous Province. The Lofoten-Vesterålen margin is located at the outer limit of the Iceland Plume influence (Breivik et al., 2017). Observations indicate that the margin exhibits characteristics of Type I narrow margins (Figure 1a). The Lofoten-Vesterålen margin is characterized by a narrow shelf and steep slope (Faleide et al., 2008). A thin high-velocity lower-crustal body has been identified at this margin, characterizing it as volcanic margin (Breivik et al., 2017;Tsikalas et al., 2005). Early oceanic crust has normal thickness of ∼6 km (Figure 1a). These characteristics can be explained by rifting of strong crust on top of normal thickness mantle lithosphere (e.g., model A1, Figures 3 and 11a).
The second example is the Iberia-Newfoundland conjugate margin system, which are known as the type example of non-volcanic margins (Péron-Pinvidic & Manatschal, 2009;Sutra et al., 2013;Tucholke et al., 2007). We highlight here a transect along the Newfoundland margin (SCREECH) (Sutra et al., 2013). The conjugate margins are characterized by deep faults in the continental crust, a narrow 'necking zone', and a large tract of a-magmatic exhumed continental mantle between the continental crust and the first normal thickness oceanic crust (Chian et al., 1999;Lau et al., 2006;Péron-Pinvidic & Manatschal, 2009;Van Avendonk et al., 2006;Figure 1b). The width of the exhumed mantle domain ranges between 30 and 120 km along the Iberia margin, while it is less well constrained on the conjugate Newfoundland margin (Péron-Pinvidic & Manatschal, 2009). Following exhumation of a-magmatic mantle a mature mid oceanic ridge spreading center was established at around 112 Ma with an average igneous crustal thickness in the range 4.3-6.0 km (Lau et al., 2006), in line with "normal" oceanic crust. The narrow necking zone, the earlier rupture of continental crust, exhumation of a-magmatic mantle, and subsequent normal thickness of oceanic crust can be explained by rifting of strong crust and counterflow of depleted lower lithospheric mantle (e.g., model C1, Figures 5 and 11c).
The central and southern segments in the South Atlantic are both very wide (>300 km) with highly extended crust, but they exhibit distinctly different magmatic activity (Blaich et al., 2013;Contrucci et al., 2004;Moulin et al., 2005; Figures 1c and 1d). While the central South Atlantic is considered to be magma-poor (Blaich et al., 2011;Péron-Pinvidic et al., 2015;Rowan, 2014;Unternehr et al., 2010;Zalán et al., 2011), the margins of the southern segment of the South Atlantic conjugate margins are characterized by significant volcanic addition (Blaich et al., 2013;Franke, 2013;Paton et al., 2017). The central South Atlantic segment is characterized by the presence of a thick (1-2 km) Aptian salt layer along both the Brazilian and West African margins (Karner & Gambôa, 2007;Mohriak et al., 2008). Below the salt layer is a wide sag basin, infilled by syn-rift sediments (Karner & Gambôa, 2007). The hyper-extended crust beneath the sag basin shows high offset low angle faults, suggesting ductile deformation during rifting (Clerc et al., 2018;Karner & Gambôa, 2007;Moulin et al., 2005;Unternehr et al., 2010). Magmatism is limited during rifting, and most of the margin is magma-poor (Péron- Pinvidic et al., 2015;Rowan, 2014;Unternehr et al., 2010). However, the nature of the most distal margin is less constrained and highly debated. Earlier studies have suggested the occurrence of exhumed mantle (Rowan, (Fernandez et al., 2020) identified breakup volcanism and argue that the Gabon-Angola margin evolved initially as a magma-poor margin and in its late stages into a magma-rich system (Figure 1c). The complex magma-poor to magma-rich transition may be explained by rifting of weak crustal and counterflow of depleted lower lithospheric mantle, where the depleted mantle ruptures slightly earlier than final crustal breakup (e.g., model C2, Figures 6  and 11d).
In contrast, the southern segment of the South Atlantic rifted margins is characterized by wide rifting with significant magmatic addition in the form of series of SDR's and high-velocity lower crust bodies, indicating voluminous magmatism during rifting (Blaich et al., 2013;Franke, 2013;Paton et al., 2017). Magmatism in the southern segment of the South Atlantic is commonly associated with the Tristan da Cunha plume (Blaich et al., 2013;Eldholm et al., 2000;White & McKenzie, 1989). However, some margins in the southern segment, such as the North Orange section (Figure 1d), are located >1500 km away from the Tristan da Cunha hotspot, far beyond the typical 2000 km diameter influence of plume head (White & McKenzie, 1989). Moreover, the initial oceanic crustal thickness of 7 km along the margins (Collier et al., 2009;Taposeea et al., 2017) is in the range of normal oceanic crustal thickness (Bown & White, 1994;Christeson et al., 2019). Based on geochemical analysis, Trumbull et al. (2007) inferred a normal mantle potential temperature of T p = 1300°C at the SPRINGBOK cross section in the Orange basin. We suggest that these conflicting characteristics may be explained by rifting of weak crust with preferential of mantle lithosphere, with normal mantle potential temperature, which produces wide margins with excess magmatism (e.g., model A2, Figures 4 and 11b).

Conclusions
We performed thermo-mechanical models to investigate melt generation during rifted margin formation with a focus on the effects of depth-dependent extension and counterflow of depleted lithospheric mantle. Models show that rifting of strong crust with standard mantle lithosphere promotes narrow margins and normal magmatic output as expected by seafloor spreading. Rifting of weak crust promotes depth-dependent extension and formation wide margins with increased magmatic accretion during rifting. Counterflow of depleted lithospheric mantle, in contrast, may delay syn-rift magmatism, and result in exhumed a-magmatic continental mantle at narrow margins. Combination of wide rifting and lithospheric counterflow may lead to underplating of depleted mantle below stretched continental crust and result in magma-poor wide margins, with varying magmatic activity at the final stage of breakup. We show that our models explain the contrasting magmatic productivity in the Lofoten-Vesterålen, Newfoundland, Kwanza Basin, and Orange Basin margins despite of their normal oceanic crust.

Appendix A: Coupling of Thermo-Mechanical System With Melt Parameterization Model
The melt parameterization model is coupled to the thermo-mechanical model through feedbacks of temperature, density, and viscosity. The latent heat feedback in considered by introducing a source term in the heat equilibrium equation (the last term in Equation 3), which represents the absorption of latent heat during melting.
A small amount of melt retained in the host rock ( ≤ ) is assumed to lead to a density reduction (Δ = ( 0 − ) , where ρ 0 and ρ m are densities for mantle and melt at reference temperature, respectively) and a viscosity reduction (χ m = exp(−aϕ), where a is an empirical constant (Nielsen & Hopper, 2004)).
As melting progresses, incompatible components, such as iron and water, preferentially enter into the melt phase, leading to a density decrease (melting depletion) and viscosity increase (dehydration strengthening) of the residual rocks. The density change due to depletion is parameterized as (Nielsen & Hopper, 2004;Scott, 1992) Δ = 0 − −1 (1 − 1∕ ) , where ρ Xref is the density of residual mantle at reference depletion X ref . The final density of rocks is given by = 0 − (Δ + Δ + Δ ) , where Δ = 0 ( − 0) is the buoyancy from thermal expansion. We use a simple linear relationship (Braun et al., 2000;Nielsen & Hopper, 2004;Simon et al., 2009) for dehydration strengthening factor in the damp melting regime as = 1 + −1 , where Ω = 5 is the dry-to-wet viscosity ratio. The final effective viscosity, which increase with more melting, is computed as

Data Availability Statement
Model inputs and outputs, including plotting scripts and the source code to calculate melt fraction, can be accessed from Pangaea Data Archiving and Publication (https://doi.pangaea.de/10.1594/PANGAEA.905111).