Lithospheric Structure of the Malawi Rift: Implications for Magma‐Poor Rifting Processes

Our understanding of how magma‐poor rifts accommodate strain remains limited largely due to sparse geophysical observations from these rift systems. To better understand the magma‐poor rifting processes, we investigate the lithospheric structure of the Malawi Rift, a segment of the magma‐poor western branch of the East African Rift System. We analyze Bouguer gravity anomalies from the World Gravity Model 2012 using the two‐dimensional (2‐D) radially averaged power‐density spectrum technique and 2‐D forward modeling to estimate the crustal and lithospheric thickness beneath the rift. We find: (1) relatively thin crust (38–40 km) beneath the northern Malawi Rift segment and relatively thick crust (41–45 km) beneath the central and southern segments; (2) thinner lithosphere beneath the surface expression of the entire rift with the thinnest lithosphere (115–125 km) occurring beneath its northern segment; and (3) an approximately E‐W trending belt of thicker lithosphere (180–210 km) beneath the rift's central segment. We then use the lithospheric structure to constrain three‐dimensional numerical models of lithosphere‐asthenosphere interactions, which indicate ~3‐cm/year asthenospheric upwelling beneath the thinner lithosphere. We interpret that magma‐poor rifting is characterized by coupling of crust‐lithospheric mantle extension beneath the rift's isolated magmatic zones and decoupling in the rift's magma‐poor segments. We propose that coupled extension beneath rift's isolated magmatic zones is assisted by lithospheric weakening due to melts from asthenospheric upwelling whereas decoupled extension beneath rift's magma‐poor segments is assisted by concentration of fluids possibly fed from deeper asthenospheric melt that is yet to breach the surface.


