The Mantle Viscosity Structure of Venus

The long‐wavelength gravity and topography of Venus are dominated by mantle convective flows, and are hence sensitive to the planet's viscosity structure and mantle density anomalies. By modeling the dynamic gravity and topography signatures and by making use of a Bayesian inference approach, we investigate the viscosity structure of the Venusian mantle by constraining radial viscosity variations. We performed inversions under a wide range of model assumptions that consistently predicted the existence of a thin low‐viscosity zone in the uppermost mantle. The zone is about 235 km thick and has a viscosity reduction of 5–15 times with respect to the underlying mantle. Drawing a parallel with the Earth, the reduced viscosity could be a result of partial melting as suggested for the origin of the asthenosphere. These results support the interpretation that Venus is a geologically active world predominantly governed by ongoing magmatic processes.

2 of 11 Kiefer & Peterson, 2003). Meanwhile, regional gravity and topography analyses showed that several Venusian highlands, mostly the plateaus associated with tessera terrains, were compensated by crustal thickness variations in contrast to the volcanic rises that have important support from deep mantle sources (Grimm, 1994;Maia & Wieczorek, 2022;Simons et al., 1997;Smrekar & Phillips, 1991).
Several studies have tested the impact of radial mantle viscosity variations on the predicted gravity and topography of Venus, either adopting a dynamic loading model (Herrick & Phillips, 1992;Kiefer et al., 1986;Pauer et al., 2006;Steinberger et al., 2010) or making use of 3D thermal evolution models (Huang et al., 2013;Rolf et al., 2018). The vast majority of these studies focused on the possibility of a viscosity jump at a depth analogous to the 660 km phase transition on Earth, which corresponds to about 730 km on Venus (Armann & Tackley, 2012), and found that the existence of such a feature was inconsistent with the gravity and topography observations. Alternatively, Pauer et al. (2006) made use of Monte Carlo inversions along with the dynamic loading model to estimate the viscosity structure of Venus. Their study showed that the viscosity of the mantle likely increases gradually with depth, and that there could be a thin low viscosity channel in the upper mantle. The moment of inertia and k 2 Love number of Venus could be used to investigate the viscosity profile as well, but they are not known with sufficiently accuracy to well-constrain relative variations with depth (see Figure 6 of Saliby et al., 2023).
In this work we used state-of-the-art inversion methods and data analysis techniques to constrain the mantle viscosity structure of Venus. We adopted the multitaper spatio-spectral localization method of Wieczorek and Simons (2007) to remove shallowly compensated regions from the analysis, and the viscosity estimations were performed using a Bayesian inference approach (Speagle, 2020). A variety of assumptions concerning boundary conditions and the density variations in the mantle were tested to assess the robustness of our results. In particular, we investigated scenarios where the density anomalies are concentrated within a single thin layer at a specific depth (e.g., Herrick & Phillips, 1992) or where they are uniformly distributed with depth in the mantle (e.g., Pauer et al., 2006). Ultimately, our study aims to contribute to a better understanding of the geodynamic and tectonic regime of Venus (e.g., Rolf et al., 2022) and to elucidate how the geologic histories of Earth and Venus diverged.

