Corona structures driven by plume–lithosphere interactions and evidence for ongoing plume activity on Venus

In the absence of global plate tectonics, mantle convection and plume–lithosphere interaction are the main drivers of surface deformation on Venus. Among documented tectonic structures, circular volcano-tectonic features known as coronae may be the clearest surface manifestations of mantle plumes and hold clues to the global Venusian tectonic regime. Yet, the exact processes underlying coronae formation and the reasons for their diverse morphologies remain controversial. Here we use three-dimensional thermomechanical numerical simulations of impingement of a thermal mantle plume on the Venusian lithosphere to assess the origin and diversity of large Venusian coronae. The ability of the mantle plume to penetrate into the Venusian lithosphere results in four main outcomes: lithospheric dripping, short-lived subduction, embedded plume and plume underplating. During the first three scenarios, plume penetration and spreading induce crustal thickness variations that eventually lead to a final topographic isostasy-driven topographic inversion from circular trenches surrounding elevated interiors to raised rims surrounding inner depressions, as observed on many Venusian coronae. Different corona structures may represent not only different styles of plume–lithosphere interactions but also different stages in evolution. A morphological analysis of large existing coronae leads to the conclusion that at least 37 large coronae (including the largest Artemis corona) are active, providing evidence for widespread ongoing plume activity on Venus. Thermomechanical modelling shows that the formation and diverse morphologies of coronae on Venus can be explained by interactions between the lithosphere and impinging mantle plumes. Some corona structures are consistent with ongoing plume activity.

