Lithospheric Control of Melt Generation Beneath the Rungwe Volcanic Province, East Africa: Implications for a Plume Source

The Rungwe Volcanic Province (RVP) is a volcanic center in an anomalous region of magma‐assisted rifting positioned within the magma‐poor Western Branch of the East African Rift (EAR). The source of sublithospheric melt for the RVP is enigmatic, particularly since the volcanism is highly localized, unlike the Eastern Branch of the EAR. Some studies suggest the source of sublithospheric melt beneath the RVP arises from thermal perturbations in the upper mantle associated with an offshoot of the African superplume flowing from the SW, while others propose a similar mechanism, but from the Kenyan plume diverted around the Tanzania Craton from the NE. Another possibility is decompression melting from upwelling sublithospheric mantle due to lithospheric modulated convection (LMC) where the lithosphere is thin. The authors test the hypothesis that sublithospheric melt feeding the RVP can be generated from LMC. We develop a 3D thermomechanical model of LMC beneath the RVP and the Malawi Rift and constrain parameters for sublithospheric melt generation due to LMC. We assume a rigid lithosphere and use non‐Newtonian, temperature‐, pressure‐, and porosity‐dependent creep laws of anhydrous peridotite for the sublithospheric mantle. We find a pattern of upwelling from LMC beneath the RVP. The upwelling generates melt only for elevated mantle potential temperatures (Tp), which suggests a heat source possibly from plume material. At elevated Tp, LMC associated decompression melts occurs at a maximum depth of ∼150 km beneath the RVP. We suggest upwelling due to LMC entrains plume materials resulting in melt generation beneath the RVP.

assumed for the base of the lithosphere with an approximate adiabatic increase in temperature below the lithosphere. The model takes into account rheological flow laws that allow for the generation of sublithospheric melts in a continental setting.
This study is part of the EarthCube project BALTO (Brokered Alignment of Long-Tail Observations), which is aimed at developing new, state-of-art cyberinfrastructure that enables brokered access to diverse geoscience datasets. One of the BALTO developments is a new plug-in for the community extensible National Science Foundation (NSF) open-source finite element code ASPECT that permits the user to access data over the internet using web services from any remote server that uses DAP (Data Access Protocol; Gallagher et al., 2004;Neumiller et al., 2020). This study is a use-case of this BALTO cyberinfrastructure, which accesses lithospheric thickness (Fishwick et al., 2010 updated) to constrain LMC and calculate melt generation beneath the RVP and the Malawi Rift.
LMC has a pattern of upwelling beneath the RVP where lithosphere is relatively thin and produces southern asthenospheric flow along the Malawi Rift. The upwelling from LMC is unable to generate melt except when potential temperatures of mantle are elevated. This result suggests that an additional heat source is present and is likely from plume materials because of the high 3 He/ 4 He detected in RVP lavas (Hilton et al., 2011). At elevated mantle potential temperatures, a significant percentage of sublithospheric melt from LMC occurs at depths of ∼130-155 km localized beneath the RVP, consistent with the location and maximum depth (<200 km) of slow P-wave velocity anomalies beneath the RVP (Grijalva et al., 2018;Yu et al., 2020). These results suggest plume material is present, which explains available geochemical and geophysical observations of the RVP and the Malawi Rift. Sublithospheric melt beneath the RVP provides a source for shallow lithospheric intrusions of magma that weaken the lithosphere (e.g., Buck, 2006;Schmeling & Wallner, 2012) thereby enabling magma-assisted rifting in the northern Malawi Rift.

