Geodynamic models suggest the interaction of buoyancy-driven mantle convection and plate tectonics. In fact, large-scale topography features on Earth depend on density variations in the upper mantle (Steinberger, 2016) and buoyancy of the lithosphere due to thermochemical heterogeneities (Kelly et al., 2003; Lee et al., 2005; Shulgin and Artemieva, 2019) contribute to the long stability of cratons over geological times (Lee et al., 2005; Eaton and Claire Perry, 2013). Despite the importance of subsurface density structure in mantle convection and lithospheric deformation, the density structure in the subsurface is poorly constrained. Gravity and seismic techniques are the two main methods employed to investigate the density structure in the subsurface (e.g., Darbyshire et al., 2000; Vacher and Souriau, 2001; Tondi et al., 2012; Deng et al., 2014; Kaban et al., 2016). However, there is intrinsic non-uniqueness in gravity modeling, thus the solutions need to be constrained, which introduces uncertainties in the final model. For example, due to the expected high heterogeneity in the crust, it is assumed that short wavelength features in the gravity field generally have a crustal origin while long wavelength anomalies have a deeper source. On the other hand, tomography models provide increasingly detailed images of the seismic velocity structure of the upper mantle.

In order to quantitatively investigate the implications of existing tomography models for the dynamics of the upper mantle, most geodynamic models use simple conversions of seismic velocity perturbations to density anomalies through the design of a conversion factor (e.g. Karato and Wu, 1993; Becker, 2006; Conrad and Lithgow-Bertelloni, 2006; Steinberger and Calderwood, 2006; Conrad et al., 2007; Karato, 2008; Conrad and Behn, 2010; Ghosh et al., 2010; 2017; Wang et al., 2015; Liu and Zhong, 2016; Adam et al., 2021). Previous studies assume a positive constant value for this factor for the entire mantle (e.g., Becker, 2006; Conrad and Lithgow-Bertelloni, 2006; Conrad et al., 2007; Conrad and Behn, 2010; Ghosh et al., 2010; 2017; Wang et al., 2015; Liu and Zhong, 2016). For example, Becker (2006), Conrad and Lithgow-Bertelloni (2006), and Conrad and Behn (2010) use a conversion factor between shear velocity anomalies and density anomalies of 0.15, while Ghosh et al. (2010, 2017) and Wang et al. (2015) assume a conversion factor of 0.25. Some studies consider the mineral physics approach of Stixrude and Lithgow-Bertelloni (2005a, 2005b, 2007, 2011) to derive depth-dependent values for the conversion factor typically varying between 0.2 and 0.4 (e.g., Steinberger and Calderwood, 2006; Adam et al., 2021). Other studies determine the conversion factor through laboratory experiments to range between 0.2 and 0.4 (Karato and Wu, 1993; Steinberger and Calderwood, 2006; Karato, 2008). An important observation is that these cited studies all derive positive values for the conversion factor, which implies that every positive velocity perturbation will convert to positive density anomalies and vice versa, which is not always true in nature. For example, the ultra-low velocity zones (ULVZs) at the core-mantle boundary have extremely low seismic velocities but have higher density than their surroundings in order to remain stable on the core-mantle boundary. Another example is given by cratonic keels, which are mostly characterized by negative density anomalies due to compositional effects (e.g., Kaban et al., 2003; 2016; this study) and by very high seismic velocities due to low temperatures.