Dynamic Loading Model
The earliest seismic tomography studies of Earth showed that lateral variations of mantle temperature are correlated with long-wavelength geoid anomalies (e.g., Dziewonski et al., 1977). These observations motivated a series of studies to develop a dynamic loading model which allows for a quantitative interpretation of gravity in terms of mantle dynamics (Hager & Clayton, 1989;Ricard et al., 1984;Richards & Hager, 1984). The model consists of computing the instantaneous viscous flow that is predicted by the imposed density anomalies in the mantle for a given viscosity structure. The flow induces deformations of the planet's surface and core, generating dynamic topography and associated gravity anomalies.
For this model, the mantle is treated as an incompressible Newtonian fluid whose viscosity varies only with depth. Moreover, time-dependent and inertial forces are neglected. With these conditions, and by making use of a spherical harmonic decomposition of variables with an angular dependence, the problem can be written as a system of ordinary differential equations for each spherical harmonic degree which can be solved analytically using a propagator matrix technique. The solution is propagated from the core-mantle boundary to the surface with defined boundary conditions, passing through an arbitrary number of constant viscosity layers. We assume a free-slip boundary condition at the core-mantle boundary, whereas for the surface we evaluate both no-slip and free-slip end-member cases.
The model consists of developing depth-and viscosity-dependent kernels as a function of spherical harmonic degree, which describes how the planet adjusts to a unitary mantle load at radius r. In particular, we are interested in the gravity kernel G ℓ and topography kernel H ℓ , and these are computed following the approach of James et al. (2013) that includes an elastic lithosphere (see Section S1 in Supporting Information S1). Once the kernels have been computed, we can estimate the dynamic gravity and topography as the convolution of the kernel with the imposed density anomalies in the mantle δρ ℓm (r), as follows: where R cmb is the radius of the core-mantle boundary, R is the mean planetary radius and l and m are respectively the spherical harmonic degree and order. Figure S1 in Supporting Information S1 demonstrates how the kernels depend upon the viscosity profile and depth of the mass anomalies in the mantle.
Since tomography models are not available for Venus, we treat the mantle density anomalies as an unknown that will be determined in our inversion procedure. To make the inversion tractable, we either assume that the anomalies are concentrated in a single thin layer at a specified depth (Herrick & Phillips, 1992;James et al., 2013), or that the density anomaly has the same value at all depths (Kiefer et al., 1986;Pauer et al., 2006). The predicted gravity for the single mass-sheet is with ϕ ℓm representing the surface density anomaly (in kg m −2 ) at radius R ϕ . The equation for the case where the density anomaly is constant with depth is given in Section S2 in Supporting Information S1. Following the approach of Pauer et al. (2006), each coefficient ϕ ℓm is computed such that it minimizes the difference between the observed and predicted gravity and topography (see Section S3 in Supporting Information S1).

Localized Bayesian Inversion
Although mantle flows play an important role in shaping the long wavelength gravity and topography of Venus (e.g., Kiefer et al., 1986;Phillips & Malin, 1983), there are major topographic highlands that are mainly a result of crustal thickness variations (e.g., Kucinskas et al., 1996;Maia & Wieczorek, 2022;Simons et al., 1997). These shallowly compensated regions are inconsistent with the assumptions of the global dynamic loading model and Pauer et al. (2006) found that the worst predictions from their inversions were for the highlands of Ishtar Terra and Ovda Regio. They attempted to remove these signals by applying a binary mask to the gravity and topography, and then computing a localized power spectrum, but binary masking procedures have well known spectral leakage problems (e.g., Wieczorek & Simons, 2005).
In order to more rigorously remove from our analysis the signals associated with the compensated highlands, we employ the multitaper spectral analysis technique as developed by Wieczorek and Simons (2005); Wieczorek and Simons (2007). Our analysis region excludes Ishtar Terra and Western Aphrodite Terra (see Figure 1), and following the approach of Simons et al. (2006) we constructed orthogonal localization windows using a specified spectral bandwidth. For a bandwidth of ℓ win = 3, there are a total of 9 windows that concentrate more than 99% of their power in the region of interest. This number of windows provides an acceptable uncertainty for the localized spectra, and the small spectral bandwidth provides a large number of uncorrelated spectral estimates (see Section S6 in Supporting Information S1).
The results of the multitaper localization are shown in Figure 1. We make use of the VenusTopo719 topography model (Wieczorek, 2015a) and the MGNP180U gravity solution (Konopliv et al., 1999), both of which are based on the final Magellan mission datasets. The map in panel (a) shows the total power of the 9 localization windows summed in the space domain with the target localization region outlined by the white contour. In panel (b) we present the global and localized spectral admittance and correlation of gravity and topography (see Section S6 in Supporting Information S1 for the definition of these quantities). The localization leads to an increase in the admittance of about 30% over the entire spectrum, which is caused by the exclusion of highland regions that have high topography and low gravity. The correlation also shows a significant increase in the long-wavelength range due to the data localization, for ℓ < 40 the average correlation increases from 0.81 to 0.89.
The localized spectra of gravity, topography, and admittance (shown in Figure S2 in Supporting Information S1) are then used to invert for the mantle viscosity structure. These observations are compared with similarly localized spectra predicted by the dynamic loading model. One important aspect of the model is that the predicted gravity and topography are only sensitive to relative viscosity variations, and that the absolute viscosity cannot be 4 of 11 constrained. When considering the case where the mantle density anomalies are modeled as a single mass-sheet (Equation 3), the depth of the density anomaly is taken as a free parameter. The model depends on other parameters such as the core radius, the core-mantle density contrast and the average elastic lithosphere thickness, which are fixed in this study (see Table S1 in Supporting Information S1). The mantle and core related parameters are from Aitta (2012) and are based on the Venus-scaled preliminary reference Earth model. The adopted values are consistent with moment of inertia estimates (Margot et al., 2021) and other interior modeling studies (Dumoulin et al., 2017) of Venus. The elastic thickness of the lithosphere was set to zero. As discussed in Section S5 in Supporting Information S1 and shown in Figures S3-S5 in Supporting Information S1, the choice of these parameters has only a negligible impact on our results. To statistically evaluate the uncertainties of our model estimations, and considering the relatively large number of free parameters in our problem, we opted for a Bayesian sampling technique that provides the posterior probability of each parameter. We made use of the DYNESTY package, which is a Python implementation of the dynamic nested sampling method (Speagle, 2020). Nested sampling estimates the marginal likelihood and the posterior distribution by sampling within nested shells of increasing likelihood. The likelihood function adopted in our inversions is described in Section S4 in Supporting Information S1.