The Malawi Rift
The Malawi Rift, which represents the southern prolongation of the Western Branch of the EAR (Figure 1), is a weakly extended rift (stretching factor of ∼1.54;  that spans ∼900 km from southern Tanzania, through Malawi, to northern Mozambique. The Malawi Rift (Figure 2a) is characterized by asymmetric half grabens bounded by curvilinear border faults with records of deep seismicity suggesting that the border faults extend to the base of the crust (Craig et al., 2011;Ebinger et al., 2019). Indeed, geophysical studies reveal a thick crust (∼38-45 km; Borrego et al., 2018;) and a relatively strong and thick lithosphere beneath the central Malawi Rift (∼115-210 km; Fishwick, 2010 updated;. Geodetic studies suggest that the rift is opening at a surface velocity of 2.2 mm/yr in the north and 1.5 mm/yr in the south due to an eastward movement of the Rovuma Plate away from the Nubian Plate (Figure 1; i.e., Saria et al., 2014;Stamps et al., 2008;Stamps et al., 2020). The rift is largely magma-poor with volcanism limited to the Pliocene-Pleistocene RVP located in the northern tip of the rift (e.g., Fontijn et al., 2012;Furman, 2007). It is possible that magmatism beneath the RVP contributes to the relatively fast spreading rate of the northern segment of the Malawi Rift.

The Rungwe Volcanic Province
The RVP is the southernmost volcanic region in the Western Branch of the EAR (Figure 1), which lies at the northern tip of the Malawi Rift at the intersection of the Rukwa Rift and Usangu Rift, and covers approximately 1,500 km 2 (Figure 2b; Ebinger et al., 1989Ebinger et al., , 1997Fontijn et al., 2012). The RVP comprises 3 large active volcanoes (Ngozi, Rungwe, and Kyejo; Figure 2b) in addition to more than 100 cones and domes (Fontijn et al., 2010;Harkin, 1960). The RVP lies at the nexus of three major border fault systems including the Livingstone fault of the Malawi Rift, Lupa fault of the Rukwa Rift, and the Usangu border faults, all of which have been active in Miocene-Recent times (Figure 2b;Fotijn et al., 2012).
The relationship between the tectonic structures in the region and magmatism remains controversial (Mesko et al., 2014;Roberts et al., 2012). Ebinger et al. (1989Ebinger et al. ( , 1993 suggest that volcanism in the RVP was synchronous with initial faulting and that local variations in the state-of-stress might result from plate bending under the load of the volcanoes. Radiometric studies of samples from the RVP using 40 Ar/ 39 Ar dating techniques suggest that magmatism in the RVP started by 19 Ma (Mesko et al., 2014(Mesko et al., , 2020 and possibly as early as 25 Ma (Roberts et al., 2012). Although the onset of rifting in the Malawi and Rukwa Rifts is poorly known, Ebinger et al. (1993) used 40 Ar/ 39 Ar radiometric dating of samples from the RVP to suggest that faulting along the Livingstone border fault in the northern Malawi Rift started ∼8.6 Ma. Moreover, U-Pb zircon ages of sediments from Lake Rukwa suggest reactivation and renewed subsidence in the Rukwa Rift at 8.7 Ma (Hilbert-Wolf et al., 2017). The ages of the RVP suggest that volcanism may predate the estimated onset of faulting along the Livingstone border fault and reactivation of the Rukwa Rift. Indeed, seismic studies by Grijalva et al. (2018) and Yu et al. (2020) image low velocity zones beneath the RVP, which suggest the thermal weakening of the lithosphere (with associated magmatism). Thermomechanical modeling by Koptev et al. (2018) suggest that the magmatism beneath the RVP may precede rift-related fault development.
The source of the magma beneath the RVP still remains enigmatic. Studies by Grijalva et al. (2018) and Koptev et al. (2018) suggest the presence of a mantle plume material beneath the RVP that generates melt. This hypothesis is supported by geochemical evidence from RVP lavas showing elevated mantle potential temperatures (Rooney et al., 2012) and elevated 3 He/ 4 He isotopic ratios (Hilton et al., 2011). However, Yu et al. (2020) suggest that the melt beneath the RVP arises from decompression melting in response to passive upwelling associated with lithospheric stretching based on relatively shallow low seismic velocity anomalies. Our study tests these hypotheses. NJINJU ET AL.

Methods
We simulate time-dependent LMC that incorporates melt generation in the sublithospheric mantle in a 3D domain using the finite element code ASPECT (Bangerth et al., 2018a(Bangerth et al., , 2018bHeister et al., 2017;Rose et al., 2017) to test the potential role of LMC in sublithospheric melt generation beneath the RVP and the Malawi Rift. Recent studies demonstrate the capabilities of modeling melt generation and magma dynamics in ASPECT (Dannberg et al., 2019;Dannberg & Heister, 2016), however this is the first study that uses present-day lithospheric structure as input to model melt generation associated with LMC.

Governing Equations
We generate LMC beneath the Malawi Rift by solving for the velocity term u in the conservation equation for momentum (Equation 1) and mass (Equation 2) for an incompressible fluid: where ε(u), p, η, ρ, and g is, respectively, the strain rate, total pressure (sum of the lithostatic and dynamic pressure), viscosity, density, and gravitational acceleration. In order to model melt generation, we also simulate changes in temperature caused by heat transfer in the model by solving the energy conservation equation (Equation 3). We apply the extended Boussinesq approximation that includes shear heating, adiabatic heating, and latent heat of melting in the heating model (Christensen & Yuen, 1985): where C p is the specific heat, α is the thermal expansivity, and k is the thermal conductivity. We assume a uniform crustal thickness of 30 km in the model with an average thermal conductivity of k = 2.5 W.m −1 . K −1 (Njinju, Kolawole, et al., 2019) and a uniform density of 2,700 kg/m 3 . For the lithospheric mantle, we assume an average thermal conductivity of 3.5 W.m −1 . K −1 (Burov, 2011;Koptev et al., 2018) and a uniform density of 3,300 kg/m 3 . We assume a thermal conductivity of 4.7 W.m −1 . K −1 for the sublithospheric mantle (Clauser & Huenges, 1995;Dannberg & Heister, 2016). However, the effective density of partially molten rocks (  eff ) in the sublithospheric mantle is given by: where F is the melt fraction (see Section 3.2), the density of solid rock  solid and the density of melt  melt vary with both temperature and pressure as follows: where β is the compressibility coefficient, ρ 0 is the reference density at reference temperature T 0 and reference pressure p 0 . ρ 0 = 3,300 kg/m 3 for solid rocks and ρ 0 = 3,000 kg/m 3 for melts. We do not model phase changes in the mantle transition zone. We normalize the pressure to a surface pressure of p 0 = 10 5 Pa where T 0 = 293 K (Yang et al., 2018). For the Earth's mantle, α = 3 × 10 −5 K −1 , C p = 1250 J.kg −1 . K −1 , and β = 4.2 × 10 −12 Pa −1 . The latent heat consumed during melting is proportional to the melting rate  and the entropy change ∆S. The latent heat of melting is incorporated, with an entropy change of ∆S = −300 J.kg −1 . K −1 (Dannberg & Heister, 2016).

Model Setup
Our model domain has dimensions of 550 × 1,000 × 660 km along latitude, longitude, and depth, respectively, for a spherical chunk geometry ( Figure 3b). However, our regions of interest beneath the RVP and the entire Malawi Rift are distant from the model boundaries so that boundary effects are limited. Based on previous tests of edge-effects on the interior of the model due to different velocity boundary conditions (free slip vs. zero velocity; ; we set the velocities at all the sides of the model to zero which exerts minimal edge-effects on the model interior from the boundaries of the model. We do not impose plate velocities because we are focused on modeling Poiseuille's flow which we call LMC in this case, rather than Couette flow due to plate motion. Moreover, we are using realistic constraints of the lithospheric structure and want to avoid deforming the lithosphere. We refine the entire model domain to a global mesh refinement of 6 such that each element is ∼8 × 15 × 10 km with 17.5 million unknowns computed on 120 cores. The model is comprised of a lithosphere and an underlying sublithospheric mantle that extends to 660 km depth. The lithospheric structure is read from the BALTO site by the BALTO-ASPECT plugin using the web services provided by the BALTO broker. The lithosphere is part of the updated lithospheric structure model of Fishwick (2010) mapped into the 3D domain for the Malawi Rift and surroundings ( Figure 3a). The lithosphere is thinnest beneath the RVP at the northern tip of the Malawi Rift (∼100 km) and also beneath the southern rift segment (∼100-125 km). The lithosphere is thickest beneath the central segment of the Malawi Rift (∼175-200 km). The lithospheric structure produces lateral variations in temperature and pressure such that relatively thin lithosphere has hotter geothermal gradients and lower overburden pressure than the thicker lithosphere, which is relatively colder and exerts a higher lithostatic pressure. The temperature and pressure variations due to variations in the lithospheric thickness lead to lateral density perturbations in the sublithospheric mantle that drive LMC (Equations 4,5).
The initial temperature structure ( Figure 3b) includes an approximation of a conductive geotherm for the lithosphere with a linear temperature gradient from the surface (T 0 = 293 K) to the base of the lithosphere. The base of the lithosphere is an isotherm defined by the mantle potential temperature T p below which NJINJU ET AL.
10.1029/2020JB020728 6 of 18 the temperature increases approximately adiabatically (0.4 K/km) until it reaches the base of the model at 660 km where the temperature is fixed over time, but laterally varying. This initial temperature condition ( Figure 5) results in lateral temperature variations throughout the model with the relatively coolest temperatures beneath thick lithosphere regions. The mantle potential temperature, T p is the temperature that the mantle would have at the surface if ascended along an adiabat without melting (McKenzie & Bickle, 1988). The T p in the RVP from eruptions occurring in the past 10 Ma ranges from ∼1,420 to 1,450°C (∼1,693-1,723 K) based on the geochemical observations of Rooney et al. (2012). We test a range of T p (1,693-1,800 K) and find that T p = 1,800 K produces a geotherm hot enough to generate melt (see Section 3.3.1). The temperature boundary conditions are given by fixed temperatures at the surface and bottom of the model with zero heat flux at the sides of the model. The surface temperature is fixed at 293 K while the temperature at the base of the model varies laterally in temperature but remains fixed with time.

Rheology
Mantle convection is highly dependent on the viscosity. Since we are interested in LMC in the sublithospheric mantle, we impose a strong, uniform viscosity of 10 23 Pa.s for the lithosphere (Figures 4a and 4b).
For the sublithospheric mantle, we use non-Newtonian, temperature-, pressure-and porosity-dependent creep laws of anhydrous peridotite. Unlike the viscosity model of Keller et al. (2013), which is given by the application of an exponential melt-weakening factor to a constant background mantle viscosity of 10 22 Pa.s, we assume the background viscosity of the sublithospheric mantle is governed by composite rheology for dry olivine material parameters (Jadamec & Billen, 2010). The composite rheology (η comp ) is the effective value of the viscosity from dislocation-creep (η disl ) and diffusion-creep (η diff ) flow laws of dry olivine and is given by: where A is the prefactor, n is the stress exponent,   is the square root of the second invariant of the deviatoric strain rate tensor, d is the grain size, m is the grain size exponent, E a is the activation energy, V a is the activation volume, p is pressure, R is the gas constant, and T is the temperature. The values for the parameters A, n, m, E a , and V a are obtained from experimental studies of dry olivine (Hirth & Kohlstedt, 2004, Table 1). The viscosity of the sublithospheric mantle (η sublith-mantle ) with porosity dependence is given by: where the exponential melt-weakening factor is experimental-constrained to 25 ≤ α Φ ≤ 30 (Mei et al., 2002). We use the average value of α Φ = 27.5. The porosity Φ is simply the melt fraction, or the volume fraction, of peridotite that is already molten. The material properties for each layer (lithosphere and sublithospheric mantle) are defined by compositional fields. The compositional fields for the lithosphere include "crust" and "lithospheric mantle," while the sublithospheric mantle is divided into two compositional fields called "porosity" and "peridotite." Partial melt in the model is tracked through the compositional field called "porosity." The viscosity at each quadrature point is calculated from the effective value of the viscosity of the compositional fields weighted by the volume fraction of each composition at the same location (Figures 4a  and 4b; Rajaonarison et al., 2020).
Setting the velocities at the bottom boundary to zero approximates the effect of the high viscosity jump across the transition zone on slowing mantle flow velocities (410-660 km; Ballmer et al., 2007;Rajaonarison et al., 2020).

Partial Melting
For efficient modeling of melt generation in the sublithospheric mantle, we model low extent melting of anhydrous peridotite, particularly lherzolite, which occurs prior to the exhaustion of clinopyroxene. We use the melting parameterization of Katz et al. (2003), which is valid for shallow upper mantle melting beneath NJINJU ET AL.
10.1029/2020JB020728 8 of 18  continental lithosphere at pressures generally less than 13 GPa, which is the pressure at which the solidus temperatures is a maximum, that is, - ). Partial melting in the sublithospheric mantle is highly dependent on the mantle potential temperature, T p ( Figure 5; McKenzie & Bickle, 1988), and will occur if the T p is such that adiabatically ascending mantle intersects the solidus ( Figure 5). The derived melt fraction F (p, T) depends on the lithostatic pressure p (Pa) and temperature T (K) and is given by: where the mantle solidus temperature T solidus and liquidus temperature T liquidus are respectively given by: where A 1 = 1,085.7  C, A 2 = 1.329 × 10 −7°C /Pa, A 3 = −5.1 × 10 −18°C /Pa 2 , B 1 = 1,475.0°C, B 2 = 8.0 × 10 −8°C / Pa, and B 3 = −3.2 × 10 −18°C /Pa 2 .
We model batch melting, which is melting of an upwelling parcel of mantle rocks without instantaneous melt extraction (Asimow & Stolper, 1999;Ribe, 1985). Batch melting assumes that the melt fraction only depends on temperature and pressure and how much melt has already been generated at a given point, but not considering movement of melt in the melting parameterization. Although we do not model melt extraction by two-phase flow, the effect of latent heat release is switched off by setting the melt freezing rate to zero such that the latent heat term in Equation 3 turns off once the melting rate ) becomes negative due to downwelling and resultant cooling in the melting region. We simulate convection and batch melting for 20 Ma to ensure that steady state is achieved.

Tests of Mantle Potential Temperatures and Lithospheric Thickness Variations in Melt Generation
Partial melting in the sublithospheric mantle occurs if the mantle potential temperature, T p , is such that an adiabatically ascending mantle intersects the solidus. A for dry olivine (Becker, 2006).

Table 1 Rheological Parameters for Dry Olivine Used in the Viscosity Flow Law of the Sublithospheric Mantle
Therefore, partial melting in the sublithospheric mantle is highly dependent on T p and lithospheric thickness variations.

Sensitivity Tests of Mantle Potential Temperature in Melt Generation During LMC
We use T p to constrain the temperature at the base of the lithosphere. The T p in the RVP ranges from ∼1,420 to 1,450°C (∼1,693-1,723 K) based on geochemical observations in Rooney et al. (2012). The geotherms for T p = 1,623-1,723 K do not intersect the solidus ( Figure 5) and so there is no decompression melt even after running the model for 20 Ma. For T p = 1,800 K (1,527°C), the geotherm crosses the solidus ( Figure 5) producing decompression melts. We further test a range of T p values that are closer to 1,800 K (1,770, 1,780, 1,790, 1,800 K; Figure 6). For T p = 1,770 K (blue lines, Figure 6), T p = 1,780 K (yellow lines, Figure 6), T p = 1,790 K (green lines, Figure 6), and T p = 1,800 K (red lines, Figure 6), the geotherms cross the solidus producing initial decompression melts (∼1.5%, ∼3.4%, ∼5.8%, and 8.5%, respectively) which decrease rapidly and cease before 2 Ma. The transient behavior of our melting model is likely due to a transient phase of adiabatic cooling while convection reaches steady state. For T p = 1,800 K (red lines, Figure 6), the initial geotherm is hot enough such that as the model evolves to steady state, LMC generates the second stage of decompression melting when the adiabatic gradients within the active convection cell must have adjusted from the initial values.

Sensitivity Test of Lithospheric Thickness Variations in Melt Generation During LMC
Partial melting in the sublithospheric mantle is also highly dependent on the lithospheric thickness. This is because a thick lithosphere serves as a mechanical barrier to adiabatic ascent of hot mantle materials. The thickness of the melt zone and the maximum extent of partial melting are limited by the lithospheric thickness (e.g., McKenzie & O'Nions, 1995). We test the sensitivity of melt generation to lithospheric thickness variations in our model by conducting simulations with varied lithospheric thickness (+10 km and −10 km) based on the model of Fishwick (2010, updated). For T p = 1,800 K, we find that when we increase the lithospheric thickness by 10 km (increase of mechanical barrier to adiabatic ascent of hot mantle materials), no melt is generated due to LMC (green line, Figure 7). However, when we reduce the lithospheric thickness by 10 km (reduction of mechanical barrier), an unrealistically high peak melt fraction (∼12% melt; blue line, Figure 7) is generated from LMC beneath the RVP. For T p = 1,723 K, reducing the lithospheric thickness by 10 km, generates a maximum instantaneous melt fraction of 0.04% (purple line, Figure 7). This test demonstrates that for LMC to generate decompression melt beneath the RVP, the mantle potential temperature must be highly elevated (T p = 1,800 K). And that decompression melt beneath the RVP will occur at lower T p , only if the lithospheric thickness is more than 10 km thinner than the lithospheric thickness model of Fishwick (2010, updated).

Lithospheric Modulated Convection
In our simulation, LMC develops spontaneously from our initial thermal conditions and forms where there is a transition in lithospheric thickness from relatively thick to thin (see Figure 3a). Figures 8a and 8b show flow patterns for T p = 1,800 K at 17 Ma (time during which the flow is steady-state) resulting from our numerical modeling of LMC at 150 and 250 km depth slices, respectively. Sublithospheric mantle upwelling occurs beneath thin lithosphere, while downwelling occurs beneath relatively thick lithosphere. Our results indicate asthenospheric upwelling beneath the RVP driven by LMC. At 150 km depth (Figure 8a), asthenospheric upwelling (∼1 cm/yr) with a diverging (∼2 cm/yr) horizontal flow occurs only beneath the RVP where the lithosphere is thin (∼100-120 km). At 250 km depth (Figure 8b), the asthenospheric upwelling beneath the RVP is faster (∼3 cm/yr). Another zone of weaker upwelling (∼0.5 cm/yr) occurs beneath the southern end of the Malawi Rift where the lithosphere is ∼160 km thick compared to the thicker lithosphere (∼180 km) in the central part of the rift. Our model suggests a southward flow of the upwelling mantle beneath the RVP toward the thick lithosphere in the central part of the Malawi Rift where the asthenospheric flow is characterized by downwelling (Figure 8b). The lithosphere, which is made rigid in the model, is not deforming.

Melt Generation
The time evolution of our melting model for T p = 1,800 K ( Figure 9) reveals two stages of melting. The first stage, which we call "the transient or unsteady melting state," occurs in the first 2 Ma of the model evolution NJINJU ET AL.
10.1029/2020JB020728   (Fishwick, 2010 updated). Blue line shows melt evolution after reducing the lithospheric thickness by 10 km. Green line shows the melt evolution after increasing the lithospheric thickness by 10 km. Although no melt is generated for a normal lithospheric thickness at T p = 1,723 K, reducing the lithospheric thickness by 10 km (purple line) generates a maximum melt fraction of 0.04% at 0 Ma.  Fishwick (2010, updated). Brown profile AA′ in Figure 8a is the profile location for Figure 10. beneath the RVP. The instantaneous (0 Ma) decompression melt (∼8.5% melt) is not due to LMC; rather the melt arises from the initial conditions, which includes relatively thin lithosphere beneath the RVP and a high mantle potential temperature (T p = 1,800 K). The endothermic melting process consumes latent heat and there is adiabatic cooling of upwelling mantle that rises to the melting region. The melting region thus experiences a net heat loss and progressively cools, such that melting sustained by intrinsic density variations decreases rapidly and ceases by 2 Ma (Ballmer et al., 2007). We suggest the second melting phase arises from LMC, which attains steady state at ∼10 Ma after which the adiabatic gradients within the active convection cell have adjusted from the initial values. During this second stage of decompression melting that arises from LMC, the melt fraction increases rapidly from 0 to <1.5% between 12 and 16 Ma and saturates to 1.5% melt above 17 Ma (Figure 9).  (Figure 10c) with upwelling focused beneath the thin lithosphere of the RVP. The similar flow patterns from 10 to 20 Ma suggest that LMC is stable between 10 and 20 Ma. As described above, we suggest that at 10 Ma ( Figure 10a) melt has yet to generate because of transient adiabatic cooling effects that occur while the convection field is reaching steady state. At 17 and 20 Ma (Figures 10b and 10c, respectively), upwelling from LMC has had enough time to transport deeper mantle materials to shallower depths, raising the sublithospheric geotherm above the mantle solidus temperature that leads to decompression melting of up to 1.5% melt fractions.
Depth slices of the melt model at 17 Ma for T p = 1,800 K (Figures 11a-11d) indicates that melt generation due to LMC is restricted to depths of ∼130-155 km beneath the RVP where the lithospheric thickness is < 120 km. The maximum melt fraction occurs at the center of the melting region (∼145 km) with melt fractions reaching ∼1.5% (Figures 11c).

Sources of Deep Melt Beneath the Rungwe Volcanic Province
The most prominent feature in our model is the isolated region of sublithospheric mantle upwelling due to LMC beneath the RVP that is unable to generate decompression melt except when the mantle potential temperature is elevated (T p = 1,800 K) suggesting that there is a heat source from plume material. In support of this hypothesis, a geochemical study of lava and tephra samples from the RVP by Hilton et al. (2011) shows significantly elevated values of helium isotope ratios ( 3 He/ 4 He) of 15 R A (R A = air 3 He/ 4 He) which far exceeds typical upper mantle values. The high 3 He/ 4 He ratios associated with the RVP could be sourced from the primordial mantle in the core-mantle boundary brought to the surface by upwelling mantle plumes (Courtillot et al., 2003). Such plume-like 3 He/ 4 He ratios, suggest that mantle plume material contributes to the magmatism beneath the RVP. Geochemical study of geothermal fluids/gases in the RVP by Barry et al. (2013) suggest mid oceanic ridge-like 3 He/ 4 He and argue that the absence of hydrothermal samples in close proximity to high (plume-like) 3 He/ 4 He phenocryst localities could explain why plume-like 3 He/ 4 He ratios (>10 R A ) are not found in the fluid and gas samples.
For T p = 1,800 K, decompression melt is generated beneath the RVP at depths of ∼130-155 km with maximum melt fractions reaching ∼1.5% (Figures 11c) which is consistent with melt fractions (1%-2% melt) beneath the Baikal Rift (Yang et al., 2018), which is a relatively magma-poor rift similar to the Malawi Rift. The melting region is spatially consistent with a pronounced low velocity anomaly (LVA) beneath the RVP imaged from P-and S-wave tomographic models developed by Grijalva et al. (2018). Similarly, Yu et al. (2020) used data recorded by seismic stations of the Seismic Array for African Rift Initiation experiment to develop a P-wave anisotropic tomography and found a LVA beneath the RVP that is mostly constrained to the sublithospheric mantle above 200 km depth. Grijalva et al. (2018) attribute the LVA beneath the RVP to the flow of warm, superplume mantle material from the southwest, upwelling beneath and around the Bangweulu Craton lithosphere. This result suggests that the LVA beneath the RVP may be a consequence of partial melts generated from plume material. Upwelling from LMC likely entrains the plume material into shallow sublithospheric mantle beneath the RVP to generate melt. The melt that reaches the base of the lithosphere may refreeze, accumulate in a deep or shallow magma reservoir, inject into the lithosphere as dikes, or erupt to create new crust; however the fate of the melt in the lithosphere is beyond the scope of this study.

Implications for Incipient Rifting
Our numerical model of LMC reveals an isolated upwelling beneath the RVP where the lithosphere is thin. This upwelling beneath the RVP is unable to generate decompression melt except when the mantle potential temperature is elevated (T p = 1,800 K). For T p = 1,800 K, decompression meltoccurs beneath the RVP at depths of ∼130-155 km which is supported by the presence of a LVA beneath the RVP that is mostly constrained to the upper mantle above 200 km depth (Grijalva et al., 2018;Yu et al., 2020). The sublithospheric melt may pond beneath the lithosphere and, subsequently, be injected into the mantle lithosphere and crust through preexisting lithospheric structures ( Note. The similarity in the structure of the mantle flow indicating steady-state LMC from 10 to 20 Ma. We include the yellow dotted lines to help visualize entrainment of deep, hot sublithospheric mantle rising to shallower depths beneath the lithosphere. Also, the melting regions in Figure 10b and 10c are zoomed-in for better visibility and shown as insets. used local measurements of Rayleigh-wave phase velocities to invert for shear wave velocities and clearly observed low velocities (<4.3 km/s) beneath the RVP at crust and upper mantle depths that are consistent with the presence of injected magma. The injection of magma into the lithosphere is an important factor in the process of continental rift initiation since magma can greatly reduce the strength of thick lithosphere and facilitates rifting (Bastow et al., 2011;Buck, 2006;Kendall et al., 2005;Kendall & Lithgow-Bertelloni, 2016;Schmeling & Wallner, 2012;Wallner & Schmeling, 2016). Recent seismic tomography models developed for the RVP and the northern Malawi Rift indicate that the lithosphere beneath the Malawi Rift may have been weakened prior to rifting (Accardo et al., 2020;Grijalva et al., 2018;Yu et al., 2020).

Conclusions
Here, we investigate the sources of melt beneath the Rungwe Volcanic Province (RVP) by developing 3D thermomechanical models of LMC and constraining parameters for sublithospheric melt generation due to LMC. We assume a rigid lithosphere, while for the sublithospheric mantle we use non-Newtonian, temperature-, pressure-and porosity-dependent creep laws of anhydrous peridotite. Our LMC simulation is characterized by an isolated upwelling beneath the RVP where the lithosphere is relatively thin. The upwelling from LMC is unable to generate melt except for elevated mantle potential temperatures suggesting an additional heat source likely from plume material as supported by the high 3 He/ 4 He detected in RVP lavas. We, therefore, conclude that sublithospheric melt generated beneath the RVP results from elevated mantle potential temperatures associated with plume material and that sublithospheric mantle upwelling due to LMC enables the entrainment of this plume material to shallower depth beneath the RVP.

Data Availability Statement
The lithospheric thickness file can be accessed from the Hyrax server through the URL http://balto.opendap. org/opendap/lithosphere_thickness/. The mantle flow models are made available through the Open Science Framework repository with doi https://doi.org/10.17605/OSF.IO/NHFU2.