In this study, we investigate the density structure beneath the Rungwe Volcanic Province (RVP; Figure 1A, B) and surroundings, which is the southernmost volcanic center in the magma-poor Western Branch of the East African Rift (Figure 1A). The RVP lies within the Ubendian-Usagaran mobile belts that circumvent the thick lithosphere of the Tanzanian and Bangweulu cratons (Figure 1B; e.g.; Corti et al., 2007; Fritz et al., 2013) and consists of three large volcanoes, Ngozi, Rungwe, and Kyejo, that were last active about 500 years ago (Fontijn et al., 2010; black triangles; Figure 1B). Magmatism in the RVP is highly localized with past eruptions covering ∼1,500 km² (Figure 1B; Ebinger et al., 1989; 1997; Fontijn et al., 2012). ⁴⁰Ar/³⁹Ar radiometric dating of samples from the RVP suggest that magmatism in the RVP started by 19 Ma (Mesko et al., 2014; Mesko, 2020) and possibly as early as ∼25 Ma (Roberts et al., 2012) which predates rifting in the northern Malawi Rift at ∼8.6 Ma (Ebinger et al., 1993). These ages suggest that magmatism in the RVP might have played an important role in thermally weakening the lithosphere, thereby facilitating rifting.

A number of previous studies have been conducted to investigate the crustal and upper mantle structures beneath the RVP and surrounding using controlled-source seismic (lake-bottom seismometers), body and surface wave tomography, receiver function stacking and from the analysis of gravity and aeromagnetic data (O’Donnell et al., 2013; Accardo et al., 2017; Grijalva et al., 2018; Njinju et al., 2019a; Njinju et al.; 2019b; Emry et al., 2019). O’Donnell et al. (2013) used Rayleigh wave phase velocities to invert for a 3-D shear wave velocity model and observed pronounced velocity lows beneath the RVP. Grijalva et al. (2018) used P and S wave seismic tomography to also image the low velocity zones (LVZ) beneath the RVP and the northern Malawi Rift and interpret it to represent the flow of warm, mantle superplume upwelling from the southwest, beneath, and around the thick lithosphere of the Bangweulu Cratonic Block. Although the source of these LVZs still remains enigmatic, their presence suggests density perturbations beneath the RVP and surroundings. Lateral and vertical variability of in situ density structure beneath the RVP may be inferred from a joint inversion of seismic and gravity data (Forte and Claire Perry, 2000; Deschamps et al., 2002; Tondi et al., 2012), however, the main limitation is the non-uniqueness of the relationship between seismic velocities and density (Barton, 1986; Deschamps et al., 2001; Lee, 2003; Brocher, 2005).

We use shear-wave velocity perturbations ( d l n v s ) from Emry et al. (2019) as a reference model to perform a constrained inversion of satellite gravity data centered on the RVP. We use the code jif3D (Moorkamp et al., 2011) with a d l n v s -density coupling criterion based on mutual information (MI; Haber and Holtzman Gazit, 2013; Mandolesi and Jones, 2014; Lösing et al., 2022; Moorkamp, 2022), which constructs a one-to-one relationship between the d l n v s and density directly from the gravity data without considering an empirical relationship. In the inversion, a joint probability distribution of the gravity data and the reference model ( d l n v s ) is estimated and its entropy minimized in order to generate a density model that is statistically compatible with the reference model. The algorithm generates a 3D density model of the lithosphere and sublithospheric mantle beneath the RVP and surroundings down to a depth of 660 km by seeking a combined d l n v s -density model that fits the gravity data with maximum correspondence between the d l n v s and density. Our results reveal a conspicuous negative density anomaly (−200 kg/m³) in the sublithospheric mantle (extending at depths from ∼100 km to ∼250 km) beneath the central part of the Malawi Rift extending to the west, beneath the Niassa Craton. This negative density anomaly beneath the Niassa Craton is coincident with locations of positive seismic velocity perturbations (+7%). We calculate a 3D model of the velocity-to-density conversion factor (f), which is defined as the ratio of the derived density perturbations to the shear-wave velocity perturbations (e.g., Root et al., 2017; Liang et al., 2019). We find negative f-values beneath the Niassa Craton which suggest the observed negative density anomaly is mostly due to compositional variations. Apart from the Niassa Craton, there are generally positive f-values in the study area, which suggest dominance of temperature control on the density structure. Although the RVP generally shows negative density anomalies and positive f-values, at shallow depths (<120 km), f ≈ 0, which suggests important contributions of both temperature and composition on the density structure beneath the RVP. We suggest that the presence of plume material at shallow depths beneath the RVP thermally reduces the upper mantle density. The plume material contributes metasomatic fluids that precipitate dense minerals that slightly increase the density due to compositional variations. We suggest the negative buoyancy of the Niassa Craton contributes to its long stability, while constituting a barrier to the southward flow of plume material, thus restricting the southward continuation of magmatism in the Western Branch of the EAR. The presence of a negative-density anomaly where d l n v s values are positive is incompatible with models based on the use of simple d l n v s to density conversion factors. Thus, these results have implications on how d l n v s models are converted to density perturbations.