Introduction
Continental rifts are elongate, tectonically induced depressions within the Earth's continents. They are underlain by lithosphere that has been modified by extension and typically characterized by volcanism, thinning of the crust and lithospheric mantle, low seismic velocities, and elevated heat flow (Olsen, 1995;Thybo & Nielsen, 2009). Numerical models for continental rift initiation highlight the importance of an actively upwelling asthenosphere followed by extensive magmatic diking that thermally weakens the lithosphere and may reduce the force required to tectonically stretch an initially thick continental lithosphere up to a factor of six (e.g., Buck, 2006;Schmeling, 2010). Numerous studies point to the important role of magma and fluids in rift initiation because tectonic stretching requires forces that are greater than is tectonically available to extend and thin the lithosphere to lead to lithospheric rupture (e.g., Buck, 2006;Stamps et al., 2010Stamps et al., , 2014. Globally, there is ample support for magma-assisted rifting provided by the occurrence of surface magmatism associated with many rift systems such as the eastern branch of the East African Rift System (EARS; Figure 1; Ebinger, 2005;Kendall et al., 2005). However, numerous examples also exist worldwide where rifting is accompanied by limited exposure of volcanic rocks, such as in the Rio Grande Rift (Wilson et al., 2005) and the Malawi Rift, which is the focus of this study and represents the southern segment of the western branch of the EARS (e.g., Ebinger et al., 1989;Fontijn et al., 2012;Heilman et al., 2019).
The EARS represents an ideal laboratory for studying magma-poor rifting because it comprises the two contrasting magmatic styles of rifting with the eastern branch being magma rich and the western branch being magma poor [ Figure 1; Koptev et al., 2015Koptev et al., , 2018 and references therein]. The western branch includes rift segments, such as the Malawi Rift, that are considered well developed in terms of morphological expression yet have only isolated zones of magmatism. Thus, our understanding of the role of magma in the development of magma-poor continental rifts remains unclear, and magma-poor rifts remain enigmatic. This is largely because there are inadequate geoscientific observations to accurately constraint the lithospheric and thermal structures beneath magma-poor continental rifts.
Geoscientific observations from different segments of the EARS suggest that preexisting structure plays an important role in localization of extension during the evolution of continental rift systems at both regional and local scales (e.g., Corti et al., 2007;Katumwehe et al., 2015;Kolawole et al., 2018;Leseane et al., 2015;Van Wijk, 2005). Preexisting structure such as Precambrian suture and shear zones are regions of mechanically weak lithosphere, and they are prone to reactivation (McConnell, 1972) by normal faulting when aligned within 30°of the trend and 10°of the plunge of the orientation of the maximum tensile stress (Morley et al.,  Saria et al. (2014). Black bold lines represent plate boundaries from Stamps et al. (2008). Black thin lines delineate international borders with the names of the main countries transect by the western branch labeled in blue colors. The inset map shows the relative location of the East African Rift System (red rectangle) on Earth.

10.1029/2019TC005549
Tectonics mutually independent stretching of the crust and the lithospheric mantle during the rifting process (Hopper & Buck, 1998), such that stretching in the lithospheric mantle accommodates most of the lithospheric extension.
In this work, we use the two-dimensional (2-D) radially averaged power-density spectrum analysis (henceforth spectral analysis) of Bouguer gravity anomalies of the World Gravity Map 2012 (WGM2012) to determine the lithospheric structure beneath the Malawi Rift and surrounding Precambrian orogenic belts and Paleozoic-Mesozoic Karoo-type rift basins. We then construct 2-D forward models of the gravity data in order to reinforce lithospheric structure results obtained from the spectral analysis. Subsequently, we use numerical experiments (3-D EDC modeling) to model asthenospheric flow induced by variations in the lithospheric thickness beneath the Malawi Rift and quantify the rate of lithospheric extension exerted by asthenospheric upwelling.

Tectonic Setting and Architecture of the Malawi Rift
The Malawi Rift (henceforth MR) (Figures 2a and 2b), which is the southern segment of the western branch of the EARS, extends for~900 km from southern Tanzania, through Malawi, to northern Mozambique. The rift is largely magma poor except for the Pliocene-Pleistocene RVP located in northern tip of the rift (e.g., Ebinger et al., 1989;Fontijn et al., 2012;Heilman et al., 2019). From north to south, the first 550 km extent of the rift is occupied by Lake Malawi, which has a width ranging between 50 and 75 km. In the north, the rift

10.1029/2019TC005549
Tectonics strikes is NNW trending and then continues with a general N-S trend before assuming a NNE trend at its southern end (Figures 2a and 2b).
The MR is bounded by curvilinear border faults, the most prominent of them is the~120-km long W dipping Livingstone border fault located at the eastern side of the northern MR (henceforth NMR) (Figure 2a). The 2009 M w 6.0 Karonga earthquake in the NMR is related to migration of tectonic activity from the Livingstone border fault onto the hanging wall faults (Biggs et al., 2010;Fagereng, 2013;Gaherty et al., 2019) where the earthquake reactivated prerift Precambrian basement structures along the hanging wall hinge zone (Kolawole et al., 2018). The~200-km long Usisya fault system dominates the central MR (henceforth CMR; Figure 2a; Contreras et al., 2000;Laó-Dávila et al., 2015), whereas the~125-km long Bilila-Mtakataka fault dominates the southern segment of the MR (henceforth SMR; Figure 2a; Jackson & Blenkinsop, 1997). The SMR terminates against the NW trending Shire Graben bounded by the Mwanza fault (Figure 2a; Castaing, 1991).
Previous studies of the MR observe a southward decrease in the amount of the basin's subsidence, thickness of the Cenozoic-Quaternary sedimentary rocks, and elevation of the topographic escarpments (Betzler & Ring, 1995;Flannery & Rosendahl, 1990;Specht & Rosendahl, 1989). These observations are used to suggest that the NMR is relatively old than its southern part and provide evidence for southward propagation of the rift and its opening in a "zipper-like" fashion from north to south (Betzler & Ring, 1995;Flannery & Rosendahl, 1990). This zipper-like N-S opening of the MR is supported by geodetic studies, which suggest that the rift is opening at a surface velocity of 2.2 mm/year in the north and 1.5 mm/year in the south due to an eastward movement of the Rovuma Plate away from the Nubian Plate [ Figure 1; i.e., Stamps et al., 2008;Saria et al., 2014]. Based on 40 Ar/ 39 Ar radiometric dating of samples from the RVP, Ebinger et al. (1993) conclude that rifting at the northern tip of the MR started~8.6 Ma, with volcanism continuing to the present. However, recent 40 Ar/ 39 Ar radiometric dating of samples from the RVP by Mesko et al. (2014) place the onset of active volcanisms at~17 Ma.

Prerift Basement Structures of the MR
The MR traverses a complex array of Precambrian orogenic belts and Paleozoic-Mesozoic Karoo-type sedimentary basins (Figures 2a and 2b). Studies by Laó-Dávila et al. (2015) suggest that preexisting lithospheric heterogeneities control a hierarchal segmentation of the MR. The NNW strike of the NMR is parallel to the tectonic fabrics of the Paleoproterozoic Ubendian Belt (Heilman et al., 2019;Laó-Dávila et al., 2015;Ring et al., 2002). In the NMR, the Mughese Shear Zone and other shear zones within the Ubendian Belt control the border and intrabasinal fault development (Daly et al., 1989;Dawson et al., 2018;Delvaux et al., 2012;Ebinger et al., 1987;Heilman et al., 2019;Kolawole et al., 2018;Ring, 1994;Wheeler & Karson, 1989). In the CMR, particularly south of the Mwembeshi Shear Zone and the Macaloge Shear Zone, the rift is at high angle to the prerift Precambrian basement structure, and the border faults are poorly developed [ Figure 2b; Laó-Dávila et al., 2015]. The Mwembeshi Shear Zone represents the boundary between the Irumide Belt in the NW and the southern Irumide Belt to the SE [ Figure 2b; Fritz et al., 2013]. The pre-Mesoproterozoic Niassa Craton is located in the southern Irumide Belt (Andreoli, 1984;Daly et al., 1989;De Waele et al., 2006;Sarafian et al., 2018). The southern Irumide Belt is sutured to the Marrupa-Unango complex by a shear zone that is subparallel to the surface expression of the MR [ Figure 2b; Westerhof et al., 2008]. In the southern termination of the MR, the Precambrian basement is comprised of the NW trending Neoproterozoic Zambezi Belt (Hanson et al., 1994;Hargrove et al., 2003), which represents the westward continuation of the Precambrian Lurio Belt (Sacchi et al., 2000).

Previous Geophysical Studies of the MR
A number of previous studies have been conducted to investigate the crustal and upper mantle structures beneath the MR using controlled-source seismic (lake-bottom seismometers), body and surface wave tomography, and receiver function stacking (Accardo et al., 2017;O'Donnell et al., 2013;Reed et al., 2016). O'Donnell et al. (2013) used Rayleigh wave phase velocity to invert for a 3-D shear wave velocity model and observed pronounced velocity lows beneath the northern and southern ends of the Lake Malawi but not beneath the central portion of the lake. However, the origin of the LVZs is not well understood.

10.1029/2019TC005549
Tectonics Reed et al. (2016) used receiver function stacking to image the 410 and 660-km discontinuity and observed a normal mantle transition zone beneath the entire MR indicating the absence of any thermal anomaly due to a deep hot asthenospheric upwelling at the mantle transition zone. However, in the CMR, Reed et al. (2016) observed an apparently uplifted zone of 20 km within both the 410 and 660-km discontinuities, which they interpreted to be due to the presence of a cold and strong cratonic lithosphere. A surface wave tomography study by Fishwick (2010) found an E-W belt of thick lithosphere (~180-210 km) beneath the CMR and SMR. Reed et al. (2016) interpreted the belt of thick lithosphere to represent the eastward continuation of cratonic lithosphere of the Congo and Kalahari Cratons. Accardo et al. (2017) used lake-bottom seismometers and onshore ambient noise and teleseismic Rayleigh wave phase velocities to determine the crustal and uppermost mantle structures beneath the NMR and found decrease in velocities at longer periods (25-60 s) consistent with the replacement of the colder higher velocity lithospheric mantle by hotter slower asthenospheric mantle. This low-velocity anomaly is to the west of Lake Malawi and persists beneath the Mughese Shear Zone and Nyika Plateau. Accardo et al. (2017) suggest that the presence of localized melting beneath the RVP may be due to thermal or compositional anomalies in the upper mantle, synrift or prerift lithospheric thinning, or a combination of both.

Data and Methods
Magnetotelluric and seismic techniques have been widely used to estimate the Moho and lithosphereasthenosphere boundary (LAB) depths (e.g., Borrego et al., 2018;Sarafian et al., 2018;Wang et al., 2019). However, the magnetotelluric and seismic methods are often used to investigate the lithospheric structure along profiles or in small regions because of the expense and time needed for data collection, thus limiting their application for investigating large regions like the entire MR. The recent availability of global coverage of satellite gravity data allows for the determination of the lithospheric structure over wide areas. Using power-density spectral analysis of satellite gravity data to determine the crustal thickness and depth to LAB offers an affordable and reliable alternative [i.e., Tselentis et al., 1988;Gómez-Ortiz et al., 2005, 2011Emishaw et al., 2017].

Gravity Data
We use Bouguer gravity data extracted from the WGM2012 Bonvalot et al., 2012) to estimate the depth to the crust-mantle discontinuity (Moho) and the LAB beneath the entire MR and surroundings. The WGM2012 model, produced by the International Gravimetric Bureau, takes into account the topography of the Earth, which is approximated as a spherical prism, by computing the global Bouguer gravity data using a spherical terrain correction . The correction involves subtracting the Earth's topography-derived gravity model from the Earth Gravity Model 2008 (EGM2008). EGM2008 comprises surface gravity measurements from land, airborne, and marine surveys and measurements from satellite altimetry and satellite gravimetry (part of the Gravity Recovery and Climate Experiment; Bonvalot et al., 2012). The WGM2012 model has a spatial resolution of~9 km. However, in regions where terrestrial gravity data are sparse, the spatial resolution of the WGM2012 model is more than 9 km. In regions where terrestrial gravity is unavailable, like over Lake Malawi, the spatial resolution is~140 km (spatial resolution of Gravity Recovery and Climate Experiment data; Pavlis et al., 2012), therefore, shorter-wavelength anomalies will be artifacts. The Bouguer gravity anomalies of the WGM2012 model are computed using the EGM2008 data set and spherical harmonic analysis of the Earth Topography 1 arc second data set with a reference density of 2.67 g/cm 3 for the crust . We gridded the Bouguer gravity data of the WGM2012 using the minimum curvature method (Briggs, 1974) with a grid cell size of 1′ × 1′ (~2 × 2 km; Balmino et al., 2012) to produce a Bouguer gravity anomaly map ( Figure 3). We used the minimum curvature method of gridding (Briggs, 1974), which is built into the Oasis Montaj MAGMAP program, because the Bouguer gravity is a continuous field and the minimum curvature surface is the smoothest possible surface that will fit the Bouguer gravity data values (Foss, 2011).

2-D Radially Averaged Power-Density Spectrum Analysis
We upward continue the observed Bouguer gravity data to 2 km to attenuate short-wavelength anomalies that would be <1 km while accentuating deep-source long-wavelength anomalies (Foss, 2011). Gravity surveys derived from dense networks of gravimeters often span regions less than 1 × 1 km, therefore, upward 10.1029/2019TC005549 Tectonics continuing to 2 km would sufficiently filter out potential artifacts generated by the gravity observations from these dense networks. We use the spectral analysis developed by Tselentis et al. (1988) to estimate the depth to the Moho and the LAB beneath the MR. The power spectrum of gravity describes how the power of the gravity response distributes over different wave numbers (Russo & Speed, 1994). Crucial to this work is the concept of radially averaged power spectrum, which is the direction-independent mean spectrum defined as the average of all possible directional power spectra in concentric rings with diameters equal to the lengths of corresponding wave number vector (Maus & Dimri, 1996;Naidu & Mathew, 1998). The radially averaged power spectrum provides a means to condense information contained in 2-D spectra in 1-D (Maus & Dimri, 1996;Naidu & Mathew, 1998;Spector & Grant, 1970), and power spectra in 1-D enable the ensemble of density anomaly parameters (e.g., average thickness and average widths) to be factored out completely for effective analysis.
We calculate the radially averaged power spectrum of the upward continued gravity data for~110 ×~110 km (1°× 1°) overlapping windows using the Oasis Montaj MAGMAP software (MAGMAP, 2007). The power spectrum is calculated in the wave number domain so we transformed the gravity data for each window from the space domain to the wave number domain by means of a Fast Fourier Transform. This transformation produces the Gibbs phenomena in which the boundaries of the windows behave as a jump discontinuity. In order to minimize the Gibbs phenomena, we overlap the windows by 75%. Using a larger window size would result in biases from unwanted regions and would also average depths to density interfaces of structures that may be heterogeneous.
We generate the spectral curves for each window by plotting ln (power spectrum) against the wave number (k; Figure 4). Linear segments and breaks in their slopes characterize the spectral curves, which represent major density discontinuities at different depths (Emishaw et al., 2017;Fairhead & Okereke, 1987;Gómez-Ortiz et al., 2005, 2011Leseane et al., 2015;Sanchez-Rojas & Palma, 2014;Tselentis et al., 1988). Each of the linear segments is associated with a range of wave numbers and provides an indication of their average depth (Gómez-Ortiz et al., 2005, 2011. The lowest wave number segment represents a sharp density contrast between the lithosphere and the asthenosphere, while the intermediate wave number portion corresponds to the sharp density contrast between the crust and mantle, and the highest wave number segment corresponds to shallow causative bodies or noise in the data. We define the depth to the crust-mantle discontinuity (Moho) as the slope of the intermediate linear segment (Emishaw et al., 2017;Fairhead & Okereke, 1987;Gómez-Ortiz et al., 2005, 2011Leseane et al., 2015;Sanchez-Rojas & Palma, 2014;Tselentis et al., 1988) and the depth to the LAB from the slope of the steepest segment (Figure 4) as previously described by Gómez-Ortiz et al. (2005, 2011. Gómez-Ortiz et al. (2011) used the spectral analysis to show that the gravity-derived LAB depth (134 km) beneath the Iberian Peninsula is in agreement with the LAB depth (110-140 km) determined from the modeling of a combination of elevation, gravity, geoid, surface heat flow, and S wave and P wave velocity data by Fullea et al. (2010). Also, the Moho depths that were derived from spectral analysis of gravity in previous studies are in agreement with those determined from passive seismic methods (e.g., Emishaw et al., 2017;Leseane et al., 2015;Ngalamo et al., 2017). For example, Moho depth estimates in the Okavango Rift Zone (Figure 1) derived from spectral analysis of gravity data by Leseane et al. (2015) are in agreement with those obtained from passive seismic studies by Nguuri et al. (2001) with a rootmean-square error of 2.40 km.  Figure 6. The white circles with black infill represent the passive seismic stations of the Seismic Arrays for African Rift Initiation experiment (Gao et al., 2013;Wang et al., 2019). Color-coded squares represent locations of passive seismic stations of the SEGMeNT experiment (Borrego et al., 2018; Table 1). Black squares represent stations in the Ubendian Belt, red squares for stations in the Usagaran Belt, white squares for stations in the Mozambique Belt, and the purple squares for stations in the Irumide Belt. Seismically derived crustal thicknesses (Borrego et al., 2018; Table 1) for stations with underlined names where used to constrain the selection of the range of wave numbers for Moho depth estimation from the spectral curves. Dotted white lines represent the outline of Paleozoic-Mesozoic (Karoo-aged) rift structures. The white triangles represent the Rungwe Volcanic Province.

Tectonics
To minimize uncertainties associated with selecting the part of the linear segments where the slope is estimated, we calculate the crustal thickness for 1°× 1°windows centered on the locations of the 57 broadband seismometers of the SEGMeNT experiment from which the crustal thickness of the NMR and the RVP was determined from receiver function analysis [ Figure 3; Borrego et al., 2018]. We used the seismically derived crustal thickness (Borrego et al., 2018) for four of the broadband seismometer stations, each located in different tectonic regions in the NMR (MBEY, NGEA, MLOW, and KIFA; underlined stations in Figure 3 and Table 1) to constrain the selection of the wave number interval for Moho depth estimation from the power spectrum method as follows. From the power spectrum generated for the 1°× 1°windows centered on the location of the four aforementioned seismic stations, we selected the range of wave numbers in the intermediate linear segment of the spectral curve for which the slope of the polynomial fit is approximately equal to the seismically derived crustal thickness (Borrego et al., 2018). The rest of the stations are used for comparison of the results of the two approaches ( Table 1).
The selection of the wave number interval for LAB estimation is clearer because each of the spectral curves show a discontinuity between the lower wave number segment and the higher wave number segments (Figure 4). We determine the error associated with selecting the part of the linear segments of the spectral curves for the Moho and LAB depth determination by using three first-order polynomial fits to produce straight lines. The average slope of the lines represents the depth to a density discontinuity. Subsequently, for the corresponding linear segment, we determine the standard deviation of each fit. We use the combination of the variations of the slope measurements from three different polynomial fits and their standard deviations to determine that the error in estimating the depth to Moho and LAB is~±2 and~±4 km, respectively. Apart from the estimated experimental fitting errors, another source of error might arise from the data quality given that the terrestrial gravity data are relatively sparse in parts of the MR. Since the Bouguer gravity anomalies from the World Gravity Model 2012 are more sensitive to long-wavelength anomalies, in order to remove short-wavelength artifacts from our lithospheric structure that is estimated at 0.25°× 0.25°intervals, we reinterpolate the crustal and lithospheric thicknesses to 0.75°× 0.75°using the minimum curvature technique in the Generic Mapping Tool (Wessel et al., 2013). Reinterpolating the crustal and lithospheric thicknesses to 0.75°× 0.75°also allow for tests of boundary effects on the numerical model (Supporting Information Figure S2).

2-D Forward Gravity Modeling
We use 2-D forward modeling of the gravity data along profiles A-A′, B-B′, and C-C′ (see Figure 3 for locations of the profiles) using the 2-D GYMSYS Oasis Montaj software (Talwani et al., 1959) to determine the lithospheric structure beneath the MR. Profiles A-A′ and C-C′ follow the passive seismic stations of the Seismic Arrays for African Rift Initiation experiment (Gao et al., 2013;Wang et al., 2019), and B-B′ is an east-west profile that cuts across the widest part of the MR and includes portions of the Luangwa Rift in Zambia. Gravity models have the limitation of nonunique solutions; however, we construct geological models by using realistic constraints on rock densities and lateral variations of rock units (e.g., Fletcher et al., 2018;Leseane et al., 2015;Mahatsente et al., 1999;Mickus et al., 2007;Simiyu & Keller, 2001). Because these constraints are not readily available in all parts of the MR and surrounding areas for the 2-D forward modeling of the gravity data, we used as initial constraints, estimates of crustal thicknesses obtained from the inversion of Rayleigh wave phase velocity maps constructed from ambient noise tomography recorded by broadband (Seismic Arrays for African Rift Initiation) seismic data sets, which range between 36-38 km for the A-A′ profile and 38-40 km for the B-B′ and C-C′ profiles (Wang et al., 2019). Similarly, we use seismically derived depth to LAB ranging between 150-170 km for profile A-A′ and 170-210 km for profiles B-B′ and C-C′ (Fishwick, 2010) as initial constraints of the 2-D forward modeling of the gravity data. The thickness and densities of the upper crust (2.76 g/cm 3 ), the lower crust (2.92 g/cm 3 ), and the lithospheric mantle   Tectonics determine the final gravity model by varying the initial crustal thickness and depth to LAB by~10% and the initial densities of the different Earth layers by~15%.

3-D Numerical Modeling of EDC
We simulate EDC in the asthenosphere due to variation in the lithospheric structure using the adiabatic boundary approach defined in the open-source finite element code Advanced Solver for Problems in Earth's ConvecTion (Bangerth et al., 2018a;Bangerth et al., 2018b;Heister et al., 2017;Kronbichler et al., 2012) to test the potential role of EDC in magma-poor rifting along the MR. We generate an EDC model by solving for the velocity term of the Stokes' law, which is the momentum and continuity equation. Our model domain has dimensions of 1,400 × 1,600 × 500 km (Supporting Information Figure S1a). We use our derived lithospheric structure to constrain the initial temperature condition by implementing the LAB in the model as an isothermal boundary (1,673.15 K). We approximate a linear gradient from the surface (273.15 K) to the LAB, below which the temperature increases adiabatically (Supporting Information Figure S1a). The temperature variations due to variations in the lithospheric thickness will lead to lateral density perturbations in the asthenosphere that drive small-scale convection (King & Anderson, 1998). We use non-Newtonian visco-plastic rheology in the crust with wet quartzite material parameters [Supporting Information Table S1; Ranalli, 1995], and the lithospheric mantle and sublithospheric mantle are respectively dominated by dislocation-creep and diffusion-creep flow laws of dry olivine (Karato & Wu, 1993;Kohlstedt et al., 1995;Koptev et al., 2015). We also test the influence of wet olivine as the dominant composition in the upper mantle [Hirth & Kohlstedt, 2004; Supporting Information Table S1] but choose dry olivine following Karato and Wu (1993), Kohlstedt et al. (1995), and Koptev et al. (2015) and because the region is far from fertile oceanic lithosphere and subduction zones where fluids may be entrained. Further, we test both free-slip (tangential flow) and fixed (zero velocity) boundary conditions at all sides of the model to assess edge effects, which are minimal on the region of study (see supporting information for details).

Moho Depth Estimates
The Moho depth estimates beneath the MR and surroundings derived from the spectral analysis range between 35 and 45 km (Figure 5a). There is relatively thin (~38-40 km) crust beneath most of the NMR compared to a consistent pattern of thicker (~40-45 km) crust beneath the surface expression of the CMR and SMR. This thick crust is unexpected since rifted regions typically have stretched crust. Outside of the rift there is thick crust (40-45 km) beneath the region slightly north of the RVP and an E-W belt of thick (~42-45 km) crust beneath the Niassa Craton. The thickest crust (~45 km) is in the southernmost region of the SMR beneath the Mulanje Mountain where the rift changes orientation from a dominantly N-S trend to a NNE-SSW trend. Generally, the crust is thinner (~35-39 km) beneath the surrounding Permo-Triassic Karoo rifts, in particular the northern Luangwa Rift valley, the Zambezi Rift, and the Ruhuhu Trough.
Beneath the N-S extent of the MR, the 2-D forward gravity models ( Figure 6) indicate a relatively flat Moho with depth ranging between 38 and 45 km. Along profile A-A′ (Figure 3), the crust is~42 km beneath the northern half of the RVP. South of the RVP, the crust is thinnest (~39 km) beneath the axis of the northern part of Lake Malawi. The crustal thickness increases progressively southward to~45 km beneath the SMR  (Borrego et al., 2018) for the underlined stations where used to constrain the selection of the range of wave numbers for Moho depth estimation from the spectral curves.

10.1029/2019TC005549
Tectonics ( Figure 6a). For the profiles across the MR, profile B-B′ indicates a uniform crustal thickness of~40 km beneath the rift axis and the rift flanks (Figure 6b). For profile C-C′, the crustal thickness beneath the rift axis and the western rift flank is~43 km and decreases to~40 km beneath the eastern rift flank (Figure 6c).

LAB Depth Estimates
Estimates of the depth to the LAB beneath the MR and surroundings based on the spectral analysis range from 115 to 210 km (Figure 5b). We find an elevated LAB underlies nearly the entire length of the MR. The depth to the LAB has an average thickness of~140 km beneath the surface expression of the MR with the thinnest (115-130 km) lithosphere underlying the NMR and SMR. The CMR is underlain by an E-W trending belt of relatively thick lithosphere that ranges from 180 to 210 km thick beneath the Precambrian terranes of the eastern and western flanks of the rift.
Consistent with the spectral analysis-derived LAB, the 2-D forward gravity modeling approach indicates thinning of the lithosphere beneath the surface expression of the MR (A-A′ profile; Figure 6a). The depth to LAB varies between 110 and 160 km along the profile. The lithosphere is thinner beneath the NMR (~110 km) and SMR (~120 km) compared to the CMR (~140 km). The models along profiles B-B′ and C-C′ that are transverse to the MR indicate symmetrical thinning of the lithosphere beneath the rift axis tõ 120 km (Figure 6b) and~143 km, respectively (Figure 6c), as opposed to~180-210-km lithospheric thickness beneath the rift flanks.  Figures 7a and 7b show lithospheric deformation and asthenospheric flow patterns resulting from our numerical modeling of EDC. Our results indicate that there may be asthenospheric upwelling (~3 cm/year) beneath the RVP and beneath the NMR and SMR driven by EDC. Minor upward (~1-cm/year vertical velocities) and diverging (<0.5-cm/year horizontal velocities) displacement of the lithospheric mantle at~150 km in the NMR and to the southeastern flank of the MR (Lurio Shear Zone) is present. Deformation patterns in the lithospheric mantle are spatially consistent with modeled EDC in the asthenosphere at~250-km depth (Figure 7b), indicating the coupling of the asthenosphere to the lithosphere. Crustal velocities are negligible, which may indicate crust-mantle decoupling. Our model is used to suggest a southwestward flow of the upwelling mantle in the NMR towards the cratonic lithosphere of the CMR. At deeper depths (~350 km), there is change in the direction of horizontal mantle flow with an inward radial flow pattern beneath the zones of mantle upwelling, indicating that the depth limit of the upwelling is shallower than 350 km (Supporting Information Figures S3 and S4).

Crustal Structure
The crust beneath the NMR is relatively thin compared to a consistent pattern of thicker crust beneath the CMR and SMR. The thickness of the crust beneath the MR thus increases progressively southward. This observed pattern of southward increase in the crustal thickness is in accord with models advancing southward propagation of the rift and its opening in a zipper-like fashion from north to south (Betzler & Ring, 1995;Flannery & Rosendahl, 1990).
Although the crust is relatively thin in the NMR, it is generally thicker than seismically derived crustal thicknesses along the MR (Borrego et al., 2018;Wang et al., 2019). The differences in our gravity-derived crustal thicknesses and seismically derived crustal thicknesses may result from the uncertainties in both techniques and/or from the quality of the data. Alternatively, the differences could result from the distinct definitions of the Moho from each technique. The spectral analysis determines the Moho from to the density contrast adjacent to the LAB, but the seismically derived Moho is defined by the first instance (top down) of the S wave velocity that reaches~4.3 km/s (Borrego et al., 2018;Wang et al., 2019). Hence, it is possible that there is mafic material constituting a magmatic underplating, which is accreted beneath the lower crust, and may also have an S wave velocity of~4.3 km/s. This leads the seismic technique to image the Moho as the top of the magmatic underplating. However, there might be density similarity between the lower crust and the magmatic underplating and greater density contrast between the base of the magmatic underplating and underlying lithospheric mantle. Hence, this leads the gravity technique to image the Moho as the base of the magmatic underplating.
The presence of magmatic underplating in the SMR is supported by the Mw 6.1 March 10, 1989 earthquake with a hypocenter of~30 km caused by a dip-slip rupture of the Bilila-Mtakataka Fault [ Figure 2a; Jackson & Blenkinsop, 1993, 1997. Because of the deep focus of the earthquake, we suggest the presence of a thick crust with a strong and brittle lower crust beneath the SMR. It has been argued, on the basis of the absence of Cenozoic magmatism, that the magmatic underplating in the SMR took place during the Mesozoic and possibly contributed to the magmatic evolution of the Cretaceous Chilwa Alkaline Province (ring complexes) including the Mulanje Mountains [ Figure 5b; Platt & Woolley, 1986].
Magmatic underplating may also be present in the NMR beneath the RVP. In the NMR, we compare our Moho depth estimates obtained from the spectral analysis with Moho depths obtained from 4 of the 57 seismic stations (CRTR, IGOM, MZUN, and NKAL) based on the Borrego et al. (2018) study ( Figure 3; Table 1). The four stations were chosen because each is positioned on volcanic extrusive of the RVP according to Borrego et al. (2018). Beneath the RVP, Moho depth estimates obtained from the spectral analysis are deeper than those obtained by Borrego et al. (2018) with an average difference of 3.1 km. Similar to the SMR, we propose that the presence of the~43-km thick crust beneath the RVP and the discrepancy in the Moho depth estimates between the spectral analysis and the passive seismic imaging is due to magmatic underplating beneath the crust of the NMR. However, different from the SMR, it appears that the magmatic underplating beneath the RVP took place during the Cenozoic (17 Ma; Mesko et al., 2014). This might be due to

10.1029/2019TC005549
Tectonics decompression melting of the asthenosphere as the result of thinning of the lithosphere providing a source of the magmatism in the RVP that predates the 8.6-Ma rifting in the NMR.
We also compare our Moho depth estimates with Moho depth estimates from the remaining seismic stations of Borrego et al. (2018; Table 1). These seismic stations are located within different Precambrian tectonic blocks and are shown with different colored seismic stations in Figure 3. We find the average difference in the Moho depth estimates from the seismic stations within the Ubendian orogenic belt and the Moho depth estimates from the spectral analysis to be -0.4 km. Additionally, we find the average difference within the Usagaran orogenic belt to be 1.4 km and that within the Irumide orogenic belt to be −2.7 km. We suggest that these differences are within the error ranges of the two techniques.

Lithospheric Mantle Structure
Unlike the crust that does not show a pattern of thinning typical of continental rifts, we observe thinning of the lithospheric mantle beneath the entire MR. The thinning suggests symmetrical stretching of the lithospheric mantle beneath the MR and has a stretching factor of 1.54. We calculate the lithospheric stretching factor by dividing the thickness of the lithosphere outside the MR (200 km) by the thickness of the lithosphere beneath the rift (130 km; Figure 5b). Thinning of the lithosphere may generate upwelling of the asthenosphere due to EDC, whereas the relatively thick zones will cause downwelling (represented by warm and cool colors, respectively, in Figure 7b). Our results reveal that the southern part of the CMR is underlain by an E-W trending belt of 180-210-km thick lithosphere. The presence of this thick lithosphere beneath the southern part of the CMR had previously been imaged by Fishwick (2010) using surface wave tomography. We propose that this thick lithosphere represents the NE continuation of the Niassa Craton-a little cratonic block in eastern Zambia and that has been identified using magnetotelluric imaging (Sarafian et al., 2018). In addition, Reed et al. (2016) found an apparent 20-km shallowing of the mantle transition zone that they attribute to high-velocity anomalies in the upper mantle. This result from Reed et al. (2016) also supports the

10.1029/2019TC005549
Tectonics existence of relatively thick and strong lithosphere in the south-central MR. It is noteworthy mentioning that the CMR is the widest segment of the rift and is characterized by border faults that have the lowest vertical throw as inferred from the height of their escarpments (Laó-Dávila et al., 2015). This can be attributed to lower rate of extension and extension diffusion due to the presence of the thick, cold, and strong cratonic lithosphere. Diffusion of extension is also observed in the Mbulu faulted domain due to the propagation of the Northern Tanzania Divergence into the Tanzanian Craton (Le Gall et al., 2008).
The zone of thin lithospheric mantle we observed beneath the NMR (Figure 5b) spatially coincides with an elevated topography within the Nyika Plateau, the Livingstone Mountains, and the RVP (Figure 2a). O'Donnell et al. (2015) and Accardo et al. (2017) found a low-velocity anomaly beneath the RVP that can be interpreted as due to the presence of anomalous thermally modified lithosphere or ascending asthenosphere. Here, we propose that the thinned lithosphere observed beneath the NMR is due to stretching and possibly thermal erosion of the base of the lithospheric mantle that triggered asthenospheric upwelling through EDC, hence, resulting in the elevated topography (Figure 2a) through rapid tectonic uplift.
We also observe that the thinning of the lithospheric mantle and the resulting asthenospheric upwelling due to EDC align with the trend of the Precambrian shear zones (Figures 2a, 2b, and 7b). For example, the broad zone of thin lithospheric mantle beneath the NMR coincides with the Mughese Shear Zone, which extends within the Ubendian orogenic belt. Similarly, the zone of thin lithospheric mantle beneath the SMR aligns with the Lurio Shear Zone. Viola et al. (2008) suggest that the Lurio Belt (which includes the Lurio Shear Zone) was a zone of intense lithospheric contraction that accompanied collision of the Marrupa-Nampula and Unango complexes during the Neoproterozoic Mozambique orogeny (Figure 2b). Sarafian et al. (2018) proposed that the observed thin lithospheric mantle beneath the Neoproterozoic Mwembeshi Shear Zone (Figures 5b) was a site of zonal lithospheric mantle delamination at the southern edge of the Bangweulu Cratonic Block, and this might have resulted in localization of extension during the onset of rifting in the Luangwa Rift. Further, Sarafian et al. (2018) interpreted the conductive signature of this shear zone to be indicative of upward movement of asthenospheric fluids to shallower lithospheric depths. Hence, similar processes can be used to explain the thinner lithospheric mantle beneath the Lurio Shear Zone.

Decoupling of the Crust From the Lithospheric Mantle
The 2-D forward gravity models (Figures 6a-6c) show that the MR has a crustal stretching factor of~1.1. This is calculated by dividing the thickness of the crust outside the rift (~43 km) by the crustal thickness beneath the rift (~40 km). In addition, the 2-D forward gravity models (Figures 6a-6c) indicate that the rift has a lithospheric stretching factor of~1.4 to 1.7. This is calculated by dividing the thickness of the lithosphere outside the rift by the thickness of the lithosphere beneath the rift, that is, 200/143 km ( Figure 6c) and 200/120 km ( Figure 6b). We observe from the 2-D forward gravity modeling thinning of both the crust and the lithospheric mantle ( Figure 6a); however, the weak EDC is not accompanied by observable surface deformation in the region. Thinning of the lithospheric mantle beneath the entire MR, the lack of crustal thinning beneath the CMR and SMR, and the insignificant crustal deformation due to the EDC indicate decoupling of the crust from the lithospheric mantle during extension. This is similar to results of 3-D gravity modeling in the Iberian Peninsula by Torne et al. (2000), which found a relatively flat Moho underlain by a significantly thinned lithospheric mantle indicative of decoupling of extension at the base of the crust.
Decoupling of extension has been previously observed in narrow rifts such as the Valencia Trough (NW Mediterranean; Gaspar-Escribano et al., 2003). It has also been observed in wide rifts such as the basin and range (Buck, 1991). Several models highlight the ductile nature of the lower crust as the main mechanism of decoupling of extension (Buck, 1991;Gaspar-Escribano et al., 2003;Handy et al., 2005;Ter Voorde et al., 1998), which is accounted for in our numerical modeling. Another mechanism assumes that subhorizontal faults, or shear zones deep in the crust, result in decoupling of extension in the upper crust from that in the lithospheric mantle (Reston, 1990).

Implications for Magma-Poor Rifting
We suggest that in areas of magma-poor rifting such as the western branch of the EARS, zones of localize crustal magmatism are characterized by coupled crust-lithospheric mantle extension, and rift segments far field of the magma centers are characterized by decoupled crust-lithospheric mantle extension. This

10.1029/2019TC005549
Tectonics suggestion provides insight into the lithospheric structure of rift segments far field of isolated magmatic centers along the western branch of the EARS (Figure 1).
We propose that the NMR underwent magma-assisted rifting that resulted in coupling of the extension between the crust and the lithospheric mantle. The presence of thin lithosphere (Figures 5b and 6a) allowed asthenospheric upwelling beneath the RVP and beneath the NMR and SMR (Figure 7b), and this resulted in the production of asthenospheric decompression melt. In the NMR, melt may intrude the crust through shear zones such as the Mughese Shear Zone (Figures 2b and 8). The presence of decompression melt in the NMR is supported by the presence of uppermost mantle LVZ at depth of 68 km (O'Donnell et al., 2013). Southwestward flow of the upwelling asthenosphere from the NMR (Figures 7b and 8) towards the CMR possibly led to thermal erosion of the base of the thick cratonic lithosphere, thereby enabling localization of extension in this part of the MR. Furthermore, rifting of the thick cratonic lithosphere beneath the south-central MR might have been additionally facilitated by asthenospheric fluids transported by the southwestward flow of the upwelling asthenosphere beneath the NMR. Typically, percolation of asthenospheric fluids leads to metasomatic refertilization of the lithospheric mantle (Schutt & Lesher, 2010). Adding fluids to peridotite of the initially dry and depleted lithospheric mantle can mechanically weaken it, hence, allowing for localization of extension (Peslier et al., 2012). Similarly, lithospheric extension in the SMR might have been partially facilitated by mechanical weakening of the lithospheric mantle through metasomatic refertilization from fluids sourced from deep asthenospheric upwelling (Figure 7b).

Conclusions
In this study, we characterize the lithospheric structure beneath the MR and surroundings using the spectral analysis and 2-D forward modeling of the WGM 2012 Bouguer gravity anomalies. We find thinner crust beneath the NMR and a relatively thick crust beneath the CMR and SMR. We also find thin lithosphere beneath the entire length of the MR with the exception of an E-W trending belt of relatively thick lithosphere across the southern part of the CMR. Further, we use our new lithospheric structure model to constrain 3-D numerical modeling of EDC beneath the MR, which indicates asthenospheric upwelling beneath the thin lithosphere of the NMR and SMR with possible generation of decompression melt. We infer, based on the

10.1029/2019TC005549
Tectonics presence of the RVP together with both thin crust and lithosphere, that the NMR experienced and is experiencing coupling of extension between the crust and lithospheric mantle. In contrast, we infer from the apparent absence of surface magmatism, together with the presence of thick crust but thin lithosphere, that the CMR and SMR are experiencing decoupling of extension between the crust and the lithospheric mantle where extension is mostly accommodated within the lithospheric mantle. We, therefore, conclude that magma-assisted rifting is an important mechanism for the development of magma-poor continental rifts, particularly in isolated magmatic zones where there is coupled crust-lithospheric mantle extension and lithospheric weakening by localized injection of asthenospheric melt. Whereas, fluids-assisted rifting, where fluids are migrating from far distance asthenospheric upwelling sources, is an important mechanism for the development of segments of the continental rifts where magmatism is apparently absent because extension in the lithospheric mantle, which is facilitated by metasomatic refertilization from deep asthenospheric fluids may be decoupled from the crust.