in the models since estimates vary on Venus 14,39,40 . Results of the numerical experiments (Extended Data Table 2) show systematic development of corona-like structures and reveal four regimes of plume-lithosphere interaction: (i) lithospheric dripping, (ii) ephemeral subduction, (iii) embedded plume and (iv) plume underplating. As summarized in Figure 3, these regimes are strongly dependent on plume buoyancy and the lithospheric and crustal configuration.  (Fig. 4c), the plume spreads radially outwards, thereby broadening the developing corona structure at the surface. The crust and lithospheric mantle initially located above the plume are further displaced and thickened. Being rheologically weak, they form a viscous lithospheric drip at the plume margins. The corresponding topographic relief is characterized by a deep annular trough and an elevated interior that slowly loses its elevation. As crustal material is brought to greater depths and is subjected to higher temperatures, it densifies due to the basalt to eclogite phase change (Methods, Extended Data When the same plume pierces through a lithosphere with a twice thinner crust (and hence lower moho temperatures, Figs. 5a-c), less basaltic material is transformed into eclogite. Instead of lithospheric drips, a radial subduction zone develops at the plume margin. This subduction zone is short-lived and followed by circular slab-detachment (Fig. 5b). The topographic evolution is similar to that of the reference lithospheric dripping model (Fig. 4), but bending of the lithosphere into down-dipping slabs produces an outer rise with higher amplitude. A similar subduction-like scenario is obtained in the model with an older (80-120 Myr) and thus colder lithosphere that promotes a more rigid behavior (models M10-M14 in Extended Data Table 2).
Lithosphere rigidity results in one or several episodes of retreating subduction that may be focused towards a certain arc, and are always followed by segmented slab break-off. As in the case of lithospheric drips, these models eventually isostatically relax to form an inverted topographic profile with an outer rim and inner depression.
The ratio of plume buoyancy over lithospheric strength 38,41 exerts a major control on plume penetration in the Venusian lithosphere (Methods, Fig. 3). When the lithosphere is sufficiently cold (Tmoho < 1100 K) and the plume is not buoyant enough relative to the lithospheric strength (Buoyancy ratio Br < 4), the plume cannot reach and pierce through the lithosphere but spreads laterally under it (plume underplating regime, Figs. 5d-f). The lithosphere is thinned above the plume and the surface displays only a dome or plateau, without any rim or trench.
Even after isostatic relaxation, the interior of the corona remains elevated. A transitional "embedded plume" regime was also detected in several experiments, in which the mantle plume partially penetrates the lithosphere but does not force crustal material to be recycled down into the mantle (Extended Data Fig. 4). The embedded plume induces magma emplacement in the crust followed by crustal flow, thereby resembling some of the mechanisms seen in prior smallscale crustal convection models of novae and coronae 26 .
Model limitations (Methods) prevent a detailed comparison with observed tectonic structures, which may also be influenced by fine-scale phenomena such as dike intrusion and elastic loading. Nevertheless, most of the models show surface strain-rate patterns consistent with concentric deformation structures at the corona margins (Extended Data Figure 5). In the most rigid models, stellate and radial surface deformation patterns can also be distinguished.
Volcanism is present in all plume penetration models as newly-formed crust (volcanic material) is emplaced in the corona interior via melt percolation (Figs. 4, 5a-c and Extended Data Fig. 4).

Morphological diversity of Venusian coronae
The temporal evolution of the radially averaged topographic profiles for models representative of each regime shows how different corona morphologies not only represent different plume-lithosphere interactions, but also different stages in evolution (Fig. 6). Our models suggest that a penetrating plume is key for lithospheric dripping (high moho temperatures) or lithospheric subduction (low moho temperatures) at the plume margin, and for lithospheric thinning of the corona interior (Extended Data Fig. 6). The suction below the lithospheric drips or slabs and bending of the lithosphere creates a circular trench and an outer rise at the surface, that are inverted during the isostatic relaxation stage (Figs. 6a-b). In our models, coronae with rims and trenches are only produced by mantle plumes that at least partially penetrate into the Venusian lithosphere. Therefore, the common occurrence of coronae displaying rims 24,25 suggests that most plumes that formed coronae were able to penetrate at least partially into the Venusian lithosphere. In addition, for coronae formed by a penetrating mantle plume, we are able to distinguish active from inactive structures: active coronae feature an outer trench and rise that imply ongoing suction above downwards-moving lithospheric material and an elevated interior supported by plume buoyancy, whereas inactive coronae show an inverted topographic profile of outer rim and inner depression linked to a thinned lithosphere ( Fig. 6).
Observational evidence for coronae activity Based on our models, we can assess the activity of individual coronae on Venus, assuming that mantle-lithosphere interactions are dominantly responsible for their formation. Most coronae formed in the numerical models range in diameter from 300 to 1000 km (Extended Data Table 2), encompassing a wide range of corona dimensions that are observed on Venus 23,42 . A direct comparison of the surface topography of representative coronae with those obtained in selected models is shown in Figure 1 and Extended Data Figure 7. Thouris Corona (6.5°S; 12.9°E) and Aruru Corona (9.0°S 262°E) display raised rims surrounding lower elevations, features that mark later corona evolution stages, implying that the underlying plume is now inactive (Figs. 1bc). In contrast, the outer rise and deep trench displayed by Aramaiti Corona (25.5°S 82°E, Fig. 1a) is evidence of lithospheric bending and suction above downgoing lithosphere, implying that the coronae is currently in an active stage that predates final isostatic relaxation into an inverted topographic profile. Based on these topographic characteristics, we suggest that the impressive Artemis corona (25°S, 135°E) also hosts ongoing activity at its south-eastern margin (Extended Data Fig. 7), which is in agreement with former interpretations of slab-retreat at the corona margin 29,30,31,43,44 .
We evaluated the possibility of present-day corona activity by systematically analyzing the surface topography of large coronae (diameter > 300 km) on Venus (see Methods for details regarding this analysis). A corona is labeled as currently "active" if it features a clear outer rise and trench, "inactive" if no outer rise but a rim and inner depression is evident, or "unclassified" if the presence of these features is ambiguous. The results of this global analysis are summarized in Fig. 2 19,46 , which would be consistent with current activity, but lack the topographic trench characteristics of dripping or subducting lithosphere (Extended Data Table 1). Lada Terra, where many of our active coronae are located, and Quetzalpetlatl Corona (68°S,357°E) have previously been proposed sites of recent volcanic formation. If regions of thin elastic lithospheric thickness 14 such as the Themis, Alpha and Eistla Regio, also have thin thermal lithospheres, they would be particularly prone to plume penetration. Indeed, several of our identified locations of enhanced activity also feature thin elastic thicknesses (Fig.2). Conversely, almost no (active) coronae are identified within the hightopography regions of the Aphrodite and Ishtar Terra. A low apparent elastic thickness there may indicate isostatic compensation, not thin thermal lithosphere 14,15 .
So far, we have been able to ascertain ongoing activity or inactivity for relatively few coronae on Venus. This is due not only to the dimensions and resolution of our numerical models (Methods), but also to the scarcity of high-resolution topography data ("Poorly resolved" coronae in Extended Data Table 1). Improvements upon numerical techniques and resolution are needed to allow for better comparisons with different types of datasets (such as gravity, heat emission, lava flows and imagery). Furthermore, any future mission that would collect more and higher

Additional information
Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints.

Competing financial interests
The authors declare no competing financial interests.

Modelling approach
A modelling tool designed to investigate the dynamics involved with plume-lithosphere interactions needs to incorporate a thermo-rheologically realistic lithosphere fully coupled to mantle dynamics in three dimensions. We meet this challenge using the staggered grid/particlein-cell viscous-plastic 3D code I3ELVIS 53,54 . This parallel implicit multi-grid code is based on a combination of the finite difference method, applied on a fully staggered Eulerian grid and a marker-in-cell technique 55 . The momentum, continuity, and energy equations are solved on the Eulerian grid, whereas physical properties are transported by Lagrangian markers that move according to the velocity field interpolated from the fixed grid. Non-Newtonian viscous-plastic rheologies are used in the model (Extended Data Table 1), which is also fully thermodynamically coupled and accounts for mineralogical phase changes, adiabatic, radiogenic and frictional internal heating sources. The model takes into account melt extraction and upward transport from the plume, rheological weakening of the lithosphere subjected to melt percolation, crustal growth by magmatic processes, and the basalt to eclogitic phase change to capture the essential geophysical plume-lithosphere interaction processes (see below) Full details on the method are provided in prior studies 55 and below. This algorithm has been thoroughly tested, both in two and three dimensions, and used for simulating mantle plume behavior and lithospheric deformation experiments in various previous studies 26,38,41,57 .

Reference model design
The reference model (M0) set-up is shown in Extended Data where is the total rock pressure, is the internal friction coefficient with ! as the initial friction coefficient, taken to be 0.2 for the weak Venusian rheology (following approaches of Gerya et al. 26,55 ), ≥ 0 the time-integrated plastic strain that is integrated as long as plastic deformation is active ( ! = 0.5 is the upper strain limit for fracture-related weakening), is the compressive rock strength at = 0, is time in seconds, ̇: ;(C02-3*1) is the plastic strain rate tensor, and @A03 A++ is the long-term melt-induced weakening factor 56 defined by: where @A03,! is the standard value for melt-induced weakening, @A03 is the local melt flux, and Rocks are considered non-molten (refractory) when the extracted melt fraction is larger or equal than the standard one (∑ AH3 I = ! ). To simulate melt extraction, an extraction threshold @2H = 4 wt.% is defined. If the total amount of melt M surpasses this threshold ( > @2H ), all melt except a non-extractable @*/ = 2 wt.% is extracted: Melt percolation is assumed to be faster than the deformation of non-molten mantle 63 . Based on this assumption, the extracted melt AH3 is instantaneously removed from its source region, moved upwards vertically, and is added to the bottom of the crust as intrusive plutons or to the surface as volcanic rocks according to a predefined plutonic-to-volcanic rock ratio (80:20%).
The effective density of the mafic magma and molten crust is computed as 60 : where M = 1000 J/kg is the heat capacity of the solid crust and N = 380 kJ/kg is the latent heat of crystallization of the crust 59 .

Phase changes
Eclogitization of dripping or subducted crust is implemented as a linear density increase with pressure from 0% to 16% in the pressure-temperature region of garnet-in and plagioclase-out phase transitions in basalt (Extended Data Fig. 3) 64 . Numerical implementations of slab dehydration and mantle hydration used in former numerical studies on 3D plume-lithosphere explorations 38,54 are deactivated in the code, so that our model represents a dry planetary interior appropriate for Venus.

Heat flux calculations
Vertical heat flows are calculated for each model just beneath the surface at a given depth. We choose a depth slightly below the surface ( -,5+ = 1 km) such to not be affected by possibly high thermal conductivity values for the sticky air. To do so, heat flow (in W/m 2 ) is calculated at the center of the neighboring cell (above and below this depth) using: with averaged thermal conductivity : for cell (equations for thermal conductivity are given in Extended Data Table 1), and the vertical thermal gradient The Eulerian grid used in most of the experiments is discretized in cells that represent a volume of 5 × 5 × 2 km 3 . Even though the many Lagrangian markers better resolve the transitions in rock types and physical properties than this resolution, the grid size precludes the analysis of small-scale localized shear zones. Higher resolution models would be needed to assess assumptions that have been made with the initial model set-up also promote a relatively rapid evolution of the models. First off, our set-up places a mantle plume right beneath the Venusian lithosphere (Extended Data Fig. 1), resulting in near-instantaneous crustal uplift. Secondly, the mantle potential temperature was not investigated in this study. Higher mantle potential temperatures lead to a compositionally more buoyant lithosphere 68 , yet also to a relatively less buoyant mantle plume, which might slow down model evolution. Finally, plume activity at the surface could be prolonged by the constant supply to the coronae by heat and magma coming from multiple plumes 69 or an active plume-tail 70 .

Comparison with natural data
Numerical results and visualizations thereof can be compared with available observational data of the Venusian surface. While a reasonable amount of topographic maps and radar images is and are hence used in previous scientific studies, which are listed in this database. In our global coronae analysis, we focus on coronae that are best resolved by the observable data and with comparable dimensions to those in our numerical experiments, and therefore analyze only coronae with a diameter of >300 km, encompassing ~40% of the database. In addition, we analyzed several coronae with smaller dimensions that are situated in formerly identified "hot spots" 16,19,46 .

Data availability
The numerical data that supports the findings of this study can be requested from the corresponding author. A KML file based on the coronae classification in this paper is available on

Code availability
The numerical code is available upon reasonable request. Requests can be made to T.V. Gerya.   Fig. 3 and Extended Data Fig. 4).