Gravity exploration has the advantages of low cost and high efficiency but its low vertical resolution hampers the inversion process. Therefore, some a priori information can be used as constraints to enhance the inversion accuracy. Here, we used seismic velocity perturbations ( d l n v s ; Emry et al., 2019) as reference models to investigate the density structure beneath the RVP and surroundings through a constrained inversion of satellite gravity data (Pail et al., 2018).

We use the code jif3D (Moorkamp et al., 2011) with a d l n v s -density coupling criterion based on mutual information (MI; Haber and Holtzman Gazit, 2013; Mandolesi and Jones, 2014; Lösing et al., 2022; Moorkamp, 2022). MI is an unsupervised machine-learning approach that quantifies the amount of information shared by two random variables. For our model, the MI coupling constructs a one-to-one relationship between the reference model ( d l n v s ) and density directly from the gravity data without requiring a pre-defined relationship. During the inversion, the joint probability distribution of the gravity data and the reference model is estimated and its entropy minimized in order to guide the inversion results towards a solution that is statistically compatible with the d l n v s -model. Although the MI technique does not presume any explicit relationship between the estimated density-model and the reference d l n v s -model, if a link exists in the gravity data, then it is highlighted in the estimation of the joint probability distribution (Eq. 1a).

We use a kernel density approach with a Gaussian kernel to calculate the MI between d l n v s and density in three main steps:

(1) We estimate the probability density distribution (pdd) for the joint parameters and the marginal distributions for each parameter. If we denote the pairs of transformed parameters in each model cell as ( x i ; y i ) where i = 1 … M and M is the number of cells used to discretize the inverse problems, the pdd is approximated as: p ρ j , d l n v s , k = 1 2 π   M ∑ i = 1 M exp x i − ρ j 2 + y i − d l n v s , k 2 2 σ 2 (1a) where ρ j and d l n v s , k are the discrete variables at which we approximate the pdd. Since we extract the Bouguer gravity data with a resolution of 0.1° (∼11 km), only subsurface density anomalies >11 km (in each dimension) will be properly resolved by the gravity data. We therefore refine the entire model domain (with dimensions of ∼1,000 x ∼890 x 660 km; Figure 3B) to a global mesh refinement of 6 such that each cell is ∼16 x 14 x 10 km. For this study, we use M = 64 (i.e., 2⁶, where 6 = global mesh refinement) evenly spaced values across the expected parameter range for each parameter and choose the standard deviation of the Gaussian, σ, as half the discretization width.

From the joint pdd we can also calculate the marginal pdd’s: p ρ j = ∑ k p ρ j , d l n v s , k (1b) and p d l n v s , k = ∑ j p ρ j , d l n v s , k (1c)

(2) We then calculate the Shannon Entropy of the marginal and joint pdd’s:

The full objective function that we minimize iteratively includes the data misfit term for the gravity anomaly data (Φ _d,grav ), the regularization term (Φ _reg,ρ ), and the mutual-information term (MI): Φ ρ = λ g Φ d , g r a v ρ + λ ρ Φ r e g , ρ ρ + λ M I M I ρ , d l n v s (4)

