The Tectonics and Volcanism of Venus: New Modes Facilitated by Realistic Crustal Rheology and Intrusive Magmatism

To explain Venus' young surface age and lack of plate tectonics, Venus' tectonic regime has often been proposed to be either an episodic-lid regime with global lithospheric overturns, or an equilibrium resurfacing regime with numerous volcanic and tectonic activities. Here, we use global 2-D thermochemical convection models with realistic parameters, including rheology (dislocation creep, diffusion creep, and plastic yielding), an experiment-based plagioclase (An$_{75}$) crustal rheology, and intrusive magmatism, to investigate the tectonics and mantle evolution of Venus. We find that surface tectonics is strongly affected by crustal rheology. With a ''weak'' plagioclase-rheology crust, models exhibit episodic overturns but with continuously high surface mobility and high distributed surface strain rates between overturns, leading to a new tectonic regime that we name ''deformable episodic lid''. On the other hand, olivine-crustal-rheology models exhibit either standard episodic-lid tectonics, i.e. with mobility that is high during overturns and near zero between overturns, or stagnant-lid tectonics, i.e. with near-zero mobility over the entire model time. Also, a combination of plagioclase crustal rheology and dislocation creep can weaken the lithosphere sufficiently to facilitate lithospheric overturns without applying plastic yielding. Internally, the composition-dependent density profile results in a ''basalt barrier'' at the mantle transition zone, which strongly affects Venus' mantle evolution. Only strong plumes can penetrate this basalt barrier and cause global lithospheric overturns. This basalt barrier also causes global internal episodic overturns that generate global volcanic resurfacing in stagnant-lid models, which suggests a new resurfacing mechanism (we name it ''stagnant episodic-volcanic-resurfacing'') that does not involve lithospheric overturns.

with continuously high surface mobility and high distributed surface strain rates between overturns, leading to a new tectonic regime that we name "deformable episodic lid". On the other hand, olivine-crustal-rheology models exhibit either standard episodic-lid tectonics, i.e. with mobility that is high during overturns and near zero between overturns, or stagnant-lid tectonics, i.e. with near-zero mobility over the entire model time. In models with plagioclase-rheology crust, magmatic resurfacings are short-lived and randomly located, leading to a rather uniform resurfacing rate. Also, a combination of plagioclase crustal rheology and dislocation creep can weaken the lithosphere sufficiently to facilitate lithospheric overturns without applying plastic yielding. Internally, the composition-dependent density profile results in a "basalt barrier" at the mantle transition zone, which strongly affects Venus' mantle evolution. Only strong plumes can penetrate this basalt barrier and cause global lithospheric overturns. This basalt barrier also causes global internal episodic overturns that generate global volcanic resurfacing in stagnant-lid models, which suggests a new resurfacing mechanism that does not involve lithospheric overturns. We name this regime "stagnant episodicvolcanic-resurfacing". Results in models with plagioclase crustal rheology and low eruption efficiencies best fit current estimates of crustal thickness, surface age, and tectonic structures based on observations of Venus. The crustal thickness in models with an olivine-rheology crust is limited by the basalt-eclogite phase transition and is higher than observation-based estimates. Higher intrusion rates better fit these estimates, as models with high eruption efficiency predict very low surface age (< 100 Myr).