Inversion Results
Even though we performed inversions for a wide range of scenarios (see Section S5 and Figures S4 and S5 in Supporting Information S1), we chose to focus our analysis on the three cases that presented the largest variations in the results. The nominal case has a no-slip boundary condition at the surface and the mantle density anomaly is parameterized by a single mass sheet. The free-slip case differs from the nominal by having a free-slip boundary at the surface that allows for tangential movement of the surface, while the δρ-constant case has constant density anomalies with depth in the entire mantle along with a no-slip boundary at the surface. For these three inversions the number of constant-viscosity layers was set to four, with each layer being specified by its viscosity η i and depth to the bottom of the layer d i . Since our model is only sensitive the relative viscosity variations, we set η 1 = 1. Assuming that the core radius is known, the viscosity structure is defined by six free parameters. The nominal and free-slip scenarios have the mass-sheet depth d ϕ as an additional free parameter. Given the lack of information about our parameters, we considered for our priors a uniform distribution for the depth-related parameters and a log-uniform distribution for the viscosities. The only strong prior we set was to assume that the viscosity of the uppermost lithospheric layer was greater than the underlying layer (i.e., log 10 (η 2 /η 1 ) < 0). Such an assumption is a natural feature of temperature dependent rheological models (e.g., Breuer & Moore, 2015) and was used previously by Pauer et al. (2006). Figure 2 presents the posterior probability distribution of each free parameter for our three scenarios. The upper four panels are for the depths of the first three viscosity layers and depth of the mass-sheet. The bottom panels correspond to the viscosity ratios of the second, third, and fourth layers with respect to the overlying layer. Positive ratios indicate an increase in viscosity with respect to the layer above while negative ratios indicate a decrease in viscosity. All parameter estimations are detailed in Table S3 in Supporting Information S1.
Our results show that all scenarios consistently prefer shallow depths for the base of the lithospheric layer, with values less than about 200 km (Figure 2a). In contrast, the viscosity decrease to the underlying layer ( Figure 2e) is relatively unconstrained, which is partially a result of the small thickness of the uppermost layer. The viscosity interface between the second and third layer is the best constrained from our inversions (panels 2b and 2f). All three model scenarios indicate that the third layer increases in viscosity by about one order of magnitude at a depth of about 245-435 km, with median values for the mass-sheet depth of 239 km for the nominal case and 329 km for the free-slip case. The change in viscosity between the third and fourth layers is quite variable and differs in all three loading scenarios (panels 2c and 2g). For the nominal case, both the layer depth and viscosity are poorly constrained, although the solutions prefer depths larger than 1,000 km. The free-slip case tends to prefer larger depths and lower viscosities for the deepest layer, with the most probable case corresponding to no change in viscosity. The δρ-constant model, on the other hand, has a well-constrained viscosity increase of about 10 times at 1,300-1,550 km depth.
In Figure 3 the viscosity profiles for all solutions are shown via 2-dimensional posterior distributions for our three loading scenarios. Since one of the best constrained aspects of our model is the increase in viscosity between the second and third layer, for a better visualization, we scaled our viscosity profiles such that the viscosity of the second layer was 1. The solid curves in these figures represent the logarithmic mean viscosity at each depth. The upper mantle structure is similar for the three scenarios, with a consistent viscosity increase occurring between the second and third layers. As for the lower mantle, the loading scenarios that use a single mass-sheet (Figures 3a and 3b) indicate an isoviscous structure, although some of the solutions suggest a basal low-viscosity layer above the core. As for the δρ-constant model (Figure 3c), we see a second viscosity jump at about 1,400 km depth.