Details of each term in the objective function is given in Moorkamp et al. (2011). The gravity data misfit (Φ _d,grav ) is defined as the l ² norms of the data residuals weighted by the uncertainties (10 mGal). The data misfit for the density model ρ with respect to the gravity data set d is given by: Φ d , g r a v ρ = g ρ − d o b s , g r a v T C d − 1   g ρ − d o b s , g r a v (5) where g ρ is the synthetic gravity data from the forward calculation for the given ρ -model, and d o b s , g r a v , is the observed gravity data. C d − 1 is the inverse of the data covariance matrix, which is included in the objective function to reduce the influence of uncertainties or errors in the gravity observations.

The regularization term (Φ _reg,ρ ) is included in the objective function in order to stabilize the inversion and yield smooth models and possibly coherent and geologically meaningful structures (e.g., Jupp and Vozoff 1975; Zhdanov 2002). Φ _reg,ρ is the ρ -model norm, which in our case is defined as the second spatial derivative of the ρ -model weighted by individual a priori model covariance matrix. In more details, the regularization term is given by: Φ r e g , ρ ρ = ∑ i α i ρ − ρ o T W i T C M − 1 W i   ρ − ρ o (6)

The model regularization term is obtained by summing over the three axis directions i = {x, y, z} and weight the contribution for each direction by a weight α i . In order to obtain a smooth ρ -model with minimum curvature between adjacent cells, we choose the matrices W i as finite difference approximations to the second spatial derivative of the ρ -model. The inversion is regulated by keeping the result ( ρ -model) close to an a priori reference model ρ o ; C M − 1 is the inverse of a diagonal model covariance matrix whose elements can be changed to limit the variation of certain parts of the model. We use an a priori model covariance that increases as the square of depth to balance the decrease in sensitivity of the gravity measurements with the distance from the causative body (Li and Oldenburg, 1998). The parameters λ _g , λ _ρ , and λ _MI are weights that control the balance between the different terms of the objective function. The choice of the regularization parameter is guided by the analysis of trade-off curves or L-curves (Hansen, 1992; Supplementary Figure S1).

The objective function is minimized iteratively using the limited memory quasi-Newton approach (L-BFGS, Avdeev and Avdeeva, 2009) for parallel forward solvers; flexible model parametrization and large-scale model optimization. Each model cell (rectangular prism with dimensions ∼16 x 14 x 10 km) is assigned a value from the starting density model ρ o and the starting model is optimized iteratively until the misfit reaches a minimum or the user prescribed number of iterations is exhausted. The gravity forward problem is based on a parallelized analytical solution of the gravitational attraction of adjacent cells (Moorkamp et al., 2010). We first ran independent gravity inversions for up to 100 iterations in order to determine the gravity data weight λ _g and the density regularization weight λ _ρ . We found an acceptable trade-off that minimizes both the data-norm and the model-norm when λ _g = 1 and λ _ρ = 10^3.5 (Supplementary Figure S1). For the constrained inversion, we start with a large weight for the MI term (λ _MI = 10⁸) and successively relax the weight until the inversion has converged to a data misfit comparable with the misfit for inversions based on the gravity data set individually, which occurred for λ _MI = 10⁶ (Supplementary Figure S2). The inversion termination criteria were the same as the independent gravity inversion termination criteria except that we ran the constrained inversion for up to 300 iterations since the minimization of the full objective function converges more slowly than the independent gravity inversion and requires a larger number of iterations to reach a comparable level of fit.