Introduction
Venus is of particular interest in planetary science because of its Earth-like size, density, distance from the Sun, and active, rocky surface, all of which suggest a comparable bulk composition between Venus and Earth. From images from NASA's Magellan mission, the observations of ∼1000 impact craters on the surface of Venus lead to an estimated surface age between 300 Myr and 1 Gyr (e.g. Hauck et al., 1998;McKinnon et al., 1997;Nimmo and McKenzie, 1998). Some studies suggest an even younger surface age of 150 -240 Myr (Herrick and Rumpf, 2011;Le Feuvre and Wieczorek, 2011) based on new analyses of the crater population and a revised crater chronology of the Solar System. These estimates suggest that Venus is the only terrestrial planet in the Solar System that has a resurfacing age comparable to that of Earth. However, the nature of the current tectonics of Venus that leads to a young and relatively uniform surface remains unclear.
The near random distribution of the impact craters on Venus' surface suggests a fairly uniform surface age for Venus, with the age of about 80% of the surface being not significantly different from the average surface age . Also most of these craters are not significantly altered by volcanism and tectonics (Strom et al., 1994). Such a uniform surface can hardly be explained by Earth-like plate tectonics. Thus, two end-member models have been proposed to explain Venus' tectonics: the catastrophic resurfacing model and the equilibrium resurfacing model. The catastrophic resurfacing model suggests that a global catastrophic resurfacing event occurred in a relatively short period of time and reworked most of the surface of Venus (e.g. Strom et al., 1994;Romeo and Turcotte, 2010). The equilibrium resurfacing model argues that the crater distribution results from numerous randomly distributed volcanic or tectonic events, and can also explain these cratering observations (e.g. Bjonnes et al., 2012). These two end-member models correspond to different evolution paths for volcanic and tectonic activities, which can be reflected in variations in the history of the resurfacing rate and surface mobility.
In the absence of plate tectonics, most of the plains on Venus' surface are believed to be resurfaced by volcanism (Stofan et al., 2005), which favors the equilibrium resurfacing model. Some observations indicate that volcanism could remain active until today, such as the high surface emissivity at topographic rises (Smrekar et al., 2010) and lava flows (Brossier et al., 2021), a locally thinned lithosphere at Aramaiti Corona (Russell and Johnson, 2021), and spikes of sulfur dioxide in the atmosphere of Venus (e.g. Marcq et al., 2013). Numerical models of coronae formation also suggest that some coronae could be currently active, possibly due to ongoing plume activities (Gülcher et al., 2020). However, although these volcanic features on Venus' surface show evidence of regional-scale resurfacing, the formation of large volcanic plains on Venus can hardly be explained by the equilibrium resurfacing models (Romeo and Turcotte, 2010).
The evolution of Venus' surface mobility is reflected by tectonic features.
From stratigraphic analyses of Venus' surface, the tectonic and volcanic units are interpreted to form at two different periods Kreslavsky et al., 2015), which suggests that Venus could have experienced both catastrophic and equilibrium resurfacing events. For most of these units, the crater densities are similar to the mean crater density of Venus' surface and are interpreted to form at a similar time with a global resurfacing event.
Some units, such as those related to the network rifting-volcanism regime in Ivanov and Head (2015) and the post-regional-plain units in Kreslavsky et al. (2015), are significantly younger. In addition, there is evidence of high lithospheric mobility preserved on Venus' surface, including tesserae (e.g. Hansen et al., 1999) and compressional structures around Lakshmi Planum (Harris and Bédard, 2015). Recent analysis of satellite images of Venus' lowland plains also suggests that there are fragmented lithospheric blocks that collide and rotate relative to each other, indicating possibly widespread surface mobility over time (Byrne et al., 2021) which is not consistent with a purely stagnant-lid regime. These observations led to the inference that Venus could have been operating in an intermediate regime between catastrophic and equilibrium resurfacing, with one global overturn followed by continuous regional-scale resurfacing and tectonic events.
These models often assume an olivine composition and diffusion-creep rheology for Venus' crust, and neglect intrusive magmatism. Furthermore, most of the two-dimensional (2-D) and three-dimensional (3-D) numerical models focused on the episodic-lid regime, with near-global lithospheric overturns and near-stagnant-lid phases between overturns. However, a near-stagnant-lid phase after the last global overturn cannot explain the presence of young tectonic units on Venus' surface. Pure stagnant-lid models also cannot explain Venus' surface age and crustal thickness, because when magmatism is purely eruptive, heat transfer is dominated by the magmatic heat-pipe mechanism, resulting in high volcanic resurfacing rates and very young surface ages that are inconsistent with estimates based on cratering density, as demonstrated by the stagnant-lid simulations of Armann and Tackley (2012). Recently, a new tectonic regime named "plutonic-squishy lid" was proposed by Lourenço et al. (2020). In this regime, intrusive magmatism warms, weakens and thins the lithosphere and facilitates surface deformation and mobility. Interpretations that Venus' lowland plains are fragmented and moving through time based on satellite images (Byrne et al., 2021), and that the lithospheric thickness and heat flow of Venus' coronae are similar to those at rifting zones on Earth , are more consistent with this tectonic regime than with episodic or stagnant-lid regimes. For a review of Venus' mantle dynamics and tectonic regime see Rolf et al. (2022).
Although it is generally accepted that Venus' crust is basaltic (e.g. Surkov et al., 1984), the rheology of this basaltic crust is uncertain. Basalt is mainly made up of two groups of minerals, plagioclase feldspar and pyroxene, which have quite different creep strengths. Compared to olivine, plagioclase is typically considered to be weaker, while pyroxene is typically considered to be stronger. Mackwell et al. (1998) measured the rheology of dry diabase (the same composition as basalt but larger grain size), and found that its strength is dependent on the plagioclase fraction, but in any case is likely to result in a high viscosity crust that is strongly coupled to the mantle. On the other hand, Azuma et al. (2014) argued that the crust of Venus is likely weaker than the mantle. The observations of crustal deformation discussed above, as well as the low to moderate (<40 km) effective elastic thickness for extensive regions of Venus' lithosphere (e.g. Jimenez-Diaz et al., 2015), are also evidence for a relatively weak crustal rheology for the planet. Using regional numerical models, Gerya (2014) and Gülcher et al. (2020) were able to obtain a good match to the observed coronae and novae structures using a plagioclase An 75 crustal rheology, rather than the olivine crustal rheology used by previous numerical models. Dislocation creep was also included in their models, which has a further weakening effect.
In this work, we aim to (1) constrain Venus' crustal rheological parameters by comparing the results of numerical experiments using different crustal rheologies to observations, and (2) understand how crustal rheology and intrusive magmatism affect Venus' tectonics and evolution.

