Effects of rheological stratification and elasticity of lithosphere on subduction initiation

The breaking and bending of rigid and elastic lithosphere was probably essential for the initiation of plate subduction, although how this occurred is still poorly understood. Here we test effects of rheological stratification and elasticity of lithosphere on subduction initiation, which are possibly resulting from thermal cracking and seawater penetration into the lithosphere. In addition to the strong influence of water on rheological behavior, the material rigidity is also sensitive to the development of crack and fluid saturation. Numerical modeling indicates that water-weakening and a low-rigidity lithosphere are essential for the initiation of plate subduction, and such conditions are likely to have arisen on the early Earth due to extensive thermal contraction of the planet. Our results indicate that the formation of thermal cracks and penetration of seawater play an important role on subduction initiation, and are likely to have operated on planets other than Earth. However, if the ocean is disappeared, fluid penetration is likely to cease and plate tectonics would have stopped due to increasing the strength and rigidity of the lithosphere.


Introduction
Plate tectonics facilitates the recycling of materials, including carbon and water, into the Earth's interior, and this has acted to stabilize the climate and nurtured biological evolution throughout geological history (e.g., Korenaga, 2013;Stern and Gerya, 2018). Although the physical mechanism that underlies plate tectonics is simple, being driven by the negative buoyancy of cold and dense lithosphere, the initiation of plate subduction requires the bending and breaking of strong lithosphere, which is difficult to achieve. The development of a weak zone within the lithospheres has been discussed in previous studies (e.g., Regenauer-Lieb et al., 2001;Hall et al., 2003;Bercovici and Ricard, 2014;Gerya et al., 2015); however, a major barrier to the initiation of subduction is achieving the flexural bending of rigid lithosphere that behaves elastically over geological timescales. Flexural rigidity is an indicator of bending capability and is controlled mainly by the elastic OPEN ACCESS EDITED BY John Hernlund, Tokyo Institute of Technology, Japan properties and rheological structure of the lithosphere (Watts et al., 2013;Turcotte and Schubert, 2014), indicating it is highly sensitive to crack development and fluid-rock interaction.
During the early history of Earth, rapid cooling due to the formation of an ocean likely resulted in the extensive formation of thermal cracks within the lithosphere. Simple viscoelastic models suggest that the accumulation of thermal stress could have created tensional cracks as deep as 50 km (Korenaga, 2007). These fractures could have facilitated the penetration of seawater, reduced the rigidity of the lithosphere, and enabled the development of weak zones that became plate boundaries ( Figure 1). As a modern analogue, the plate-bending faults at close to trench cause fluid infiltration and serpentinization of the oceanic lithosphere that are inferred from low seismic velocity and high electrical conductivity (e.g., Fujie et al., 2013;Naif et al., 2015;Grevemeyer et al., 2018). A significant reduction of plate rigidity is suggested by the flexure modeling that simulates the bathymetric profiles toward the trench (Contreras-Reyes and Osses, 2010). However, the present-day hydration mechanism is a consequence of plate subduction, which cannot apply to the onset of subduction on the early Earth. Such that the extensive thermal contraction is one of the plausible mechanisms to reduce the strength and rigidity of the lithosphere. In this study, we test this hypothesis based on rock physics and thermomechanical modeling, examining how plate rigidity and the style of mantle convection are influenced by thermal cracking and fluid penetration.
2 Plate strength model 2.1 Development of thermal cracks and its influence on elastic properties Thermal cracks form when the thermal stress exceeds the material strength. Such cracks can develop over a wide range of scales, from the grain scale as a result of the different thermal expansivities of minerals (Kranz, 1983), to the geological scale due to cooling of the lithosphere (Turcotte, 1974). Given a thermal expansivity of 10 -5 K −1 , a bulk modulus of 100 GPa, and temperature difference of 1,000°C, thermal stress can be created as high as 1 GPa, which is sufficiently large to fracture geological materials. Although such stress could be dissipated by viscous relaxation, the cold parts of lithosphere are strong enough to accumulate thermal stress according to the temperature-dependent viscosity (Korenaga, 2007). Such development of thermal cracks can change the elastic properties of the lithosphere.
The effective medium theory indicates that elastic properties of rocks are strongly dependent on crack density and crack geometry (e.g., Kuster and Toksöz, 1974;O'Connell and Budiansky, 1974;Kachanov et al., 1994). Figure 2 shows how Young's and shear moduli change with crack porosity and crack aspect ratio. Laboratory experiments indicate the aspect ratio of~10 -3 for thermal cracks developed in granitic rocks (Wang et al., 2013), suggesting that even small amounts of thermal crack can drastically decrease the elastic properties of materials ( Figure 2). Although a large extrapolation is needed for the scale of cracks between laboratory and geological events, application of the effective medium theory indicates that a marked decrease in elastic properties is expected in a thermally damaged lithosphere.

Calculation of the plate strength
The plate strength is mainly controlled by brittle deformation at shallow depths, in which the strength is determined from the resistance to frictional sliding. The shear stress at which frictional sliding begins is proportional to the frictional coefficient and the normal stress applied to the fault plane. If faults are filled by aqueous fluids, pore fluid pressure reduces the effective normal stress acting on a fault plane, resulting in a decrease of the mechanical strength of the lithosphere. In the calculation of wet model, we assumed a hydrostatic pore pressure (pore fluid pressure ratio of λ = 0.3). Frictional coefficient is nearly constant for various materials under dry condition; however, hydrous phyllosilicates showed a markedly low frictional coefficient under fluid-bearing conditions (Moore and Lockner, 2007;Katayama et al., 2015;Hirauchi et al., 2016). Such that we calculated the brittle strength using frictional coefficient of µ = 0.6 for the dry model and µ = 0.1 for the wet model. These two effects, pore fluid pressure and low frictional coefficient, cause a marked decrease in the brittle strength of the lithosphere at fluid-saturated environment.
At deeper depths, plastic deformation becomes dominant as temperature increases with depths, and the plastic strength is calculated using flow laws. A commonly used flow law is the power-law relationship between strain rate ( _ ε) and stress (σ) as, where A is a pre-exponential factor, d is grain size, H is activation enthalpy (H = E + PV, where E is activation energy, V is activation volume and P is pressure), T is temperature, and R is the gas constant. If deformation is controlled by diffusion creep, the stress and grain-size exponents are n = 1 and m = 3, respectively, whereas dislocation creep is insensitive to grain-size and has the stress exponent of n~3 (e.g., Kohlstedt et al., 1995). This relationship is commonly used to calculate the plastic strength of lithosphere; however, the power-law creep has limitations under certain conditions (Tsenn and Carter, 1987). Under high stress, the stress dependence of activation enthalpy becomes important, and the appropriate flow law is changed to where σ P is the Peierls stress. This mechanism, known as Peierls creep, becomes dominant at depths close to the brittle-plastic transition (e.g., Katayama, 2021). Note that the stressdependence of activation enthalpy given by this relationship is only applicable to dislocation in a forward motion, whereas dislocation in a backward motion is needed to extrapolate to lower stresses (Poirier, 1985), such that the Peierls creep can be rewritten as,

FIGURE 2
Effect of crack on elastic properties of rocks. Changes in Young's modulus (A) and shear modulus (B) as a function of crack porosity (%) with different aspect ratio α, calculated from the effective medium theory of Kuster and Toksöz (1974).
(3) Figure 3 shows the plastic strength of each creep mechanism including the effect of backward motion in the Peierls creep. Plastic deformation is generally controlled by the weakest constituent mineral, and we used flow laws of plagioclase and olivine to calculate the plastic strength of crust and mantle, respectively. Given that water has a strong influence on plastic deformation, the flow parameters determined under dry and fluid-saturated conditions were applied for the dry and wet models (Table 1). To calculate the strength profile, we assumed a constant strain rate of 10 -15 , 10 -16 , or 10 -17 s −1 and a grain size of 1 mm, noting that these values have little influence on the lithospheric strength ( Figure 4). We also applied the flow parameters determined by the Bayesian statistics (Jain et al., 2017), which gives a slightly lower yield stress, but the overall shape of the strength envelope does not change.

Calculation of the flexural rigidity
On the geological time scales, the lithosphere behaves as an elastic layer. However, a purely elastic behavior, i.e., linear stress distribution, is unrealistic because the constituting rocks are deformed by frictional sliding in the upper parts and by plastic creep in the lower parts as illustrated in the strength envelope of lithosphere ( Figure 4). Such that the flexure of the lithosphere has been discussed using a rheologically layered model; in this case flexural rigidity is obtained by numerical integration of the bending moment in the lithosphere (e.g., Goetze and Evans, 1979;McNutt, 1984). The bending moment must be balanced between tensional and compressional forces in the lithosphere, and is given by, where T m is the mechanical thickness of the lithosphere, σ(z) is the strength at depth z, and z n is the depth of the neutral stress plane where the stress field switches from tension to compression or vice versa. The analytical solution for the bending moment of lithosphere is expressed as, where E is Young's modulus, ] is Poisson's ratio, K is flexural curvature, and h is elastic thickness (McNutt, 1984). This relationship is used to calculate the flexural rigidity D, which is expressed as the bending moment divided by the curvature (D = M/K). Because the bending moment of lithosphere is strongly dependent on the mechanical strength at depths close to the brittle-plastic transition, the incorporation of Peierls creep makes a large contribution to the estimate of flexural rigidity. Mechanical strength in both brittle and plastic regimes is highly dependent on the presence of water, such that a marked change in the bending moment is expected between the dry and wet models. In addition to strength weakening by the penetration of seawater, the presence of thermal cracks can markedly change the elasticity of the lithosphere. Consequently, the flexural rigidity is strongly reduced in the wet model ( Figure 5), indicating that the flexure of an elastic layer would be substantially enhanced in a thermally damaged lithosphere, thereby promoting the initiation of plate subduction.

FIGURE 3
Relationship between stress and temperature for wet olivine calculated using a constant strain rate of 10 -15 s −1 (A). Peierls creep dominates at low temperatures. Relationship between strain rate and stress for wet olivine calculated using a temperature of 700°C (B). Calculation without backward direction results in a crossover of Peierls at low stress.
Frontiers in Earth Science frontiersin.org 04 3 Numerical modeling of plate subduction We have tested how strength weakening and flexural rigidity of the lithosphere influence the style of plate subduction using thermomechanical numerical simulations. The computational domain of the numerical model was filled by a fluid with visco-elasto-plastic rheology in 2-D Cartesian geometry (e.g., Liao and Gerya, 2014;Yoshida et al., 2020). The width and depth of the model were 750 km and 310 km along the horizontal (x) and vertical (z) directions, respectively ( Figure 6). The thermo-chemical convection of the fluid was solved using the Rheological parameters are mainly from Katayama (2021), but dislocation creep for plagioclase is referred to Rybacki and Dresen (2000) and Peierls creep for plagioclase was recalculated.

FIGURE 4
Strength profiles of lithosphere for dry (A) and wet (B) models. Plastic strengths for variable thermal gradients calculated using different strain rates of 10 -15 , 10 -16 and 10 -17 s −1 . A flexural curvature of 5×10 -7 m −1 was assumed, and positive and negative stresses correspond to tension and compression, respectively. The bending moment was calculated via a balance between tensional and compressional forces within the flexed lithosphere (shaded area). The depth of Moho was assumed to be 20 km, in which the crustal thickness could be thicker during Archean than that of the present-day oceanic crust.
Frontiers in Earth Science frontiersin.org 05 finite-difference and a marker-in-cell methods. The grid interval of the computational grid was uniform through the entire domain and set as 2 km for the x and z directions. The number of tracer markers in the entire domain was 750 × 310 for the x and z directions, respectively, which is equivalent to an average tracer density of 1 km in both directions.
The boundaries of the computational domain were set as being shear-stress free and impermeable. Isothermal boundary conditions were imposed on the top and bottom surface boundaries (0°C and 1700°C, respectively), and a symmetric condition was imposed on the side boundaries. Considering that the reference temperature of the Archean mantle is about 200°C higher than that of the present-day mantle (Korenaga, 2013), the potential temperature at the top surface was set to 1,550°C. The initial condition of the thermal state above the bottom of the lithosphere was based on the half-space cooling model (Turcotte and Schubert, 2014), and the adiabatic thermal gradient in the upper mantle was set to 0.5°C km −1 .
We set up a calculation domain containing shallower and deeper crustal layers, the lithosphere, and the upper mantle with a localized fracture zone caused by extensive cracking and water penetration ( Figure 6D). A sticky water layer with a thickness of 10 km and a constant viscosity of 10 18 Pa s was included in the topmost part of the computational domain to enable free surface deformation (e.g., Schmeling et al., 2008). The viscosity except for the sticky water layer was determined by the flow law parameters shown in Table 1 and was truncated between η um and 10 26 Pa s, where η um is a free parameter that controls the upper mantle viscosity. The physical properties of each rock type in the numerical model are listed in Table 2.
For the visco-elasto-plastic rheology used in the numerical modeling, the deviatoric strain rate _ ε ij is the sum of the three components as (e.g., Gerya, 2019), where _ ε viscos ij , _ ε elastic ij , and _ ε plastic ij are the viscous, elastic, and plastic strain rates, respectively. The viscos strain rate is given by the sum of strain rates under independent creep mechanism as, where the dislocation power-law creep and Peierls creep are both related to the dislocation motion and do not operate simultaneously; therefore, the faster one is a rate-controlling mechanism (e.g., Frost and Ashby, 1982). The shallower crust is assumed to be fluid-saturated conditions with µ = 0.1 and λ = 0.9, whereas the deeper crust and mantle are assumed to be dry conditions with µ = 0.6. The cohesion of the lithosphere and upper mantle was set to 50 MPa. As seawater can penetrate along the fracture networks easily, the wet rheology with a low frictional coefficient was used for the fracture zone. In addition to the water weakening, we tested how the style of subduction is affected by the elasticity of the lithosphere, using different shear moduli G/G 0 = 1 to 0.1, where G 0 is the reference shear modulus. Note that both G and E are both sensitive to the development of cracks, and the changes in G show a similar trend as in E (Figure 2). Figure 7 shows the results of numerical modeling for the wet model with different plate rigidity. The upper mantle viscosity η um was assumed to be 2×10 20 Pa s. Even though water weakening in the fracture zone is considered in the models, bending of the lithosphere does not occur in the model with strong rigidity (G/G 0 = 1); consequently, plate subduction is not initiated even after a long duration ( Figure 7A). However, for the models with G/G 0 = 0.5, the viscoelastic instability overcomes the barrier of flexural bending, resulting in subduction initiation  Figure 4. In the wet model, the shear modulus was assumed to be half that of the dry model. Solid lines represent calculations using a strain rate of 10 -16 s −1 , and dashed lines show variations in strain rates between 10 -15 and 10 -17 s −1 .

Frontiers in Earth Science
frontiersin.org 06

FIGURE 6
Initial condition for the numerical simulation. (A) Composition field (B) temperature field (C) viscosity structure, and (D) a close-up view of the area enclosed by the dashed red rectangle in (A). The computational domain comprised shallower and deeper crustal layers with the same thickness of 10 km and the underlying lithosphere with a variable thickness: the depths of the bottom of the right and left lithospheres were set to 100 km and 50 km, respectively. White contours represent isotherms with a contour interval of 400°C.

Frontiers in Earth Science
frontiersin.org 07 within the weak fracture zone ( Figure 7B). Subduction is further facilitated in the model with a lower rigidity (G/G 0 = 0.2) ( Figure 7C). Note that for the model with G/G 0 = 1 ( Figure 7D), no further initiation of subduction was observed even after a sufficiently long time (~50 Myr), during which secondary downwelling plumes are fully developed owing to the thermal instability at the base of the stagnant lithosphere (e.g., Solomatov and Moresi, 2000).
Here, to see the effect of wet condition in the shallower crust on the subduction initiation, we have performed additional models in which the shallower crust is composed of the dry rheology same as the deeper crust, keeping the fracture zone under wet condition (Figure 8). The upper mantle viscosity η um was assumed to be 2×10 20 Pa s same as the wet model. In the case that the shallower and deeper crusts have the dry condition, plate subduction did not occur even for G/G 0 = 0.2 ( Figure 8C), which suggests that water penetration into the shallower crust is essential for plate subduction. In such a dry model, the behavior of the stagnant lithosphere had little effect on G/G 0 , because little deformation of the fracture zone occurs; consequently, plate subduction did not occur after a long duration even for the model with G/G 0 = 0.2 ( Figure 8D).
On the other hand, in the wet model, the subduction style is roughly classified into three regimes as a function of G/G 0 , and the upper mantle viscosity η um (Figure 9). Our numerical results demonstrate that the possibility of subduction depends also on the upper mantle viscosity because it is achieved by a balance between the negative buoyancy of subducting plate and the viscous resistance acting on the edge of subducting plate. As the rigidity of lithosphere is reduced, the range of the upper mantle where subduction occurs tends to increase. The upper mantle viscosity for realizing the plate subduction is 3×10 20 Pa s or lower for the models with G/G 0 ≥ 0.2, which is comparable to the average viscosity of the upper mantle obtained from a recent study on glacial isostatic adjustment (Argus et al., 2021).

Possible existence of damaged lithosphere
Our results indicate that plate subduction is facilitated by strength weakening and low flexural rigidity, which can be accomplished by the development of thermal cracks in the lithosphere and the penetration of seawater into the fracture zone. Thermal cracks were possibly widespread on the early Earth due to a large thermal contrast in the lithosphere; however, it is difficult to prove the existence of damaged lithosphere in the past. The present-day oceanic lithosphere exhibits intraplate seismicity depths down to~50 km (e.g., Wiens and Stein, 1983). These events are away from the plate boundaries, and are possibly related to contraction cracks that relieve the thermal stress in the cooling lithosphere (Turcotte, 1974). In the ophiolite sequences where the paleo-oceanic plates were emplaced, the developments of microcracks in the deeper sections were likely associated with anisotropic thermal contraction during cooling (e.g., Nicolas et al., 2003). These crack networks would be the main mechanisms favoring pervasive seawater penetration into the deep oceanic crust possibly down to the mantle. In the regions close to the trench, low seismic velocity and high electrical conductivity were observed in the relatively old oceanic lithosphere, suggesting the seawater penetration into the mantle depths along the fracture zones (e.g., Fujie et al., 2013;Naif et al., 2015;Grevemeyer et al., 2018). Although these are associated with the tensional stress due to the plate bending, the elastic properties can be influenced by the fracture developments and water infiltration. In fact, the significant reduction of flexural rigidity is indicated at close to the trench, most likely caused by the weakening of the oceanic lithosphere due to the hydration along the bendingrelated faults (Contreras-Reyes and Osses, 2010). ρ is density kgm −3 , k is thermal conductivity, W m −1 K −1 , H is the radiogenic heat production rate per unit mass, W kg −1 , k 0 is the bulk modulus, GPa; G 0 is the shear modulus, GPa. Thermal conductivity was calculated as k=[k 0 +k a /(T+77)](1+k b P), where T is temperature (K) and P is pressure (MPa) (Clauser and Huenges, 1995;Hofmeister 1999). The density of the crustal layers is assumed to be 3200 kg m −3 , based on the composition of Archean komatiite (Arndt et al., 2008). The values of H, K 0 and G 0 from Turcotte and Schubert. (2014) and those of k 0 , k a , and k b are from Gerya. (2019). For all rock types, the coefficient of thermal expansivity α is 3.0 × 10 −5 K −1 and the specific heat at the constant pressure c p is 1250 J kg −1 K −1 . For the sticky water layer (Stacey and Davis, 2008), ρ = 1025 kg m −3 , α = 1.5 × 10 −4 K −1 , c p = 3990 J kg −1 K −1 , and k = 0.59 W m −1 K −1 . Gravitational acceleration is 9.81 m s −2 and the gas constant is 8.314 J mol −1 K −1 .
Frontiers in Earth Science frontiersin.org 08

Subduction initiation on early earth and potential on other planets
It is widely accepted that the key process for subduction initiation is the formation of a weak plate boundary, although the weakening mechanism remains controversial (e.g., Regenauer-Lieb et al., 2001;Hall et al., 2003;Bercovici and Ricard, 2014;Gerya et al., 2015). Plate tectonics operates only on Earth, which differs from the other terrestrial planets in having liquid oceans, suggesting that water could play an important role in subduction initiation. However, we must avoid explaining the initiation of plate subduction by a hydration process that itself requires the operation of plate tectonics (Lu et al., 2021). In this sense, thermal cracking is the most plausible mechanism of deep fluid penetration and weakening of the lithosphere (Korenaga, 2007). Our results indicate that weakening coupled with a low flexural rigidity is essential for the initiation of plate subduction; in contrast, most previous models have considered the weakening effect but not changes in plate rigidity (e.g., Regenauer-Lieb et al., 2001;Hall et al., 2003;Bercovici and Ricard, 2014;Gerya et al., 2015). Thielmann and Kaus (2012) considered that

FIGURE 7
Results of numerical modeling for wet condition with different values of plate rigidity. (A-C) Snapshot of simulation results using (A) G/G 0 = 1 (B) G/G 0 = 0.5, and (C) G/G 0 = 0.2 at 20 Myr after the start of the simulation. From top to bottom, the composition field, the temperature field, and the viscosity structure are shown. Plate subduction occurs in cases (B) and (C,D) Snapshot of simulation results using G/G 0 = 1 at 30, 40, and 50 Myr after the start of the simulation. The upper mantle viscosity η um was assumed to be 2×10 20 Pa s.

Frontiers in Earth Science
frontiersin.org 09 subduction initiation was caused by buckling and shear heating associated with a low shear modulus, although these authors did not discuss the mechanism that leads to reduce plate elasticity.
We propose that both thermal contraction and seawater penetration are required for the operation of plate tectonics, and that these conditions were possibly achieved in the other terrestrial planets during their early stage of history. On Mars, the former presence of a global ocean has been inferred from geomorphology (Baker et al., 1991) and widespread occurrence of clay minerals (Ehlmann et al., 2011). The hemispheric dichotomy of the Martian topography might be attributed to the ancient operation of plate tectonics (Zuber, 2001), which might have ceased as the ocean disappeared. A similar situation might have occurred on Venus, where an ancient ocean is expected to have existed, if the initial water endowment was similar to that of Earth (Kasting, 1988). Thus, the plate tectonics may be a ubiquitous feature of rocky planets in their early history, but its operation timescales are dependent on the presence of liquid ocean.

FIGURE 8
Results of numerical modeling for dry condition in both shallower and deeper crust. (A-C) Snapshot of simulation results using (A) G/G 0 = 1 (B) G/G 0 = 0.5, and (C) G/G 0 = 0.2 at 20 Myr after the start of the simulation. Plate subduction does not occur in any cases (D) Snapshot of simulation results using G/G 0 = 0.2 at 30, 40, and 50 Myr after the start of the simulation. The upper mantle viscosity η um was assumed to be 2×10 20 Pa s same as the wet model.

Frontiers in Earth Science
frontiersin.org Frontiers in Earth Science frontiersin.org