The shear-wave velocity anomalies beneath the RVP and surroundings from Emry et al. (2019) that we use as our reference model exist only for the upper mantle (>33 km depth down to 660 km). Thus, the density structure in the crust (<33 km) is only from the gravity inversion without a seismic constraint. However, an extraction of the crustal density structure at 20 km and 30 km depths (Supplementary Figures S6A, B) reveals geologically meaningful structures like negative density anomalies (<−50 kg/m³) beneath rift basins, possibly due to the low density of rift basin sedimentary covers and a positive density anomaly (∼+200 kg/m³) beneath the Mughese Shear Zone. We note again that only subsurface density anomalies >11 km can be properly resolved by the gravity data, thus more heterogeneous smaller wavelength (<11 km) crustal density structures will be unresolved by our technique. We focus our interpretation of the vertical and lateral variability of in situ density structures beneath the RVP and surroundings at depths >33 km. Figures 4A–D respectively show 100, 150, 200, and 250 km depth slices of density anomalies resulting from our constrained inversion of gravity data with d l n v s -model as the reference.

Our results reveal a conspicuous negative density anomaly (−200 kg/m³) in the sublithospheric mantle (with depth extending from ∼100 to ∼250 km) beneath the central part of the Malawi Rift extending to the west, beneath the Niassa Craton. This negative density anomaly beneath the Niassa Craton is coincident with locations with positive seismic velocity perturbations (+7%). At depths >100 km, there is a region of negative density anomalies (−50 to −10 kg/m³) that extend from the Eastern Branch, through the Usangu Basin (Figure 1B), to focus beneath the RVP, which is coincident with locations with low seismic velocity anomalies (−7%). This result suggests a possible linkage of the Western Branch and the Eastern Branch through the Usangu Basin (Figure 1B; Mulibo, 2022).

Within the low-density region beneath the RVP, we highlight a denser anomaly that is highly focused beneath the volcanic centers (denoted by red triangles in Figure 4B). In order to better visualize the density anomaly at depth, we extract a cross-section along a N-S profile that cut across the RVP (Figure 5) and observe that the denser anomaly beneath the RVP is confined at shallow depths within the asthenosphere (∼110–150 km) where the lithosphere is thinnest (∼110 km; Fishwick, 2010; Njinju et al., 2019a). Figure 5 also indicates the fit between gravity observations and the forward calculated gravity and shows that the Niassa Craton extends to depths reaching ∼350 km, however, the buoyant cratonic nucleus is restricted at a depth range of 100 to 250 km. We suggest that the denser anomaly in the lithospheric mantle beneath the RVP is possibly due to compositional heterogeneity from metasomatization, while the negative buoyancy of the Niassa Craton contributes to its long stability. See Section 5 in which we discuss a possible origin of these density anomalies.

We present a 3D model of the upper mantle density structure beneath the Rungwe Volcanic Province (RVP) and surrounding regions derived from a constrained inversion of satellite gravity anomalies with shear-wave velocity perturbations from Emry et al. (2019) as a reference model. For the inversion, we used the mutual information coupling technique in which a joint probability distribution of the gravity data and the reference model is estimated and its entropy minimized in order to generate a density model that is statistically compatible with the reference model. Our results reveal a conspicuous negative density anomaly beneath the Niassa Craton and less pronounced negative density anomalies beneath the RVP. We further determined a 3D model of the velocity-to-density conversion factor (f) and found negative f-values beneath the Niassa Craton which suggest the observed negative density anomaly is mostly due to compositional variations. Apart from the Niassa Craton, there are generally positive f-values in the study area, which suggest dominance of temperature control on the density structure. Although the RVP generally shows negative density anomalies and positive f-values, at shallow depths (<120 km) the f-value is near zero which suggests important contributions of both temperature and composition on the density structure beneath the RVP. We suggest that the presence of plume material at shallow depths beneath the RVP thermally reduces the upper mantle density. The plume material contributes metasomatic fluids that precipitate dense minerals that slightly increase the density due to compositional variations. The negative buoyancy of the Niassa Craton contributes to its long stability, while constituting a barrier to the southward flow of plume material, thus restricting the southward continuation of magmatism in the Western Branch of the EAR. The presence of negative-density anomaly where d l n v s are positive is incompatible with models based on the use of simple d l n v s to density conversion factors. These results have implications on how d l n v s models are converted to density perturbations.