Numerical Model and Boundary Conditions
The thermo-chemical evolution of the interior of Venus is modelled using the finite-volume code StagYY (Tackley, 2008) in a 2-D spherical annulus geometry (Hernlund and Tackley, 2008). The physical parameters and the numerical model are based on those of Armann and Tackley (2012) and Lourenço et al. (2020), with parameters listed in Table 1. The conservation equations for a compressible fluid with an infinite Prandtl number are solved on a fully staggered grid for 4.5 Gyr model time, i.e. from 0 Gyr to the present day at 4.5 Gyr. The model domain is resolved by 512 × 96 cells with radial grid refinement up to a factor of two at the surface and the core-mantle boundary (CMB), giving ∼ 13 km grid spacing at the surface and ∼ 25 km at the CMB in the radial direction. The grid spacing in the azimuthal direction ranges from ∼ 38 km at the CMB to ∼ 74 km at the surface. About half a million tracers (about 10 tracers per cell) are used to track the non-diffusive advection of composition, temperature, and heat-producing elements (HPE) within the domain. Solid and molten tracers can be created or merged during partial melting and freezing such that arbitrarily small melt fractions can be generated or frozen even with very few initial tracers per cell. Tests indicate that increasing the number of tracers makes no qualitative difference to the results (see Supplementary Material).
Free-slip mechanical boundary conditions are applied both on the upper and lower boundaries of the domain. Thermally, a constant temperature of 740 K is applied at the top boundary, while the bottom boundary is isothermal and cools with time due to removal of heat from the core. The parameterisation of core cooling is similar to that in Tackley (2004, 2010), which takes into account gravitational energy released by inner core growth and latent heat release. Internal heating is produced by HPE that initially have a uniform concentration over the whole domain but subsequently partition between solid and melt. The partition coefficient for HPE is set to be 0.001, which means that the concentration of HPE in the melt is 1000 times that in the coexisting solid. The initial temperature profile is adiabatic with a potential temperature of 1800 K, thermal boundary layers 50 km thick at top and bottom, and random temperature perturbations of peak amplitude 20 K throughout the domain. The initial CMB temperature of 4025 K is identical or similar to that used in previous models (Armann and Tackley, 2012;Gillmann and Tackley, 2014;Rolf et al., 2018;Uppalapati et al., 2020); tests with higher values of up to 6025 K indicate that the initial CMB temperature only influences the early part of the evolution and not the long-term dynamics (see Supplementary Material), consistent with previous models for Earth    (2012). R CMB is from Stevenson et al. (1983), C p from Gülcher et al. (2020); heatproducing isotopes of K, U and Th are lumped together into a single "HPE" component with D HPE being the order-of-magnitude partition coefficient (e.g. Xie and Tackley (2004)), t half an average half-life Nakagawa and Tackley (2005) and the heating rate H is based on Bulk Silicate Earth models (Palme and O'Neill, 2007).

Composition and Phase Transitions
The parameterisation of composition is based on mineralogy: material is divided into olivine (ol) and pyroxene-garnet (px-gt) phase systems. The composition varies between two end members, basalt and harzburgite, which are linear mixtures of the ol and px-gt mineral systems. Basalt is assumed to be 100% pyroxene-garnet, and harzburgite is assumed to be 75% olivine and 25% pyroxene-garnet. An initial composition of pyrolite (20% basalt and 80% harzburgite) is assumed for the entire mantle. Both the ol and px-gt systems undergo pressure-and temperature-dependent solid-solid phase transitions.
The depth, reference temperatures, density jumps, and Clapeyron slopes for the phase transitions (see Appendix, Table A

Rheology
A visco-plastic rheology is assumed (see Appendix, Table A.2). Viscous deformation includes both diffusion creep and dislocation creep, with the composite viscosity (η viscous ) given by The diffusion creep viscosity (η diffusion ) is governed by an Arrhenius-type law.
where T and P are temperature and pressure, the subscript i represents different phases, η 0 is the reference viscosity at T 0 = 1600 K and zero pressure, ∆η i is a phase-dependent viscosity offset (calculated at phase change T and P to give the relative viscosities listed in Table A.2), E a and V are the activation energy and activation volume for diffusion creep, and R = 8.314 J mol −1 K −1 is the gas constant. The dislocation-creep viscosity, on the other hand, is calculated as: where E a,i , V i , and n i correspond to the activation energies and activation volumes for dislocation creep, and stress exponents for different phases. σ is the second invariant of the deviatoric stress tensor and σ 0 is the crossover stress between diffusion creep and dislocation creep: σ 0 = 0.3 MPa for upper mantle phases (based on the below-cited laboratory measurements and for a grain size of 0.825 mm), and σ 0 = 3 MPa for lower mantle phases (slightly higher than typical convective stresses based on the common view that Earth's lower mantle is dominated by diffusion creep, e.g. Ammann et al. (2010)) and consistent with experimental results (Tsujino et al., 2022).
The activation volumes for perovskite phases are pressure-dependent (Tackley et al., 2013) and follow: For the shallowest phase in the px-gt system, which represents basalt before the basalt-eclogite transition, a plagioclase (An 75 ) rheology from Ranalli (1995) is applied, as in Gülcher et al. (2020). For other phases in the ol and px-gt systems, the rheological parameters are similar to the parameters used in Lourenço et al. (2016), which are based on data for olivine from Karato and Wu (1993) for the upper mantle and bridgmanite (Yamazaki and Karato (2001); Ammann et al. (2009)) for the lower mantle. Two viscosity jumps are applied in our models (Table A.2): a viscosity jump of 3.0 for the transition from upper-mantle phases to lower-mantle phases (similar to that in current Earth viscosity models (Čížková et al., 2012) and compatible with geoid and admittance ratios on Venus (Rolf et al., 2018), and a factor of 0.0133 for the transition from eclogite to plagioclase when the plagioclase rheology is applied to crust (calculated from the above-cited laboratory measurements).
These viscosity jumps are applied at the relevant phase transition (T,P) and result in different ∆η i in Equation (2) applicable at the reference (T,P). A reference viscosity of 1.0 × 10 20 Pa · s is used for all our models. In a grid cell that contains different phases, the overall viscosity of the cell is calculated as the geometric mean of the viscosities of the individual phases; the geometric mean is used because it is in-between the limiting cases (Yamazaki and Karato, 2001) of arithmetic mean (strong phase dominates by forming a load-bearing framework) and harmonic mean (weak phase dominates by forming interconnected weak layers).
Plastic yielding is also included in some of our models and can act to weaken and mobilise the lithosphere (Moresi and Solomatov, 1998;Tackley, 2000). The yield stress in our models follows a Drucker-Prager yield criterion: where φ is the friction coefficient and C is the cohesion. The effective viscosity η effective is calculated as: where˙ is the second invariant of the strain rate tensor. Finally, the viscosity is truncated between the lower and upper limits of 10 18 and 10 26 Pa · s, respectively.

Melting
The parameterisation of melting and crust production in the model is extended from that in Lourenço et al. (2020). Depth-dependent solidus and liquidus curves are used to calculate melt production, with the solidus applying to pyrolite (20% basalt) and increasing linearly as basalt is depleted. The solidus is fitted to Herzberg et al. (2000) for the upper mantle and Zerr et al. (1998) for the lower mantle. Melt that forms at depths shallower than the depth of neutral density (∼330 km on Venus) is either erupted or intruded.
A fraction of the produced melt given by the "eruption efficiency" input parameter is placed on the surface and the rest is intruded below the crust, the In reality, the viscosity of molten rocks is much lower than corresponding solid rocks, for example 10 2 − 10 13 Pa · s for hydrous leucogranitic melts (Hess and Dingwell, 1996). However, due to computational limitations, the viscosity of melts in our models is artificially set to be the minimum viscosity in the domain (10 18 Pa · s). Melts that are not erupted or intruded, like those deeper than 330 km, are advected with the flow and may later freeze.

Diagnostics
In this study, we computed the surface age and crustal thickness to compare our results with inferences from satellite observations and previous nu- We also calculate surface mobility to characterise the tectonic activity at the surface. The surface mobility M is defined as the ratio of the root mean square (rms) of surface velocity to the rms velocity averaged over the entire domain, in the no-net-lithospheric-rotation frame of reference (Tackley, 2000): For stagnant-lid tectonics, M is near zero, due to negligible surface velocities. For mobile-lid and plutonic-squishy-lid tectonics, M is larger than zero, reflecting either a moving lithosphere mostly pulled by subduction (Nakagawa and Tackley, 2015) or lithospheric movements facilitated by intrusions (Lourenço et al., 2020). Episodic-lid tectonics, on the other hand, is characterised by occasional lithospheric overturns with high surface velocities (mobile-lid phase, M > 1), and low surface velocities between overturns (stagnant-lid phase, M ≈ 0). Many numerical studies focussing on the catastrophic resurfacing model suggest that Venus is in the episodic-lid regime (e.g. Armann and Tackley, 2012;Rolf et al., 2018;Uppalapati et al., 2020) or transitioning from mobile-lid to stagnant-lid (Weller and Kiefer, 2020), whereas Byrne et al. (2021) suggests that Venus could be in a plutonic-squishy-lid regime (Lourenço et al., 2020) from satellite images of lowland plains. For a review see Rolf et al. (2022).

Parameter Study
In this study, the two main properties that we tested are crustal rheology and eruption efficiency. As discussed in the introduction, previous studies (e.g. Gülcher et al., 2020;Jimenez-Diaz et al., 2015) suggest that a relatively weak plagioclase (An 75 ) rheology (Ranalli, 1995) may be appropriate for modelling Venus' crust. Therefore, in terms of rheology, we tested (1) the influence of this "weak" crustal rheology compared to the previously-used olivine-based "strong" rheology (e.g. Armann and Tackley In addition, different eruption efficiencies (E) are tested in our models.
Eruption efficiency is defined as the percentage of melt that is erupted to the surface, with the remaining melt being intruded below the crust. Previous geodynamic models that included melting often assumed that all melt erupts to the surface (i.e., E = 100%). However, for Earth's magmatism, the estimate of eruption efficiency is only 10 -20% on average, depending on tectonic setting (Crisp, 1984). The highest estimate for eruption efficiency in Crisp (1984) is about 60%. Also, values of 20 % extrusion efficiency have been found to better match observations for the early Earth (e.g. Rozel et al., 2017), which is thought to be an analogue for modern-day Venus. Therefore, pure eruptive magmatism would not be a realistic approximation. Recent studies (Lourenço et al., 2018(Lourenço et al., , 2020 found that a low eruption efficiency (high intrusion efficiency) could lead to a new tectonic regime, the plutonicsquishy-lid tectonic regime, which could explain several Venus' observations (Byrne et al., 2021;Smrekar et al., 2022). Therefore, we vary the eruption efficiency from 20% to 100% to test how it influences volcanic and tectonic activity in our models.

Results
Here we present 14 numerical models with different rheological parameters and eruption efficiencies (summarised in Table 2). We first present the reference model, E20P03Pl, which has a "weak" plagioclase crustal rheology, dislocation creep, standard yielding parameters similar to Gülcher et al.
(2020), and 20% eruption efficiency (80% intrusion). We then examine the influence of differences in rheology and the effect of eruption vs. intrusion of melt on the surface and mantle dynamics of Venus.

Behaviour of the Reference Model
The evolution of the reference model E20P03Pl can be seen in Figure 2 and is characterised by global overturns, but with continuous high surface mobility in between and high surface strain rates through 4.5 Gyr ( Figure 3).  . For some models (marked with asterisks ( * )), the number of global overturns is hard to count as there could be a series of regional overturns that reset the average surface age to zero in a relatively short period of time. The results for surface age and crustal thickness include the average values and standard deviations at the end of simulation (4.5 Gyr), apart from model E20P50Pl (which is before the last global overturn).
the positive buoyancy of basalt in the depth range between the olivine-system phase transition to perovskite and the pyroxene-garnet system phase transition to perovskite, as previously found in other simulations (Ogawa, 2003;Tackley et al., 2005;Nakagawa and Buffett, 2005; and named the "basalt barrier" (Davies, 2008;Papuc and Davies, 2012). Between overturns, both surface mobility and surface strain rate ( Figure   3) indicate continuous surface deformation, which is contrary to the near-zero surface mobility between overturns in previous episodic-lid models. The average surface strain rate between overturns in model E20P03Pl (∼ 10 −16 s −1 ) is higher than the strain rate estimated from impact craters (10 −18 − 10 −17 s −1 ) on Venus' surface (Grimm, 1994). At 4.5 Gyr, we observe small old blocks separated by regions with low surface age and high strain rates (Figure 3 (a) and (c)), which could agree with the fragmented lowland basins observed on Venus' surface (Byrne et al., 2021).
Magmatism also occurs between overturns, as indicated by regions with high strain rates and near-zero surface age in Figure 3 (a) and (c). How-   and (d) are for E20P03Ol, which is the same but with a strong olivine crustal rheology.
When a global lithospheric overturn occurs, the temperature difference between upper and lower mantles is effectively homogenized. The basalt fraction above the CMB also significantly decreases during global overturns. The basaltic layer at the base of the MTZ forms within the first 0.5 Gyr and exists for the whole 4.5 Gyr model time. ever, as these features do not persist over time, there are no persistent mantle upwellings in model E20P03Pl, but rather randomly-located, short-lived magmatism.
The evolution of surface age and crustal thickness of model E20P03Pl is shown in Figure 3 (b). The average surface age decreases to zero when global lithospheric overturn occurs, but the average crustal thickness remains nonzero for some global overturn events, such as the one at 2.1 Gyr. This is because these global overturns occur as a series of regional overturn events that cover parts of the surface over several million years. The increase in crustal thickness after overturns indicates an increase in melting and crustal production during and immediately after these overturns. Between global overturns, there are regional magmatic resurfacings shown by low-age regions in Figure 3 (c). As radiogenic heat production in the mantle and CMB heat flow decrease over time, the interval between global overturns and the occurrence of regional magmatism decreases. As a result, the surface age distribution of model E20P03Pl at 4.5 Gyr is bimodal, with ∼ 20% of the surface younger than 100 Myr, and the rest being around 400 Myr old, corresponding to the last global overturn at ∼ 3.9 Gyr ( Figure 6).

Influence of rheology
3.2.1. "Weak" vs "Strong" Crust To investigate the impact of crustal rheology, we compare the model E20P03Ol (with "strong" olivine crustal rheology) to the reference model E20P03Pl (with "weak" plagioclase crustal rheology). The differences can be seen in Figure 4 and by comparing Figures 3 and 7. As in the reference case, the "strong" crust model E20P03Ol is also characterised by global overturns and a basalt barrier at the base of the MTZ; however, the basalt fraction at the base of the MTZ is lower when compared to the reference case E20P03Pl, while crustal thicknesses are larger. In addition to the global lithospheric overturns observed in the reference model E20P03Pl, regional resurfacing events are observed in E20P03Ol (Figure 8). Magmatism associated with these regional resurfacing events is generally weaker and restricted to hemispheric or sub-hemispheric extent. Such regional resurfacing events are indicated by small peaks in the average crustal thickness, or periods with high surface mobility but small changes in mantle temperature (e.g. ∼ 1.9 Gyr and ∼ 3.0 Gyr in Figure 7 (b)). The existence of such regional resurfacing events is in contrast to some previous episodic-lid models in which lithospheric overturn occurs over the whole surface (e.g. Uppalapati et al., 2020). In our simulations, plastic yielding still propagates through the entire lithosphere (Figure 8), however, at the boundaries of mantle plumes, the viscosity of the lithosphere and the underlying mantle asthenosphere is further reduced by dislocation creep due to high strain rates in this region.
As a result, crustal production by mantle plumes is balanced by crustal delamination, rather than destabilising the lithosphere and triggering a global overturn. Regarding surface mobility, unlike the continuous surface mobility observed in reference case E20P03Pl, the surface mobility in the "strong" crust case E20P03Ol is similar to that in previously-published episodic-lid models (e.g. Armann and Tackley, 2012): high during global overturns and nearzero between overturns (Figure 7 (d)). The interval between the overturns is also longer than that of E20P03Pl, possibly due to the higher viscosity of the crust. Due to the long intervals between overturns, the average surface age in case E20P03Ol is higher than that in reference case E20P03Pl. The average crustal thickness (60 -75 km between overturns) is also higher than that in E20P03Pl, and higher than the current estimated crustal thickness for Venus.
In terms of magmatism and surface deformation, the mean surface strain rate between overturns decreases to ∼ 10 −18 s −1 , or two orders of magnitude lower than that of E20P03Pl. Unlike E20P03Pl, there are regions with longlived relatively high strain rates and low surface ages in E20P03Ol ( Figure 7).
As these regions indicate magmatic crust production, the mantle upwellings that produced these melts could be stationary for hundreds of millions of years in E20P03Ol. Such persistent mantle upwellings lead to non-uniform resurfacing rates at the surface, which are also observed in episodic-lid models from other studies (Uppalapati et al., 2020).

Dislocation Creep
As shown by the fractions of deformation accommodated by different rheological mechanisms shown in Figure    compared with E20P30Pl and E20P30Ol. Surface mobility depends mostly on crustal rheology ("strong" vs. "weak") rather than the inclusion of dislocation creep. Both models with "weak" plagioclase crustal rheology have continuous surface mobility, whereas both models with "strong" olivine crustal rheology have episodic-lid surface mobility (Figure 9 (d) and (e)). Regional resurfacing is also observed in E20P03Ol D at 0.8 Gyr (indicated by high surface mobility and crustal thickness but low surface age in Figure 9). The mobility peaks at 3.9 Gyr are caused by several small-scale mantle upwellings starting at 3.7 Gyr that cause crust delamination at 3.9 Gyr.

Plasticity
As in previous geodynamic models of Venus, plastic yielding is applied in our models as a weakening mechanism for the lithosphere. Plasticity represents brittle failure/fault slip, and in this study we use a friction coefficient that is fixed at 0.2 and three different values of cohesion: C = 0.3 MPa, C = 50 MPa, and infinity (no plastic yielding). This "brittle" yield stress formulation, as commonly used in the lithospheric dynamics community, is different from the formulation of "ductile" yield stress commonly used in the mantle convection community in that the latter is constant or increases only slowly with pressure (depth). Previous global numerical models often used both formulations (e.g. Armann and Tackley, 2012), but in this study only brittle yield stress is applied in our models, which is similar to the regional models of coronae formation (Gülcher et al., 2020). Nevertheless, both formulations lead to similar dynamical regimes in global models (Moresi and Solomatov, 1998;Tackley, 2000). How different "ductile" yield stresses or "brittle" friction coefficients affect the overturn frequency and crustal thickness has been discussed in previous numerical studies (e.g. Nakagawa and Tackley, 2015;Rolf et al., 2018;Uppalapati et al., 2020). Here, the two models with C = 50 MPa exhibit similar tectonics compared to models with C = 0.3 MPa. Therefore, we focus on the two models without plastic yielding E20P00Pl and E20P00Ol, which differ in having a "weak" vs. a "strong" crust.
The model without plasticity and a "strong" crust (E20P00Ol) is the only model presented in this study that exhibits stagnant-lid tectonics, indicated by near-zero surface mobility for the whole 4.5 Gyr model time (Figure 11 (d)). This shows that the weakening from dislocation creep and intrusive magmatism is less effective than plastic yielding when the crustal creep rheology is "strong" (olivine-like).
Interestingly, this model exhibits three global resurfacing events, as indicated by zero surface age and a sudden increase in crustal thickness ( Figure   11 (b), e.g. at ∼ 1.8 Gyr). However, these are not due to lithospheric overturn -they are instead due to internal overturn of the mantle ("avalanches") modulated by the "basalt barrier" mechanism. Hot lower-mantle material enters the upper mantle and causes a global pulse of magmatism. Heat transfer through the lithosphere is dominated by the magmatic heat-pipe mechanism (O'Reilly and Davies, 1981;Solomon and Head, 1982;Armann and Tackley, 2012) and new eruptions to the surface lead to zero surface age and an increase in crustal thicknesses. Between these global internal mantle overturns, there are also regional magmatic resurfacing events ( Figure   11). As a result, although there is no lithospheric overturn, the evolution of surface age in this stagnant-lid model is similar to in episodic-lid models, and reaches around 550 Myr at present-day. On the other hand, the average crustal thickness in E20P00Ol remains nonzero throughout the model time.
After mantle overturns, the crustal thickness first increases to larger than 100 km, as global mantle overturn triggers global magma production beneath the existing crust ( Figure 12). Then, the crustal thickness starts to decrease and reaches ∼ 75km within several hundred Myr due to crustal delamination, which is similar to in other models with olivine crustal rheology and possibly limited by the basalt to eclogite phase transition (reference depth is 66.4 km, see Table A.1). With a "weak" plagioclase crustal rheology (model E20P00Pl), global lithospheric overturns are observed despite the lack of plastic yielding ( Figure   13), indicating that the combination of a plagioclase crustal rheology, dislocation creep and intrusive magmatism is sufficient to weaken the lithosphere and facilitate overturns. The evolution of this model is generally similar to other models with "weak" crustal rheology, with global overturns that reduce the average surface age to near zero, and regional resurfacing between global overturns. The number of global overturns is less than in models with plastic yielding (see Table 2), and the surface mobility is higher than zero through

Influence of Eruption Efficiency
In order to systematically investigate the effect of varying eruption efficiencies on the tectonic regime, models with 20%, 40%, 80%, and 100% eruption efficiencies are tested in this study. For models with a plagioclase crustal rheology, the surface mobility is continuously high. On the other hand, for models with olivine crustal rheology and > 80% eruption efficiency, the intervals between overturns are too short for stagnant-lid phases. Therefore, the surface mobility is also continuous over the 4.5 Gyr model time.
In models with a strong (olivine rheology) crust, higher eruption efficiency leads to more frequent overturns. These models also have lower ratios between the conductive and magmatic heat flow, or R C/M (Figure 14). As the heat-pipe mechanism is less efficient at cooling down the mantle over time than intrusive magmatism (Lourenço et al., 2018), the heat transfer efficiencies in mantle are lower with higher E, resulting in more frequent overturns to release heat. On the other hand, for models with a "weak" plagioclaserheology crust, the numbers of global overturns are similar between models with different E, indicating that crustal rheology might be more important than heat transfer efficiency in terms of overturn dynamics for these models.
The average surface ages are lower with higher eruption efficiency, due to extensive partially molten zones ( Figure 15). The average R C/M is similar between different models, and is larger than 1 for most of the model time ( Figure 14). When comparing models with the same eruption efficiency and different crustal rheology, R C/M is generally larger for models with plagioclase rheology, suggesting that heat transfer by intrusive magmatism, which facilitates conductive heat transfer by forming a thin lithosphere (Lourenço et al., 2018), is enhanced by a weaker crustal rheology.
Due to more frequent overturns and extensive volcanism, models with higher eruption efficiencies have a lower average surface age ( Figure 16).
Unlike previous numerical studies, in which increasing E leads to a thicker  crust (Lourenço et al., 2018), the crust is thinner with higher E in our models.
Also, there are partially molten zones beneath the crust. With higher E, the upper mantle is hotter and more depleted (Figure 17), resulting in thinner crust at the end of the model.

Tectonics and Volcanism
Observations of Venus' surface suggest that its resurfacing history could be a combination of two end-member models: global or near-global lithospheric overturns reworking most of the planet's surface and leading to a unimodal crater distribution, with small-scale tectonic and magmatic resurfacing events that created some of the highly deformed regions and young  is near zero. However, in our current models with composite crustal rheology, dislocation creep, and intrusive magmatism, we observe both global overturns and regional resurfacing, continuous surface mobility, as well as short-term magmatic events. This study thus illustrates the importance of applying a composite crustal rheology and intrusive magmatism to numerical models that focus on the tectonics of Venus.
The introduction of more realistic rheologies and intrusive magmatism changes two main aspects of the tectonic regime: (1) continuous surface mobility, and (2) the occurrence of regional resurfacing. The surface mobility is affected mainly by crustal rheology: models with "weak" plagioclase crustal rheology exhibit continuous surface mobility in-between overturns, while models with "strong" olivine rheology and low eruption efficiency have episodic-lid surface mobility. Regarding regional resurfacing, plastic yielding still controls the viscosity in the uppermost lithosphere, which can lead to near-global overturns that are similar to those in previous numerical models (e.g. Crameri and Tackley, 2016). However, as dislocation creep dominates in the upper mantle, the high strain rates around the mantle plumes and lower crustal eclogite (due to its high density) could facilitate the formation of small-scale convective cells, which can break the lithosphere above them and cause regional resurfacing that covers less than a hemisphere when small upwellings occur. These regional resurfacing events could explain how the local tectonic features on Venus survived lithospheric overturns, one example being the "drifting" of Lakshmi Planum (Harris and Bédard, 2015).
Models with different brittle cohesion values (C = 0.3 MPa and C = 50 MPa) display similar tectonic evolutions: the number of overturns and crustal thicknesses are similar between these models (Table 2) Regarding volcanism, another new observation in these models is that crustal rheology can affect the duration of local magmatic resurfacing. With olivine crustal rheology, upwellings in the upper mantle are stationary between global overturns and can be active for tens of millions of years. Resurfacing from magmatic crustal production is focused in the regions above these plumes, and the remaining surface lacks tectonic and volcanic activity, resulting in a non-uniform resurfacing rate. However, with a "weak" plagioclase crustal rheology, the time span of magmatic resurfacing is much shorter, possibly due to surface movements. Regions with higher-than-average strain rates are active for only a few million years. There are numerous shortlived, randomly located volcanic events between overturns, resulting in a more uniform resurfacing rate for Venus' surface, which is consistent with the cratering observations. The tectonics just described are also similar to the plutonic-squishy lid regime (Lourenço et al., 2020), but with some differences as discussed later.
Another observation is that in models with olivine rheology, there are more regional events than in models with plagioclase rheology, which could be counterintuitive. This can be explained by the fact that a plagioclaserheology crust is (relatively) weak everywhere, leading to lithospheric blocks that are not well defined because, due to their relatively fast movement, weak zones caused by intrusions cannot survive for a long time. Therefore, a plagioclase-rheology crust leads to a lithosphere that can be considered as a mobile lid for most of the model time. Regarding overturn events, multiple regional resurfacing events tend to occur within a short period and cover the entire surface (similar to but less vigorous than a global overturn) and reset the age of the entire surface. In contrast, in models with an olivine-rheology crust, as the surface is mostly immobile during the near-stagnant-lid phase between overturns, the weakening from intrusions forms a compositionally continuous but rheologically discontinuous lid, which can lead to regional overturns.

Crustal Thickness, Surface Age, and Topography
Estimates of the crustal thickness of Venus vary between 20 -60 km, depending on the geological province (Breuer and Moore, 2015;Wieczorek, 2015). As the crustal thickness is estimated by inverting gravity and topogra-phy observations, the result is highly dependent on the assumptions made, as the analysis by Orth and Solomatov (2012) shows. ing with 100% eruption is too simplified for magmatic crustal production for Venus because it often results in very thick (>100 km) crust, especially in stagnant-lid cases as crustal production is enhanced by higher mantle temperature (e.g. Armann and Tackley, 2012). More recent modelling studies, in which only basalt can melt, applied different eruption fractions in Earth-like models, and obtained thinner crusts with lower eruption fractions (Lourenço et al., 2018(Lourenço et al., , 2020. However, if Venus' interior is hotter than Earth's, harzburgite could also melt, as its solidus is 150-200K higher than the solidus of lherzolite (Maaløe, 2004), resulting in komatiitic-composition magma. In this study, a partially molten zone below the crust is often ob-  et al., 2021). However, in other deformation models, the formation of felsic tessera requires a much steeper geotherm (> 40 K / km) so that deformed craters are not preserved when folds are generated (Resor et al., 2021). A recent study also supports a thin lithosphere and high surface heat flow for Venus based on the topographic flexure at 65 coronae .
Such a steep geotherm could cause major melting below the base of the lithosphere (Ghail, 2015). Recent analysis of Venusian crustal plateaus also indicates that the melting of deep crustal material could occur when crustal plateaus formed (Maia and Wieczorek, 2022).
Previous numerical studies suggest that increasing eruption efficiency results in a thicker crust (Lourenço et al., 2018). In the current simulations, however, models with high eruption efficiencies display thin crusts and partially molten zones beneath the crust. The reasons are: (1) the temperature in the upper mantle is higher for models with higher eruption efficiency, causing more magmatic crustal production ( Figure 17); then (2) mantle overturns bring the thick basaltic crust to the CMB, where much of it remains. Therefore, in models with higher eruption efficiency, the basaltic crust is quickly produced and removed from the upper mantle, resulting in a more depleted upper mantle. Thus, at the end of these models, the melting process can only create a thin crust, with large, partially molten zones beneath the crust containing harzburgitic melt that cannot erupt ( Figure 15). For models with lower eruption efficiencies and plagioclase crust rheology, as crust dripping is easier with a "weak" crust, the maximum crustal thickness is thinner and better fits estimates (e.g. 36.7 km for E20P03Pl). As for the average surface age, in models with a weak plagioclase-rheology crust, the average surface age is mainly controlled by the time of last overturn. Although a wide range of surface ages between 150 − 1000 Myr has been suggested for Venus(e.g. Hauck et al., 1998;Herrick and Rumpf, 2011), only models with a low eruption efficiency (E < 40%) can fit the estimates. For models with high eruption efficiency, the average surface age is less than 100 Myr due to rapid magmatic crustal production.
As for topography, it is not straightforward to compare the model re-

Interior Evolution
It is now well established that including compositional variations and realistic composition-dependent density changes in numerical models has a first-order effect on the model behaviour. This can be seen in the simulations presented here, previous simulations for Venus (Armann and Tackley, 2012) and Earth (e.g. Nakagawa and Tackley  include such density variations using other numerical models (e.g. Papuc and Davies, 2012;Ogawa and Yanagisawa, 2014;Vesterholt et al., 2021). Compared to purely thermal models (e.g. Solomatov and Moresi, 1996;Huang et al., 2013), applying this composition-dependent density leads to significantly different mantle structure and dynamics, including:  (Bissig et al., 2022).
This "basalt barrier" effect (Davies, 2008) leads to compositional and thermal stratification between the upper and lower mantle, which can only be penetrated by strong mantle plumes. Similar dynamics is observed in models of Venus by Ogawa and Yanagisawa (2014), with mantle bursts that destroy the basalt barrier, strongly stir the mantle, and result in widespread magmatism at the surface. The episodic variation of upper mantle temperature and spikes in basaltic crustal thickness that are observed in our simulations, were also observed in Papuc and Davies (2012), although none of their models had a proper stagnant lid because they applied a viscosity cutoff at 100 times the reference viscosity. The strength of the basalt barrier is controlled by the ratio between compositional density differences and thermal density differences in the MTZ. Thus, this effect is expected to be stronger in a stagnant-lid mode than in a mobile-lid mode because compositional density differences (between basalt and harzburgite or pyrolite) are independent of tectonic mode, whereas temperature variations (and therefore thermal density differences) inside the mantle, are lower in stagnant-lid mode because the lithosphere does not participate in convection. Indeed, the stagnant lid case (E20P00Ol) does exhibit stronger layering, with only three internal overturns, than the episodic-lid cases.
2. Basalt settling above the core-mantle boundary. Such a basaltic layer is also observed in Earth models and its structure is strongly dependent on the tectonic mode (Nakagawa and Tackley, 2015;Lourenço et al., 2020). The present simulations obtain a continuous layer of basalt around the CMB, which is consistent with previous episodic-lid models; more Earth-like continuous subduction can result in isolated piles instead. The basalt layer accumulated above the CMB strongly reduces the CMB heat flux because this layer is enriched in HPE Tackley, 2005, 2014).
Global overturns cause rapid thinning of this layer thickness and pulses of heat from the core. Therefore, it is possible that Venus' magnetic field could be episodic over time and be recorded by ferromagnetic minerals in its crust (O'Rourke et al., 2019). The amount of basalt settling above the CMB is strongly sensitive to its density contrast with the ambient mantle, which is moderately uncertain; for example, the relevant density contrasts depend strongly on the assumed basalt composition (e.g.  and are smaller when calculated using a recent mineralogical database (Stixrude and Lithgow-Bertelloni, 2022) than when calculated using an older version of the database (Xu et al., 2008). Rolf et al. (2018) and Uppalapati et al. (2020) included magmatic crustal production in their models; however, the above compositional effects are not observed in their simulations. This is because their models did not include realistic composition-dependent density contrasts, as the density differences are only about 10 kg/m 3 between basalt and harzburgite in most of the mantle. Also, the density changes around the upper mantle -lower mantle transition zone are greatly reduced in their studies. Thus, the compositional density difference is too low to support a basalt barrier and form a separate basalt layer above the CMB. When applying a more realistic density profile to numerical models, for example by calculating phase assemblages and densities by using free energy minimisation (e.g. Nakagawa et al., 2009Vesterholt et al., 2021), the aforementioned strong compositional effects are observed.

New "deformable episodic lid" tectonic regime
Models with plagioclase crustal rheology and low eruption efficiencies show different behaviors compared to previously proposed stagnant-, episodic-, or plutonic-squishy-lid regimes, leading to a new regime that we name "deformable episodic lid". This regime is characterised by a generally high surface mobility plus episodic lithospheric overturns (in one model (E20P00Pl) even without plastic yielding). Between overturns this tectonic regime shares some similarities with a plutonic-squishy lid (Lourenço et al., 2020) regime: downwellings are dominated by crustal delamination and/or drips rather than subduction. However, these two regimes are different in that in a plutonicsquishy lid the lithosphere is split into small rigid plates or blocks separated by weak boundaries due to intrusions, while in the "deformable episodic lid" lithospheric deformation is not localized into "plate boundaries" but instead high almost everywhere, making the lithosphere rheologically continuous for most of the time. Models with olivine crustal rheology are more similar to the plutonic-squishy lid regime (Lourenço et al., 2020) and show more long-lived volcanism when compared to models with a plagioclase rheology crust. These differences arise because the primary weakening mechanism in a plutonic squishy lid is magma intrusion, while the primary weakening mechanism in a deformable episodic lid is a relatively weak and non-linear crustal rheology.

Global volcanic resurfacing events
Here we find that global volcanic resurfacing events may occur in a stagnant-lid regime due to the "basalt barrier" effect at the MTZ ( Figure   11) discussed above. This volcano-tectonic regime can be described as stagnant episodic-global-volcanic-resurfacing. Similar resurfacing behavior is observed in Ogawa and Yanagisawa (2014) and Papuc and Davies (2012) and Yuen, 1985). The viscosity contrast between upper and lower mantles chosen for our models here is consistent with modern viscosity models of Earth's mantle (e.g. Čížková et al. (2012) and is also more consistent with the geoid and the correlation between geoid and topography of Venus (e.g. Rolf et al., 2018).

Limitations of this study
Several parameters and assumptions in this study are subject to significant uncertainties and discussion. Regarding rheology, the diffusion-creep rheology in this study assumes a uniform grain size giving η 0 = 1×10 20 Pa · s for all models. However, a wide range of grain sizes is expected in both Earth's and Venus' mantles and grain size evolution can have important effects on the dynamics (e.g. Rozel, 2012;Dannberg et al., 2017;Schierjott et al., 2020). The evolution of the grain sizes for mantle minerals could also affect the development of lithospheric shear zones through strain weakening and shear localization (e.g. Bercovici and Ricard, 2014). Additionally, the influence of the viscosity contrast between upper and lower mantles (here assumed to be small for diffusion creep and large for dislocation creep, based on a combination of Earth-based constraints and Venus' admittance ratios) is worth exploring in the future.
The parameterisation of melting in our models is based on melting processes on Earth. Venus' mantle solidus and HPE concentration could be different due to compositional variations (e.g. Kiefer et al., 2015;Shah et al., 2022). The melting and eruption behaviour also needs further analysis, as it strongly affects magmatic crustal production.
Finally, no three-dimensional (3-D) model is included in this study due to their computational expense. In 2-D, mantle plumes are infinite sheetlike features, which may exaggerate their resurfacing effects. On the other hand, models in Uppalapati et al. (2020) comparing 2-D and 3-D geometry show similar crustal thickness and number of major overturns, indicating similar tectonics between models with different model geometries, as did the Earth models of . Internal overturns modulated by the "basalt barrier" mechanism might also manifest differently in 3-D than in 2-D, perhaps being more regional as are "avalanches" driven by the ringwoodite to bridgmanite+ferropericlase transition (e.g. Tackley et al. (1993).

Conclusions
Using 2-D mantle convection models, we carry out a systematic modelling study of Venus' tectonics and interior evolution, considering the effects of a weak, experiment-based crustal rheology and intrusive magmatism. Unlike previous 2-D or 3-D models, where an episodic-lid or stagnant-lid tectonic mode were obtained, the results from this study show that Venus' tectonics can be a combination of global overturns, regional overturns, and shortterm, randomly-located magmatic resurfacing events over the entire surface.
Global resurfacing events can be caused either by lithospheric overturns or by internal overturns driven by breakdown of a "basalt barrier" between the upper and lower mantles, what leads to hotter lower-mantle material flooding the upper mantle and causing widespread volcanism.
The effect of applying a weak crustal rheology is both reflected in surface mobility and in upwelling activity. Models with a "weak" plagioclase crustal rheology exhibit both relatively high, continuous surface mobility and episodic lithospheric overturns, a new mode that we term deformable episodic lid, while models with a "strong" olivine crustal rheology display standard episodic or stagnant-lid mobility. For models with with a "strong" olivine crustal rheology, surface motion between global overturn events can sometimes occur due to regional overturns, rather than having a completely rigid lid. The duration of regional magmatic resurfacing is also affected by crustal rheology. For models with an olivine rheology, upper mantle upwellings are near-stationary between overturns and can last for tens of millions of years, resulting in non-uniform magmatic resurfacing rates at the surface. However, for models with a "weak" plagioclase rheology, magmatic resurfacing occurs more randomly and exists only for a short time, leading to a more homogeneous resurfacing rates in these models.
Regarding lithospheric weakening mechanisms, a combination of weak crustal rheology, including dislocation creep and intrusive magmatism, can weaken the lithosphere sufficiently for global lithospheric overturns to occur without plastic yielding. When included, plastic yielding still dominates in the top-most lithosphere, but dislocation creep can further weaken the lithosphere around mantle plumes, resulting in regional overturns.
The interior dynamics is strongly affected by the composition-and phasedependent density and dislocation creep. Composition/phase-dependent density leads to a strong "basalt barrier" at the base of the MTZ -stronger than occurs in Earth-like models, creating thermal and compositional stratification between the upper and lower mantle, plus a layer of basalt above the CMB. The formation of this basalt barrier can also cause internal global mantle overturns in stagnant-lid models, resulting in global magmatic resurfacing events without lithospheric overturn. Tectonics is strongly influenced by dislocation creep, as this creep mechanism dominates in the upper mantle.
Dislocation creep lowers viscosity, resulting in more efficient heat transfer, lower temperatures, and less magmatic resurfacing compared to identical models without dislocation creep.
Because of the different interpretations of tectonic features on Venus' surface, it is unclear which tectonic mode Venus experiences, but in general our models in a deformable episodic lid regime better match current observations and estimates, including crustal thickness, surface age, and the presence of both global resurfacing events and regional tectonic and volcanic activity between these events. Also, the episodic-volcanic resurfacing mechanism described for stagnant-lid models suggests a new global resurfacing mechanism for Venus, where global episodes of eruption create new crust that covers the whole planet. This mechanism can match the same observations and estimates as episodic-lid models with olivine crustal rheology, but without the need for global lithospheric overturns.

Acknowledgements
JT was partially supported by the ERC project NONUNE (Grant agree-   Table A.2: Rheological parameters for olivine and pyroxene-garnet phase system that are used in our models. The basalt-eclogite phase transition is marked with asterisks ( * ). If the plagioclase rheology is applied, its rheological parameters will replace parameters for basalt facies in the pyroxene-garnet system. Depths of phase changes have been scaled from Earth according to Venus' slightly lower gravity. We performed tests for our reference case (E20P03Pl, initial CMB temperature = 4025 K) with different initial CMB temperature (T CMB = 4525 / 5025 / 5525 / 6025 K). While near the beginning, models with higher initial T CMB show more frequent overturns, after ∼ 2.5 Gyr the dynamics are similar to the reference model in terms of surface mobility, frequency of global overturns, occurrence of regional resurfacings, and average surface age and crustal thickness between global overturns. Therefore, the initial T CMB only affects the early evolution, not the long-term behaviour. The figures below may be compared to Figure 3 in the main paper.