Upper Mantle Low Viscosity Zone
Our inversions indicates the presence of a zone beneath the lithosphere characterized by viscosity values that are roughly 10 times lower than the underlying mantle. Its thickness is about 150-300 km, starting at the base of the lithosphere down to a depth of 268-435 km. This low viscosity zone can be interpreted as an asthenosphere-like layer. The asthenosphere of Earth is a mechanically weak layer starting beneath the lithosphere and that extends to the top of the transition zone at about 400 km depth. It is considered to be a key ingredient for plate tectonics (e.g., Rolf et al., 2022) and its existence has been supported by several geophysical methods. On Earth, the region is characterized by high electrical conductivity, low seismic velocities, and strong seismic attenuation (e.g., Shankland et al., 1981). In oceanic regions, seismological observations have shown that the lithosphere-asthenosphere boundary occurs sharply at about 70 km depth, while for the sub-continental mantle the seismic signature of this boundary is weaker and deeper, at about 200 km depth (Karato, 2012).
Gravity investigations considering postglacial rebound and/or dynamic loading commonly indicate that the asthenosphere of Earth is also characterized by low viscosities, although the published estimates present some discrepancies (see reviews by King, 2016;Richards & Lenardic, 2018). Estimations of the low viscosity zone (LVZ) are generally associated with reductions in viscosity of one to three orders of magnitude with respect to the underlying mantle. Some studies indicate that the LVZ is fully contained within the asthenosphere (Forte & Mitrovica, 2001;Hager & Clayton, 1989), others suggest that the zone extends to the base of the upper mantle at 660 km depth (King & Masters, 1992;Liu & Zhong, 2016), and others find a low viscosity asthenosphere along with a thin low-viscosity channel at the 660 km transition (Mitrovica & Forte, 2004). Several factors could be responsible for these differences in interpretation. In particular, there is a well-known trade-off between the thickness and viscosity of the low-viscosity zone (Richards & Lenardic, 2018) and gravity studies have difficulties in accounting for lateral viscosity variations. Strong lateral viscosity variations could exist on Earth as a result of subducted slabs and differences in mantle structure beneath oceanic and continental crust (Čadek & Fleitout, 2003). However, such variations are likely to be less important on Venus given its lack of plate tectonics and different mode of heat transport, which is probably associated with regional-scale delamination scattered throughout the planet (e.g., Davaille et al., 2017;Gülcher et al., 2020;Lourenço et al., 2020;Smrekar & Stofan, 1997).
Gravity and topography studies of Venus have mostly argued against the existence of a low-viscosity zone in the upper mantle (e.g., Nimmo & McKenzie, 1998). However, most of those studies either did not perform inversions and limited the analysis to a few models representative of Earth-like scenarios (e.g., Herrick & Phillips, 1992;Kiefer et al., 1986;Kiefer & Hager, 1991;Steinberger et al., 2010). Moreover, gravity investigations prior to 1992 used data from the Pioneer Venus mission whose resolution was limited to degree 18. From a different perspective, studies by Huang et al. (2013) and Rolf et al. (2018) estimated geoid anomalies for Venus using three-dimensional mantle convection models. They showed that the presence of a thick LVZ, ranging from the base of the lithosphere down to the ringwoodite-bridgmanite phase transition at 730 km depth, was inconsistent with the observations. These results are in fair agreement with our study that predicts a LVZ thickness that is about half of the upper mantle thickness. In addition, we note that Monte Carlo inversions performed by Pauer et al. (2006) showed that many models were consistent with the presence of a LVZ with a thickness of a couple hundred kilometers.
The origin of such a low viscosity layer on Venus is arguable. As a starting point, we may draw a parallel with Earth and evaluate the mechanisms that have been proposed to explain the existence of its asthenosphere. However, this is a heavily debated topic with several proposed hypotheses. Some studies proposed that the asthenosphere results from the presence of small amounts of partial melts in the upper mantle (Anderson & Spetzler, 1970;Chantel et al., 2016;Debayle et al., 2020;Hua et al., 2023), while others consider that the region is better explained by a subsolidus regime associated with rheological weakening of mantle rocks under temperatures that are close to the solidus (Karato, 2012;Takei, 2017). For the latter hypothesis, it has been suggested that dissolved water in olivine could effectively reduce the viscosity (Hirth & Kohlstedt, 1996). Alternatively, low viscosities in the uppermost mantle could be caused by the predominance of dislocation creep, with larger depths being dominated by diffusion creep with higher viscosities (Semple & Lenardic, 2021;Van Den Berg & Yuen, 1996). Finally, viscosity interfaces in the mantle could be linked to mineralogic phase transitions (Meade & Jeanloz, 1990).

