Dynamics of fault motion and the origin of contrasting tectonic style between Earth and Venus

Plate tectonics is one mode of mantle convection that occurs when the surface layer (the lithosphere) is relatively weak. When plate tectonics operates on a terrestrial planet, substantial exchange of materials occurs between planetary interior and its surface. This is likely a key in maintaining the habitable environment on a planet. Therefore it is essential to understand under which conditions plate tectonics operates on a terrestrial planet. One of the puzzling observations in this context is the fact that plate tectonics occurs on Earth but not on Venus despite their similar size and composition. Factors such as the difference in water content or in grain-size have been invoked, but these models cannot easily explain the contrasting tectonic styles between Earth and Venus. We propose that strong dynamic weakening in friction is a key factor. Fast unstable fault motion is found in cool Earth, while slow and stable fault motion characterizes hot Venus, leading to substantial dynamic weakening on Earth but not on Venus. Consequently, the tectonic plates are weak on Earth allowing for their subduction, while the strong plates on Venus promote the stagnant lid regime of mantle convection.


Strength of rocks
The strength of the lithosphere is controlled either by the stress needed to deform a rock plastically or by the stress needed to fracture a rock. In the deep lithosphere where temperature and pressure are high, the strength is controlled by plastic flow. In this regime, the stress σ = σ 1 − σ 3 ( ) (σ 1 : the maximum compressive stress, σ 3 : the minimum compressive stress) needed to deform the lithosphere depends strongly on various variables as 1,2 , where T is temperature, P is pressure, ! ε is strain-rate, d is grain size and C W is water content. The functional form, F T ,P, ! ε,d,C W ( ) , depends on the mechanisms of plastic deformation. Three mechanism of deformation are considered 2 and the source of the data are summarized in Table S1: diffusion creep, power-law dislocation creep and the low-temperature plasticity (the Peierls mechanism). To a good approximation, plastic rheology of a rock is isotropic and the strength given by equation (S-1) depends only on the magnitude of deviatoric stress (σ = σ 1 − σ 3 and does not depend on the stress state (compression or tension). The strength depends on grain-size in diffusion creep regime, but the strength is independent of grain-size for other mechanisms.
Recently, Kumamoto et al. 3 suggested that the strength in the lowtemperature plasticity may depend on the size such as grain-size or the size of the indentation in such a way that the strength corresponding to low-temperature plasticity for a large geological grain-size might be lower than previous model.
However, we consider that the experimental basis for this hypothesis is not strong.
For example, the data point by Druiventak et al. 4 on a coarse grained sample plays a key role in order to justify this model. However, the experimental study by Druiventak et al. 4 was conducted in the semi-brittle regime, and consequently, their strength likely under-estimates the true strength corresponding to low-temperature plasticity.
In addition, the data by Evans and Goetze 5 on a single crystal was plotted for the "size" of ~3 µm assuming that in their experiments, the relevant "size" was the size of indentation. This is highly questionable. In these micro-indentation experiments, the size of indentation is measured as a function of load, and a broad range of indentation size is explored. The results show that the hardness (strength) is independent of the load (size).
Kumamoto et al. 3 also used the "agreement" between Druiventak et al. 4 and the results of dislocation velocity measurements to support their model. The calculation of low-temperature strength from the dislocation velocity measurements involves a number of assumptions and is not straightforward. One of the important factors is the influence of the use of a very thin sample. The stress acting on a dislocation in a very thin sample includes the image stress caused by the presence of a free surface 6 , and consequently, applications of such measurements to estimate the "true" dislocation mobility are complex and the uncertainties are large.
For these reasons, we will not include the "size effects" on low-temperature plasticity suggested by Kumamoto et al. 3 (see also 7 ).
The strength corresponding to diffusion creep depends strongly on grain-size.
We consider the grain-size of 1, 10, 100 microns (for a discussion on grain-size see grain size section). Strength in this regime is generally sensitive to water content, but since the oceanic lithosphere is depleted with water (see Water in the oceanic lithosphere section), we use the flow law parameters for dry conditions. The results summarized in Fig. 3 show that even an unusually small grain-size of 1 micron, the strength is still high in the shallow, low temperature regions.
In the shallow lithosphere, deformation occurs by brittle fracture. In the brittle regime, deformation is localized and occurs along the fault plane. A common assumption is that there are many pre-existing faults in the lithosphere, and the strength in this regime is controlled by the resistance against the motion of preexisting faults.
Experimental studies show that a fault moves when the shear stress on the fault plane exceeds a critical value given by 8,9 , where σ n is the normal stress to the fault, τ is shear stress (τ = σ 1 −σ 3 2 ), µ is friction coefficient, and P pore is the pore pressure (it is assumed that water fills the fault plane and is connected to the surface) and we ignored the cohesive strength.
The shear stress (τ ) on the given fault plane depends on the applied stress (σ 1 ,σ 3 ) and their orientations relative to the fault plane. Consequently, the resistance for the fault motion depends on the stress state: fault motion is more difficult for the compressional stress state (τ thrust : thrust fault) than for the tensional stress state (τ normal : normal fault) (τ transform : transform fault is in between). A simple model shows 10 , where P is the (lithostatic) pressure, and hence, τ thrust > τ transform > τ normal . The higher value of strength for compression (thrust fault) than those for tension (normal fault) implies that it is more difficult to initiate subduction than to maintain subduction.
Equation (S-2) does not specify the velocity at which sliding occurs. As far as the shear stress on a fault exceeds the value given by equation (S-2), sliding starts.
Sliding will start with the velocity imposed by the geological boundary conditions (i.e., slow initial velocity). In some cases, sliding is stable but in other cases sliding is unstable. For a stable slow sliding, sliding velocity is low and the friction coefficient has a nearly universal value, µ ≈ 0.6 11 . However, if sliding is unstable, then the sliding velocity becomes high and the friction coefficient will be reduced to less than µ ≈ 0.1 12 when sliding velocity exceeds a critical value (~ 1 m/sec).
If one uses a canonical value of static friction coefficient ( µ = 0.6), then the average strength of the lithosphere substantially exceeds the limit for plate tectonics (~100 MPa) although the magnitude of strength in the brittle regime depends on the stress state (τ thrust versus τ normal ). If there are mechanisms to reduce the friction coefficient to ~0.1, then the lithosphere would be weak enough to allow plate tectonics.
The most plausible process to cause grain-size reduction is dynamic recrystallization, i.e., deformation-induced refinement of grain-size 19,20 . In such a case, the grain-size is primarily controlled by the differential stress 19,21 (Fig. S1).
Typical tectonic stress on Earth's lithosphere is 10-100 MPa 22-24 that would lead to the grain-size of 20-500 microns (typical stress in the hot asthenosphere is 0.1-1 MPa that would lead to the grain-size of 1-10 mm). Grain-size of rocks in the lithosphere is determined mainly when rocks were deformed at hot asthenospheric conditions 25 . Consequently, grain-size in most of the lithosphere is several mm [24][25][26] .
The density contrast caused by phase transformations in the subducted slab also contributes to the stress, but its magnitude can be larger than the above estimate 22,28,29 .
Rocks with grain-size less than 10 micron are very rare, and found only in areas with intense collision at the later stage of deformation where stress is concentrated in a narrow zone 14,30 . Even with extremely small grain sizes (e.g., 1 micron), the strength is still large in the shallow regions where temperature is low ( Fig. 2). Therefore we conclude that it is difficult to reduce the grain-size to the degree that the oceanic lithosphere is weak enough to bend and subduct.
Finally, we note that grain-size model for rheological weakening and our model are

Scaling the thermal weakening behavior
Laboratory high-speed friction experiments show a marked drop in the resistance for slip when slip velocity exceeds a critical value 12,31 . Laboratory experiments also show that this reduction in resistance occurs only beyond a certain displacement.
These laboratory experiments cover a broad range of slip velocity ranging from 10 -8 to 10 m/s corresponding to the rate of slip for slow earthquakes to normal earthquakes and there is no need for extrapolation in terms of slip rate 32 . However, the normal stress in these studies is limited to ~40 MPa 33 . When we apply these results to the deep lithosphere, say ~30 km depth, the normal stress (confining pressure) goes to ~1000 MPa. Consequently, one needs to address the validity of applying these laboratory data to deep lithospheric conditions with particular emphasis on the influence of normal stress (confining pressure).
Three parameters play a key role in the thermal weakening: the threshold velocity (Vc) and the threshold slip distance (Dth) for thermal weakening, and the "steady-state (final)" friction coefficient ( µ ∞ ). Let us address the scaling issues one by one.

Threshold velocity (Vc)
Weakening (a substantial reduction in the friction coefficient) at high velocities is likely due to temperature increase. Indeed, evidence of melting is often observed in samples where a substantial reduction in friction coefficient is observed 31,34 . However, some other processes that occur at high temperatures may also be responsible for the reduction in the friction coefficient including some thermally activated chemical reactions that form materials reducing the friction coefficient (tribochemical reactions) 12 .
In the following, we use the results by 12 , and derive a scaling relationship for the critical distance D th assuming that weakening occurs when temperature is high enough. When the velocity of sliding is high, large work is done on the fault plane and substantial heating will occur. This issue can be analyzed by solving a problem of temperature increase caused by the mechanical work done on the fault plane. A simple one-dimensional model predicts the temperature rise as 12 , where ΔT ( ) c the temperature increases needed for melting ( ΔT is the critical temperature for weakening and T o is the initial temperature),τ e is the effective stress, V c is the threshold velocity for melting, ρ is density, C P is specific heat, κ is thermal diffusivity and D th is the threshold slip distance for thermal weakening. Using the relationsτ e ≈ 0.64τ o ≈ 0.4σ n and D th = a ⋅σ n Hence, The critical temperature for frictional weakening is unlikely to change much with pressure in the pressure range that we consider (P<1 GPa) 35 . Therefore we conclude that the threshold velocity for thermal weakening decreases with pressure (depth), i.e., thermal weakening becomes easier with pressure (depth). This is because the amount of work done by friction increases with pressure (normal stress).
In the above model, frictional heating is assumed to occur in an infinitesimally thin layer. This approximation is valid when thermal diffusion distance ( πκV c D th ) exceeds the thickness of the fault plane (h) where shear deformation occurs 36 . If the thickness of fault plane exceeds thermal diffusion thickness ( h > πκV c D th ), then equation (S-4) should be replaced with (see 36 In such a case, threshold velocity will be inversely proportional to the effective stress and hence depth. We conclude that at a deeper region where deformation is more diffuse (in the semi-brittle regime), critical velocity for thermal weakening is lower. In an extreme case where shear is very diffuse, then in a plausible effective stress, temperature rise is not enough to cause thermal weakening. This will be the case in the deep region (~40 km in the 60 Myrs old oceanic lithosphere (see Fig. 2)) where the dominant mode of deformation gradually changes to ductile flow (e.g., 37 ). In this study, we focus on relatively shallow lithosphere (~10-30 km), and therefore the influence of shear zone thickness will not be important.

Threshold slip distance for thermal weakening (Dth)
The thermal weakening occurs only after a finite slip. The displacement due to the fault motion (D) must far exceed the threshold slip distance, Dth (D>>Dth).
Again this is due to the fact that one needs to produce a certain amount of heat to melt the rock. At a typical laboratory condition where normal stress is ~10 MPa However, the critical distance for thermal weakening decreases strongly with the normal stress (pressure), and thermal weakening will occur more easily at deeper depths. To show this, we use an empirical rule summarized by Di Toro et al. 12 , viz., where a and b are the parameters that depend on the rock type. For all rocks b is positive (for peridotite, a=78, b=1.24 (unit for D th is m, and unit for σ n is MPa)), and this means that the threshold slip distance decreases with confining pressure. Again this is due to the fact that the mechanical work done by friction increases with confining pressure. Niemeijmer et al. 33 conducted high-velocity friction experiments to the normal stress of 40 MPa and obtained slightly different scaling law from that of equation (S-7), but both the relation (S-7) and the one obtained by Niemeijmer et al. 33 predict that shear heating is more efficient at higher normal stress and the critical distance decreases with depth.
Using this scaling, we find that for σ n = 1,000 MPa (depth of ~30 km), D th~1 cm for peridotite. The slip distance (D) for a typical Mw=8 earthquake is on the order of 1 m (for the Sanriku earthquake of 1933, D~3 m 38 ). Therefore we conclude D ≫ D th and thermal weakening occurs easily at a depth deeper than ~10 km.

Scaling for the final friction coefficient ( µ ∞ )
In most of laboratory studies, the thermal weakening is characterized by the reduction in the friction coefficient from µ o ≈ 0.6-0.7 to µ ∞ ≈ 0.1 12 . Thermal weakening could involve various processes including sub-solidus mechanisms as well as melting 12 . Since melting is detected in most cases where thermal weakening is observed 12,31,34 , we will consider a case where the reduction in friction coefficient is caused by melting.
The physical basis for "friction" is (nearly) point contact at the rough surface (e.g., 2,11 ). One can easily derive a linear relation between shear stress and normal stress by assuming that both are controlled by the slip rate-independent "yield strength" of a material, and such a model explains a near universal value of (static) friction coefficient that is independent of temperature (e.g., 2,11 ).
However, when a melt layer is present then the physical basis for using the where A is a constant that depends on the latent heat etc. that is nearly independent of normal stress, R is the characteristic length of melt migration before melt goes to injection vein and W is the characteristic velocity given by W = implying that the resistance for shear decreases with normal stress (pressure).
Therefore the resistance for shear after melting ( µ ∞ ) in the deep lithosphere will be even smaller than the resistance observed in the laboratory at the lower normal stress. Niemeijmer et al. 33 obtained slightly different scaling law from that of equation (S-8) and (S-9), namely τ ∞ ∝σ n 0.5 and µ ∞ ∝σ n −0.5 , but both results predict that shear heating is more efficient at higher normal stress (at σ n = 1,000 MPa, µ ∞ will be ~20-30% of µ ∞ at σ n = 40 MPa).
The reason for the decrease in the "friction coefficient" with normal stress is that the resistance force for fault motion depends on the thickness of the melt layer that decreases weakly (less than linear) with normal stress (pressure). The relation Applicability of the above model to subsolidus mechanisms of thermal weakening is not clear. It is clear, however, that any thermally activated processes likely occur more pronouncedly at higher normal stress, and consequently, it is likely that the reduction in friction coefficient at high normal stress is larger at higher normal stress.
In summary, we conclude that thermal weakening is likely more pronounced in the deep lithosphere and therefore the resistance for high velocity fault motion will be smaller in the deep lithosphere than in the laboratory experiments at low pressures. However, high velocity fault motion occurs only at relatively low temperatures, i.e., in the relatively shallow lithosphere on Earth but nowhere on Venus.

Venus structure
The temperature of the surface of Venus is ~470 C 39 . The atmosphere has very little water 40 although isotopic observations (D/H ratio) suggest that Venus had a lot more water in the past 41 .
The crust of Venus is made mainly of basaltic rocks 42  There are not many constraints on the composition of the mantle. We simply assume that the mantle of Venus has the same composition of Earth. The only difference is water content. The atmosphere of Venus has very little water. And if the crust has as much water as crustal rocks on Earth, its strength would be too low to explain the topography of Venus 44 . And dry crust would imply that mantle would also be dry. So we assume that both the crust and mantle of Venus are dry.
The temperature distribution in Venus is not known mainly because the nature of heat transfer and the degree to which radioactive elements are segregated into the crust are unknown 45,46 . We use the thermal model by Nimmo 47 .
In this paper, we are concerned with the long-wavelength tectonics such as (the absence of) plate boundaries in the stable stage after the catastrophic over-turn, and therefore we use the temperature profile with relatively low thermal gradients (dT/dz)=6 K/km corresponds to the Nimmo-McKenzie model 45 .

Water in the oceanic lithosphere
The oceanic lithosphere is formed at mid-ocean ridges as a residue of partial melting. Consequently, it is essentially dry 48 , i.e., depleted with incompatible elements including hydrogen (water). Indeed, the water content in minerals in the oceanic lithosphere is much less than the saturation limit 49,50 , indicating a low water fugacity.
The hydrothermal activity near mid-ocean ridges hydrates the shallow regions of the oceanic lithosphere. The depth extent to which this hydration occurs can be estimated from the direct sampling by deep ocean drilling 51 and/or using the observed heat flow patterns 52 . Both approaches suggest a depth extent of hydrothermal alteration to ~5 km.
Various models were proposed to suggest deeper hydration down to ~20 km depth or more [53][54][55][56][57] . Some proposed that normal faulting near the trench might bring water into the oceanic lithosphere 53

Dynamics of and stress evolution in a heterogeneous fault
The large reduction in friction coefficient by thermal weakening provides an explanation for why the lithosphere on Earth is weak but the lithosphere on Venus is strong. Also, the large contrast between static and dynamic friction also explains the puzzling observation of large difference between the magnitude of stress drops (~10 MPa or less 64 ) and the magnitude of tectonic stress measured by the in-situ stress measurements (~100 MPa 65 ).
However, such a model contains some subtleties. For example, if the initial strength corresponds to static friction exceeds tectonic stress on the lithosphere, how does a fault motion starts that would lead to thermal weakening? We believe that a key to address this question is a concept that the stress in the lithosphere is heterogeneous caused by the mechanical heterogeneity such as the presence of strong regions (asperities) surrounded by relatively weak regions 66 . Such a concept has strong support from the observations on the source processes of earthquakes 67,68 as well as the statistical nature of seismic activities, e.g., the Gutenberg-Richter (or the Ishimoto-Iida) law 9,69 . In these cases, even if the average stress is lower than the critical stress corresponding to the static friction, there will be regions where the stress is higher than the average stress, and fault motion could start in these regions when the local stress reaches the critical stress needed to overcome static friction.
We consider a fault in which strong asperities are present in the weak background (Fig. 6). Weak regions may be the regions with small asperities or the fault plane covered by a soft material. Strong regions may correspond to large asperities (as will be shown below, a large asperity is more difficult to break than a small one). In this simplified model, slip or creep in the weak regions leads to the stress concentration in the large asperities (strong regions), and the stress in the large asperities increases with time. When the stress in the large asperities reaches a critical value, they will break. This will release the stress there, and from that time on, a similar sequence of stress evolution will occur (e.g., 70,71 ).
The above process can be translated into a semi-quantitative model as follows. Deformation of a large asperity is initially elastic, and the elastic strain in the asperity is produced by deformation of surrounding weak regions. In such a case stress at the asperity (σ asp ) is related to strain ( ε ) as σ asp = C ⋅ε where C is the elastic modulus of the asperity. The strain in an asperity is the ratio of the displacement and the characteristic dimension of an asperity in which deformation occurs. We assume that the characteristic dimension is the same as the size of the asperity, L asp . Then the elastic strain of the asperity is given by ε = υ⋅t L asp where υ is the average velocity of fault motion which is proportional to the velocity of plate motion, υ = ξ ⋅υ plate , where ξ ≤ 1 ( ) is a constant related to the orientation of slip relative to the relative plate motion and t is the time (measured after the last failure of the asperity). Therefore the stress on a large asperity (σ peak ) evolves with time as, When this stress reaches the threshold stress to overcome static friction, fault motion starts, and the stress will be reduced. Using σ asp = 300 MPa, υ plate = 10 −9 m/s (3 cm/y), C=100 GPa 72 , L asp =100 km, the local stress on a large asperity reaches the threshold stress to overcome static friction after ~10 4 years assuming ξ = 1. This agrees with the observations on seismicity of near trench earthquakes 73 .
The relation (S-10) implies that a large asperity is more difficult to break than a small one because elastic stress due to slip (or creep) in the surrounding regions grows more slowly for a large asperity than for a small one. Consequently large earthquakes occur less frequently than small ones. However, a model summarized by equation (S-10) contains a few simplifying assumptions. One of which is the assumption that stress evolution is quasi-static: deformation in the weak region occurs by creep and deformation in the strong region is elastic until threshold stress is reached. In this case, stress evolution at a large asperity is (approximately) linear with time, and one would expect a large stress drop. If deformation in the weak region involves a number of small unstable slips (small earthquakes), then dynamic interactions will play a role that would likely reduce the threshold (initial) stress at which a strong asperity will yield, and hence reduces the stress drop as is shown below 74,75 .
A key in the dynamic interaction is the stress concentration at a propagating crack tip. Yamashita 75 investigated the dynamics of faulting in a region with distributed strength. He assumes that slip in the weak regions occur by small-scale unstable fault (crack) propagation. In such a case, stress concentration at the tip of propagating small faults enhances the stress at a large asperity. As a small fault (crack) gets close to a large asperity, stress on the large asperity increases rapidly with time until the local stress reaches the threshold stress against static friction (σ peak ) (Fig. S2) (a more detailed study was made by Noda et al. 74 ).
In this model, a large asperity will be broken when the initial stress reaches ( ) represents the role of stress concentration) rather than σ initial = σ peak . This stress amplification occurs only when the tip of a small fault approaches close to a large asperity. In these cases, the stress at a large asperity is near the peak stress only for small displacement, and hence the peak stress does not contribute to the seismic moment very much. Therefore, the stress drop observed by seismology is not σ peak − σ final but Δσ = σ initial − σ final = σ peak ζ − σ final 69,75 . This explains relatively small stress drops (less than ~10 MPa) 64,76 for earthquakes compared to σ peak − σ final . In this case, the recurrence time of earthquakes will also be modified.
To understand how such a model works, we performed a numerical simulation of stress evolution of a heterogeneous fault (Fig. S3). A 60º dipping normal fault is loaded at a constant rate of 10 -9 m/s with unstable velocityweakening properties down to 20 km depth. This contains a weak region with a static friction coefficient of 0.2 sandwiched by two strong regions with a static friction coefficient of 0.6. Both regions are assumed to follow velocity-weakening behavior. Following a standard formulation, we characterize the state-and-rate friction by two parameters "a" and "b" (see, e.g., 9 ), and a-b=-3x10 -2 for a strong region with strong dynamic weakening, while a-b= -4x10 -3 for weak regions with normal dynamic weakening.
We simulate the evolution of fault slip using the boundary integral method with the radiation approximation (e.g., [77][78][79]. As we neglect the radiation of seismic waves, our estimates of strong weakening are conservative 80 . The evolution of stress on such a fault is shown in Fig. S3. This model provides a reasonable explanation for how the asperities with the high static friction coefficient can be broken at a small overall stress as well as how often these events (earthquakes) could occur, and is consistent with the statistical nature of earthquakes. This model also explains the relatively small stress drops observed for the intra-plate earthquakes despite the high static friction of rocks in the laboratory.
However, details of the interaction between faults with a variety of asperity sizes are not fully understood, including the magnitude of the stress amplification factor, ζ = σ peak σ initial , and its dependence on the nature of small-scale faulting. This problem is also related to the long-standing issue of discrepancy between the tectonic stress estimated from the in-situ stress measurements 65   Due to the stress accumulation caused by slip (or creep) in the weak region, stress increases at a large asperity. When the local stress at an asperity reaches a critical value, faulting (a large earthquake) occurs. If slip in the weak region occurs through small-scale unstable slip (crack propagation), then stress evolution on a large asperity is complex (see the inset; after Yamashita 75 (see also Noda et al. 74 )): When the tip of a small fault approaches a large asperity, stress increases from σ initial to σ peak rapidly due to the stress concentration at the tip of the propagating crack. The stress drop observed by seismology is Δσ = σ initial − σ final 69,75 that is much smaller than σ peak − σ final .  Table S1. Flow laws used in constructing Fig. 3 and Fig. 5 deformation mechanism flow law (1) references, notes Diffusion creep Power-law dislocation creep Peierls mechanism (2) 44,85,86 87 (1) ! ε : strain-rate, A 1,2,3 : constants, σ : differential stress (σ 1 − σ 3 ), L: grain size, n, m, q, r: constants, E 1,2,3 * : activation energy, V 1,2,3 * : activation volume, P: pressure, T: temperature, R: the gas constant, σ P : the Peierls stress (2) The activation volume for diffusion creep has not been well constrained. We assume a tentative value of 5 cc/mol (this does not affect the results compared to the uncertainties of grain size).