of 11
Even though the surface and atmosphere of Venus are dry, the water content of its interior is unknown. Both Venus and Earth should probably have accreted similar amounts of volatiles (e.g., O'Brien et al., 2018), but the abundance of these volatiles in Venus's interior is poorly known (e.g., Gillmann et al., 2020Gillmann et al., , 2022Way & Genio, 2020). The existence of partial melt in Venus's mantle seems plausible, given the indications that the planet is still volcanically active, particularly in regions associated with active mantle plumes (e.g., Gülcher et al., 2020;Herrick & Hensley, 2023;Mueller et al., 2008;Shalygin et al., 2015;Smrekar et al., 2010). Moreover, recent studies have shown that large amounts of magmatic intrusions could play a primary role in the mobility of the lithosphere and crustal recycling, representing an efficient mechanism for heat loss (Lourenço et al., 2020;Smrekar et al., 2023;Tian et al., 2023). Lastly, tomography investigations by Debayle et al. (2020) indicate that hotspots on Earth are associated with high melt content in the upper mantle.
Following the well-established experimentally-derived relation between melt fraction and viscosity (e.g., Hirth & Kohlstedt, 2003) we estimate that the low viscosity layer we find on Venus could be associated with 5%-11% melt. On the other hand, more recent experiments have shown that even very small interconnected fractions of melt can have a significant impact on the viscosity (Holtzman, 2016;Takei & Holtzman, 2009). From this perspective, our viscosity reduction estimations could be associated with as little as 0.05% melt (see Figure S6 in Supporting Information S1 for more details concerning these calculations). Finally, we note that the low viscosity zone and the possibly associated partial melt does not necessarily need to be uniform throughout the planet. In fact, considering that the largest gravity and topography signatures come from the volcanic rises it is likely that their signatures dominate our analysis and that our results are mostly representative of these geologically active regions.

Consequences of the Mantle Load Parameterization
Due to the lack of information regarding density heterogeneities in the Venusian mantle we were required to make some simplifications in our inversions. To assess the importance of these assumptions, we investigate further the two density anomaly parameterizations used in our study. Herrick and Phillips (1992) argued that since mantle plumes likely dominate the density anomaly distribution in Venus, most of the anomalies should be concentrated in a relatively thin and horizontal layer associated with the plume head beneath the lithosphere. In this scenario, a single mass-sheet parameterization would be more appropriate. On the other hand, Pauer et al. (2006), made a comparison to Earth's mantle density patterns arguing that plumes and subducted slabs penetrate the mantle more or less vertically, indicating that a depth-independent distribution of density anomalies could be a reasonable first approximation.
Although the two different scenarios have a significantly different depth-dependence of the density anomalies, our results consistently predict a low viscosity zone in the upper mantle. For the deep mantle, however, clear differences in the viscosity structure were obtained. The single mass-sheet scenarios (Figures 3a and 3b) tends to accept a wide range of solutions for the deepest viscosity layer, with a large number of models consistent with an isoviscous mantle below the low viscosity zone. Nevertheless, we note that about 35% of the models present a viscosity decrease of over one order of magnitude from layer 3 to 4, which would correspond to a basal low viscosity zone. If confirmed, these regions could be analogous to the large low shear velocity provinces on Earth (e.g., French & Romanowicz, 2015) or indicate the presence of partial melting (O'Rourke, 2020). In contrast, the depth-independent case (Figure 3c) requires an increase in viscosity in the mid mantle between the third and fourth layers at 1,330-1,550 km depth, which is significantly deeper than the ringwoodite-bridgmanite phase transition. Interestingly, a comparable viscosity jump has been suggested for Earth by several dynamic loading investigations (e.g., Forte & Peltier, 1991;Rudolph et al., 2015Rudolph et al., , 2020, indicating an increase in viscosity at depths of about 800-1,200 km that corresponds roughly to 960-1,330 km on Venus. In any case, our results should motivate future studies to explore more realistic density anomaly distributions in Venus based on mantle convection simulations, or by using a more complex statistical distribution of mass anomalies as in Steinberger et al. (2010).
The predicted density anomaly distribution estimated for the nominal loading model with the largest likelihood is shown in Figure S7 in Supporting Information S1. As expected, at large volcanic rises, such as Atla and Beta Regiones, we observe negative density anomalies, associated with positive buoyancy and commonly interpreted as regions of mantle upwellings, hotspots and active volcanism (e.g., James et al., 2013;Kiefer & Hager, 1991;Smrekar et al., 2010;Stofan et al., 1995). These anomalies are generally also correlated with regions of high heat 9 of 11 flow (Smrekar et al., 2023) and with regions where coronae could be still active (Gülcher et al., 2020). On the other hand, positive density anomalies correlate with volcanic plains which can interpreted as regions of mantle downwellings. These results are consistent with the density anomaly distributions found in previous studies (Herrick & Phillips, 1992;James et al., 2013;Pauer et al., 2006). Moreover, we note that the different loading scenarios investigated here have comparable horizontal density distribution patterns, although for the δρ-constant case the shorter wavelengths present relatively higher amplitudes (Herrick & Phillips, 1992).

Conclusions
Our study used a Bayesian approach to investigate the mantle viscosity structure of Venus. We employed a dynamic loading model to predict the planet's long-wavelength gravity and topography and compared these predictions to the observations. Using a range of model scenarios, we consistently found that Venus presents a low viscosity zone in the uppermost mantle, a layer that could be interpreted as an Earth-like asthenosphere, potentially resulting from partial melt in the upper mantle. This interpretation supports previous studies that proposed that Venus is currently an active volcanic world (e.g., Gülcher et al., 2020;Herrick & Hensley, 2023;Rolf et al., 2022;Smrekar et al., 2010). Moreover, our inversions disfavor a viscosity jump associated with the ringwoodite-bridgmanite phase transition at 730 km depth.
One aspect that merits further investigation is how more realistic distributions of density anomalies in the mantle would affect the predicted viscosity profile. For example, one could use geodynamic simulations to estimate the density distribution based on temperature anomalies predicted by these models. Improvements on deep interior constraints for Venus could be achieved in the future by coupling dynamic loading to tidal deformation investigations. In particular, an integrated analysis would allow for an assessment of the absolute mantle viscosity profile. This approach will be particularly powerful once future missions obtain precise estimates of tidal Love numbers, tidal quality factor, and moment of inertia factor (Cascioli et al., 2021;Rosenblatt et al., 2021).

Data Availability Statement
A Python routine to estimate the gravity and topography predicted by the dynamic loading model can be found in Maia (2023). The localized spectral analysis was performed using the open-source package Pyshtools (Wieczorek & Meschede, 2018), while the posterior parameter distributions were estimated using the Dynesty package (Speagle, 2020). The spherical harmonic model of the gravity field used in this study can be found in Sjogren (1997), and the topography data set is from Wieczorek (2015b). The perceptually uniform color maps used in this work are from Crameri (2018).