How Hydrothermal Cooling and Magmatic Sill Intrusions Control Flip‐Flop Faulting at Ultraslow‐Spreading Mid‐Ocean Ridges

“Flip‐flop” detachment mode represents an endmember type of lithosphere‐scale faulting observed at almost amagmatic sections of ultraslow‐spreading mid‐ocean ridges. Recent numerical experiments using an imposed steady temperature structure show that an axial temperature maximum is essential to trigger flip‐flop faults by focusing flexural strain in the footwall of the active fault. However, ridge segments without significant melt budget are more likely to be in a transient thermal state controlled, at least partly, by the faulting dynamics themselves. Therefore, we investigate which processes control the thermal structure of the lithosphere and how feedbacks with the deformation mechanisms can explain observed faulting patterns. We present results of 2‐D thermo‐mechanical numerical modeling including serpentinization reactions and dynamic grain size evolution. The model features a novel form of parametrized hydrothermal cooling along fault zones as well as the thermal and rheological effects of periodic sill intrusions. We find that the interplay of hydrothermal fault zone cooling and periodic sill intrusions in the footwall facilitates the flip‐flop detachment mode. Hydrothermal cooling of the fault zone pushes the temperature maximum into the footwall, while intrusions near the temperature maximum further weaken the rock and promote the formation of new faults with opposite polarity. Our model allows us to put constraints on the magnitude of two processes, and we obtain most reasonable melt budgets and hydrothermal heat fluxes if both are considered. Furthermore, we frequently observe two other faulting modes in our experiments complementing flip‐flop faulting to yield a potentially more robust alternative interpretation for existing observations.


Introduction
Detachment faulting is a unique spreading mode observed at magma-poor slow-to ultraslow-spreading midocean ridges (MOR).It is a highly asymmetric process (MacLeod et al., 2009;Tucholke et al., 2008), where the hanging wall side consists of volcanic crust with moderate normal faulting, while the footwall side is purely tectonic with a long-lived detachment fault rooting in the weak near-axis lithosphere.During plate separation, the breakaway-the tip of the footwall, where the fault first cuts the seafloor-migrates off-axis (Buck et al., 2005).With increasing displacement on the fault plane, the uplifted footwall starts rotating and bending under its own weight (MacLeod et al., 2009;Smith et al., 2006).The domed corrugated shear planes of such detachments, known as oceanic core complexes, emerge to the seafloor at shallow angles where they expose lower crust and mantle rocks (e.g., Escartín et al., 2008Escartín et al., , 2017;;Hayman et al., 2011;Zhao et al., 2013).
In this study, we focus on an end-member type of detachment faulting observed along almost amagmatic sections of ultraslow MORs, such as the narrow corridors of highly serpentinized seafloor between 62°E and 65°E at the Southwest Indian Ridge (SWIR).In the absence of melt supply, extension in the hanging wall of a detachment cannot be accommodated by magmatic accretion.The entire fault zone of the active detachment therefore migrates across the ridge axis into thicker lithosphere until the formation of a new on-axis fault in the footwall is energetically more favorable.This new fault is most likely dipping in the opposite direction (i.e., is of opposite polarity) and is oriented along weak zones in the footwall resulting from strain induced by its rotation and flexure.Reston and McDermott (2011) proposed this so-called "flip-flop" detachment mode to be responsible for large areas of exhumed, extensively serpentinized mantle rocks at magma-poor margins and ultraslow-spreading MORs.This extension mode produces regions of unusually smooth seafloor, large ridges that are oriented parallel to the spreading axis, and a pattern of off-axis dipping fault zones that is symmetrical on a large-scale (Cannat, Sauter, et al., 2019;Cannat et al., 2006;Sauter et al., 2013).
A recent numerical study by Bickert et al. (2020) explores the processes causing the flip-flop kinematics.They find that fault zone weakening through mantle serpentinization and grain size reduction in mechanically stressed regions in the lower lithosphere play fundamental roles in stabilizing active tectonic zones and triggering new faults of opposite polarity.A crucial assumption in the aforementioned model is that the near-ridge lithosphere is in a thermal quasi-steady state with prescribed temperature maximum at the ridge axis and fixed depth-dependent temperature profiles.However, gaining further insights into the multifaceted thermo-mechanical interplay of the processes involved in this spreading mode requires a dynamically evolving temperature field.In this study we therefore investigate which processes predominantly shape the near-ridge thermal structure and how the evolving temperature field, in turn, affects the faulting sequence through rheological changes.
The thermal state of (ultra)slow amagmatic ridges is not as well understood as at faster spreading ridges with strong magmatism.A nearly continuous axial melt lens with minor depth variation, for example, at the fastspreading East Pacific Rise (Detrick et al., 1987;Marjanović et al., 2018), indicates a rather stable thermal regime in both space and time.The hot, weak axial region serves as a thermal and mechanical anchor that stabilizes the ridge axis position and focuses deformation.In contrast, amagmatic ridge sections without a significant melt budget are in a more transient thermal state that is controlled, at least partly, by its faulting dynamics and history (e.g., Behn & Ito, 2008).Here, in the absence of strong magmatic processes, advection of hot mantle material into the footwall of detachments and subsequent cooling by heat conduction and hydrothermal circulation become the main heat distribution mechanisms.
Indirect evidence for hydrothermal cooling is the more than 15 km thick axial brittle lithosphere at ultraslowspreading ridges (e.g., Grevemeyer et al., 2019), which is incompatible with solely conductive cooling (see also our model without hydrothermal cooling; Section 3.2.1).Mapping of thermal plumes in the water column point at a reduced yet still significant hydrothermal activity along sections of the ultraslow Gakkel Ridge and SWIR (Baker et al., 2004).Direct observations of hydrothermal venting in amagmatic settings are limited to the recently discovered low temperature Old City hydrothermal field in the eastern magma-poor corridor of the SWIR (Cannat, Agrinier, et al., 2019;Lecoeuvre et al., 2021).Additional evidence of potentially widespread hydrothermal activity is the inactive ultramafic-hosted Tianzuo hydrothermal field (Ding et al., 2021), located on a detachment surface in a smooth seafloor area between the two corridors investigated by Sauter et al. (2013).Further indicators are the seismically imaged ∼5 km thick layer of partially serpentinized material exhumed by the successive detachments (Corbalán et al., 2021;Momoh et al., 2017Momoh et al., , 2020;;Sauter et al., 2013) and conclusions of other petrological studies (Bickert et al., 2023;Patterson et al., 2021;Vieira Duarte et al., 2020).Observed correlations between fault and vent field positions at more magmatic ridge sections (e.g., TAG, de Martin et al. (2007); Logatchev-1, Petersen et al. (2009); Longqi-1, Tao et al. (2020)) suggest that hydrothermal circulation presumably occurs along highly permeable fault zones (McCaig et al., 2007).Potential driving forces for hydrothermal activity at magma-poor ridges are exothermic heat from serpentinization reactions, magmatic heat from sporadic intrusions into the footwall (indicated in recent seismic data; Momoh et al., 2020), and possibly heat mined from deep hot rocks through thermal contraction cracks (Lister, 1974).Fan et al. (2021) suggest that the observed lithosphere thickness results from a thermal balance between deep periodic sill intrusions and hydrothermal circulation in the upper lithosphere driven by their heat.More details on spreading-related processes including hydrothermal activity can be found in the excellent reviews by Olive (2023) and Früh-Green et al. (2022).

10.1029/2023GC011331
Our study is motivated by these recently published data and models, and associated advances in the understanding of the structure and thermal regime of this complex end-member tectonic setting.The key question is, how can the interplay of tectonic, magmatic and hydrothermal processes explain a thermal structure that will lead to flip-flop detachment faulting during ultraslow, magma-poor oceanic spreading.To address this question, we present a thermo-mechanical numerical model for amagmatic oceanic spreading that includes refined implementations of serpentinization and grain size evolution, parametrized hydrothermal cooling (enhanced along active fault zones) and the thermal and rheological effects of sill intrusions into the detachment footwall.Starting with a setup similar to the model by Bickert et al. (2020) with an imposed temperature structure, we increase complexity by calculating the dynamically evolving thermal structure including the effects of hydrothermal fault zone cooling and magmatic intrusions.

Numerical Method and Constitutive Equations
The numerical model used in this study applies the finite element method on unstructured meshes to solve for the thermal and mechanical evolution of lithosphere and underlying mantle.The program M2TRI_vep (Hasenclever & Glink, 2024, see Data Availability Statement for more information) is a successor of the mantle convection code M2TRI (Hasenclever, 2010), into which elastic and plastic deformation behavior, serpentinization of mantle rocks and a simple hydrothermal cooling parametrization have been implemented as discussed in previous studies (Hasenclever et al., 2017;Rüpke & Hasenclever, 2017).For the present study, we have further advanced the model by improving the visco-elasto-plastic rheology formulation, which now considers diffusion creep as well as non-linear dislocation creep mechanisms and a viscoplastic regularization (Duretz et al., 2021).Furthermore, we have implemented processes that are critical for lithosphere and upper mantle deformation such as strain softening, deformation-controlled serpentinization, dynamic grain size evolution, periodic emplacement of magmatic intrusions and an advanced hydrothermal cooling parametrization that takes into account the deformation state of the lithosphere.In the following, the relevant equations are written using index notation and Einstein summation convention.All parameters and variables introduced in this section are listed in Table 1.
We solve the conservation equations for mass, Equation 1, and momentum, Equation 2, for an incompressible medium using the Boussinesq approximation.The solution variables are displacement velocity components u i and total pressure P: where x i are the spatial coordinates, τ ij are the deviatoric stresses, ρ is the density and g i are the components of the gravitational acceleration.
Density in the Stokes Equation 2 is assumed to vary with temperature T and degree of serpentinization s.The latter is used to calculate the volumetric mean between two end-member rock types: mantle peridotite (for which we assume olivine properties; subscript ol) and serpentine (subscript serp).This averaging is done for all rock parameters for which two values-one for olivine and one for serpentine-are given in Table 1.
Densities ρ c,0 , with c ∈ {ol, serp}, are defined at the reference temperature T 0 = 0°C and α c are the thermal expansion coefficients.Note that volume changes associated with serpentinization as well as other volume changing processes such as the formation of a fracture network (cf.Section 2.2.3) and magmatic sill emplacement (cf.Section 2.4) are disregarded in this incompressible formulation.

10.1029/2023GC011331
The conservation equation of heat is solved to obtain the temperature T:  where C p is the specific heat capacity and κ is the thermal conductivity.H-terms represent the heat sources and sinks introduced in the next sections, that is, viscous dissipation H vd , exothermic heat of serpentinization H serp and latent heat of crystallization H lat .Radiogenic heat production is neglected in this study, since for the case of amagmatic spreading investigated here, we assume the lithosphere will mainly consist of mantle rocks with a low concentration of radioactive elements.

Visco-Elasto-Plastic Behavior
On geological time scales, the mechanical deformation of lithosphere and upper mantle rocks occurs through three different mechanisms: reversible elastic deformation, irreversible viscous flow and irreversible plastic deformation simulating brittle failure.All three mechanisms act at strongly variable quantitative proportions depending mainly on the ambient temperature-pressure conditions and stresses.A numerical model for lithosphere dynamics has to be able to resolve these processes, and different theoretical descriptions of visco-elastoplastic rheology have been suggested (e.g., Andrés-Martínez et al., 2019;Gerya, 2019;Kaus, 2010;Lavier & Buck, 2002).In line with these studies we follow the approach by Moresi et al. (2003) and assume elastic, viscous and plastic (superscripts e, v, and p, resp.)deformation processes to be active simultaneously so that the respective strain rates are additive.This allows to adopt an additive decomposition of the total deviatoric strain rate, εij : where η v is the effective viscous creep viscosity, G is the elastic shear modulus, λ is the plastic multiplier and Q is the plastic flow potential.The stress-strain rate relation is based on the visco-elastic Maxwell-model: where η eff is an effective visco-elasto-plastic viscosity (discussed below), τij are the Jaumann-rotated old deviatoric stresses (for a detailed description see de Montserrat et al., 2019, and references therein), and Δt is the model time step.The second term on the right side represents the elastic deformation memory.
Viscous deformation includes diffusion (subscript dif) and dislocation creep (subscript dis) and we assume both creep mechanisms to be simultaneously active.The effective creep viscosity η v is thus given by All rocks are assumed to be isotropic, which allows use of a scalar measure for the strain rate magnitude in the dislocation creep law.Following Gerya (2019) we use the second invariant of the viscous strain rate tensor, εI , to describe the non-Newtonian behavior of dislocation creep.Diffusion and dislocation creep viscosities are defined as The pre-exponential factor A, grain size d, grain size exponent m, activation energy E and activation volume V depend on rock type and creep mechanism.The factors in front of A result from scaling triaxial and uniaxial experiment parameters (Gerya, 2019).Mantle peridotite forms the largest fraction of lithosphere at amagmatic ridges, and olivine is the most abundant and weakest mineral of this rock type.We therefore use a dry olivine rheology in our model (Hirth & Kohlstedt, 2003).B melt incorporates the effect of partial melts in intrusion regions on the rock's viscosity (see Section 2.4).
To include the plastic deformation into the viscous formulation we adopt the Prandtl-Reuss flow rule which, upon yielding, reduces the rock's stress state to the yield stress.The standard plasticity formulation is known to cause a mesh resolution dependent shear band width of ∼3 elements (Lavier et al., 2000).To reduce this unwanted numerical effect we implemented a viscoplastic regularization following Duretz et al. (2021).This allows the model to build up overstress to a certain level controlled by the viscoplastic viscosity η vp .We implemented a yield criterion τ yield similar to a Drucker-Prager criterion, however, we use the local lithostatic pressure P ls instead of the total pressure P from Equation 2. We prefer this so-called "depth-dependent von Mises" criterion over the Drucker-Prager criterion because the latter is known to potentially cause numerical convergence problems (Spiegelman et al., 2016).The resulting yield function is given by with C and Φ being the cohesive strength and the friction angle, respectively.The rock's stress state is quantified using the second invariant of the deviatoric stress tensor τ II .
In summary, the effective viscosity including the effect of plastic yielding can be expressed by Using Equations 8-11 with the parameters of the dry olivine rheology and a geothermal gradient of 0.05°C m 1 , we find that the transition from elastic-brittle to ductile behavior occurs at ∼750°C at a depth of ∼15 km.This is in agreement with the temperature assumed for the lower end of the seismogenic zone at (ultra-)slow spreading ridges (Grevemeyer et al., 2019) and we refer to the T bd = 750°C isotherm from here on as the base of the brittle lithosphere.For numerical stability reasons, viscosity is limited between 10 18 and 10 23 Pa s.

Strain Weakening
Strain weakening of fractured rocks is a critical process that needs to be considered in numerical models for lithosphere faulting.It parametrizes the weakening effects of pore fluids, gouge materials and mineral transformations in regions of strong plastic deformation, that is, shear zones.Strain weakening should thus scale with the deformation that a rock has experienced.It is incorporated in our model by reducing rock cohesion by ΔC C 0 = 0.9 over the range of accumulated plastic strain ɛ II,p from 0 to ɛ p,sw = 1.0: To calculate plastic strain accumulated in yielding regions, the plastic strain rate, which corresponds to the total strain rate minus the visco-elastic strain rate, is integrated over time, following

Serpentinization
To investigate the lithosphere-scale influence of the reaction of olivine and water to form serpentine at temperatures below 350°C, we adopt the temperature-dependent kinetic rate used by Rüpke and Hasenclever (2017): Geochemistry, Geophysics, Geosystems with S 0 = 10 13 s 1 (Rüpke & Hasenclever, 2017) and a = 808.3,b = 3,640 K, T 0 = 623.6K, and c = 8,759 K (Malvoisin et al., 2012).Note that we do not incorporate grain size effects on the kinetic rate (Malvoisin et al., 2012).This rate describes the relative change of olivine mass per rock volume ρ ′ ol , thus allows to calculate the serpentinization degree s of the rock by Furthermore, this representation directly allows to calculate the exothermic heat released during the reaction (Iyer et al., 2010) by: where Q serp is the heat produced per unit mass of olivine of 2.9 × 10 5 J kg 1 (value for forsterite from Iyer et al. ( 2010)).
Since we do not explicitly resolve the dynamics of hydrothermal fluid flow through a porous medium, we follow the approach of Bickert et al. (2020) to asses the availability of water for the reaction.The product of stress and accumulated strain, representing the volumetric work done to the rock by deformation, is used as a proxy for the formation of a fracture network giving water access to the rock.The threshold for this work, above which temperature-controlled serpentinization is activated, is set to 10 8 J m 3 .
Serpentine properties are listed in Table 1 and are used to derive effective rock properties by calculating the volumetric mean between olivine and serpentine properties using the serpentinization degree (cf.Equation 3).Mechanical weakening due to serpentinization is incorporated by reducing cohesion and friction angle (Escartín et al., 1997).The occurrence of serpentine in the oceanic spreading context is limited to shallow and cold regions.We therefore refrain from adjusting the viscous flow law parameters for serpentine, as serpentinized regions are expected to deform predominantly elasto-plastically.

Dynamic Grain Size Evolution
Grain size reduction occurs through dynamic recrystallization at sufficiently high applied stresses and results in a decrease of the diffusion creep viscosity (Equation 9).We adopt the dynamic grain size evolution model by Ruh et al. (2022), which is based on the palaeo-wattmeter (Austin & Evans, 2007): Equation 19 describes the change of grain size r under the influence of grain growth driven by surface energy (first term) and grain size reduction through dynamic recrystallization (second term).K g is the growth rate constant, f H 2 O is the water fugacity, E G and V G are activation energy and volume, respectively, and p is the growth exponent.γ is the grain boundary energy, the geometrical factor π results from the assumption of spherical grains.Ψ = τ ij εij,dis is the mechanical work rate from dislocation creep, of which the portion λ goes into grain size reduction.This reduces the heat from viscous dissipation, which takes the form of Geochemistry, Geophysics, Geosystems 10.1029/2023GC011331 Experimentally derived model parameters are from Speciale et al. (2020) and listed in Table 1.For numerical stability we limit the grain size from 1 μm to 10 mm consistent with observations (e.g., Bickert et al., 2021).

Hydrothermal Cooling
To incorporate the large-scale cooling effect of fluid circulation, we implemented an enhanced thermal conductivity.The enhancement factor is commonly referred to as the Nusselt number Nu in this context-not to confuse with its classical definition-as it up-scales conductive heat transfer to account for the missing cooling effect from the advection of hydrothermal fluids (e.g., Gregg et al., 2009;Phipps Morgan et al., 1987;Theissen-Krah et al., 2011).Here we adapt the approach of Gregg et al. (2009), in which thermal conductivity additionally depends on temperature and depth.We have modified this formulation by using lithostatic pressure instead of depth and applying a linear taper of the conductivity-scaling across the brittle-ductile transition between 600 and 750°C.The latter avoids the temperature dependence in the exponential term and reduces non-linear feedbacks.Thus, we obtain: with Nu = 1 + θ(Nu 0 1) exp (β(1 ) κ 0 is a rock type dependent reference conductivity, Nu 0 is a reference Nusselt number, θ is the temperaturedependent taper, β = 0.75 a smoothing factor, and P ht = 330 MPa the scaling lithostatic pressure for hydrothermal activity.For a flat topography, P ht corresponds to 10 km depth.Since our models always form an axial valley and lithostatic pressure laterally balances at relatively shallow depths, this value corresponds to approximately 6 km depth below the ridge axis, which is a typical assumption for the maximum depth of vigorous hydrothermal convection (e.g., Gregg et al., 2009).
Motivated by the approach of Lavier and Buck (2002), who enhance cooling around active faults, we have incorporated the factor K FZ into Equation 21b.We use this factor to account for enhanced hydrothermal fluid circulation in damaged regions, that is, fault zones, where highly permeable pathways are likely to form.K FZ is a function of the product W = τ II ⋅ ε II,p ⋅ εI I of the second invariants of stress τ, accumulated plastic strain ɛ p and strain rate ε.K FZ follows a logarithmic increase between the activation threshold of W 1 = 4 × 10 6 Js 1 m 3 and the upper limit of W 2 = 10 3 Js 1 m 3 : Analogously to the activation of serpentinization, K FZ parametrizes the formation of a connected fracture network (τ • ɛ p ) but additionally favors active fault zones ( ε) , where processes clogging the tectonically opened pore space, such as mineral precipitation from hydrothermal fluids, mineral transformations in the rock and consolidation of fractured rocks, have not yet closed the fracture network again.The activation threshold W 1 is derived from the serpentinization threshold of 10 8 Nm 3 and a typical fault zone strain rate of 4 × 10 14 s 1 , while the upper limit W 2 has been defined empirically to have maximum cooling in the center of the most active shear zone.F max is the exponent of maximum cooling intensity, which we vary in our parameter study to test the influence of hydrothermal fault zone cooling.

Magmatic Intrusions
Magmatic intrusions are periodically emplaced as sills of variable size and periodicity.Based on the strategy by Fan et al. (2021) and consistent with thermal models by Chen et al. (2022), emplacement is assumed to occur instantaneously below the shallowest point of the 1000°C isotherm, representing the basaltic solidus T sol .This corresponds to a mean depth of about 18 km, below seafloor, however, the exact positions of individual intrusions depend on the model dynamics.The emplacement temperature is a basaltic liquidus temperature T liq of 1200°C.Upon cooling, latent heat of crystallization is assumed to be evenly released between liquidus and solidus temperature: where ρ is the density, DT Dt is the total derivative of temperature with respect to time, and Q lat = 335 kJ kg 1 is the latent heat per kilogram of melt.Numerically this is implemented by increasing the specific heat capacity in intrusion regions between liquidus and solidus temperatures (C p in Equation 4).To account for the effect of partial melts on the viscosity of crystallizing intrusions, we adopt the exponential formulation by Hirth and Kohlstedt (2003) and multiply an additional factor B melt to the pre-exponential flow law parameters (A dis , A dif in Equations 9 and 10).For simplicity, we do not track melt content, but assume melt fraction ξ to be a linear function from 0 to 1 between solidus and liquidus temperature.This yields We use a moderate pre-exponential factor k = 35 so that the intrusion viscosity remains at the lower cut-off viscosity of 10 18 Pa s until intrusions have crystallized to about 25%.From then on, viscosity increases until complete crystallization.For simplicity we do not consider rheological differences between solidified basaltic intrusions and other lithospheric rocks so that any rheological effect of intrusions vanishes once they have cooled below the solidus temperature.Further mechanical aspects of intrusions such as volume changes and related effects on the stress field are not considered here as they require a compressible mechanical model (de Montserrat et al., 2019).

Code Structure
We use the 2-D Galerkin finite element method (FEM) with unstructured triangular elements to solve the system of equations in a Lagrangian reference frame (i.e., the mesh is advected according to the calculated deformation).All meshes are generated using the mesh generator Triangle by Shewchuk (2002).Mechanical and thermal parts of the model are solved sequentially over every time step.The mechanical solver uses Crouzeix-Raviart elements with continuous quadratic order shape functions to approximate the deformation field and linear, discontinuous shape functions to approximate the pressure.We employ the so-called penalty method (Donea & Huerta, 2003, and references therein) to derive an incompressible deformation field by iteratively solving Equations 1 and 2.
Next, we solve for heat diffusion and heat sources using a fully implicit FEM with quadratic order (6-node) triangles.The converged flow field is then converted to a displacement field using the current time step to advect the nodes of the finite element mesh.This step accounts for advective heat transport.Last, the second-order variables such as accumulated strain, serpentinization degree, grain size, lithostatic pressure and rock properties are calculated for the new configuration.
Between two time steps, mesh quality is checked to decide whether or not a remesh procedure is necessary.If so, a new mesh is created that preserves the high resolution regions near the ridge axis and that is adaptively refined along the active fault zones and in regions where intrusions are emplaced.To save computation time, these high resolution mesh regions are coarsened again once the intrusion temperature falls below the solidus temperature.After a remesh, all important variables are transferred from the old to the new mesh by interpolation using the finite element shape functions.The transfer of variables only stored at integration points requires an intermediate

Geochemistry, Geophysics, Geosystems
10.1029/2023GC011331 mapping step from integration points to nodes before interpolating to the location of integration points in the new mesh (see Figure 3 in de Montserrat et al. ( 2019)).
The time step size dynamically adjusts to the model dynamics, for example, flow field and temperature changes, between 0.5 and 5 Kyr.Upon intrusion emplacement, the thermal solver sub-iterates with time steps down to 100 years in order to properly resolve the heat transfer from the intrusion to its surroundings where high temperature gradients exist.The model is written in MATLAB (ver.R2020a, www.mathworks.com) and uses a vectorized element assembly for an improved performance as suggested by Dabrowski et al. (2008).In addition we make use of the libraries SuiteSparse (Davis & Hager, 2009) and MUTILS (Krotkiewski & Dabrowski, 2010).Data visualizations are done using MATLAB and Tecpot 360 EX (ver.2021 R1, www.tecplot.com).

Model Setup
The model setup is illustrated in Figure 1.The model domain is 200 × 80 km, discretized with a triangular finite element mesh and variable mesh resolution.The mean node spacing in the region of the axial lithosphere is 375 m and increases to up to 1,250 m toward the bottom boundary.Additionally, the mesh is adaptively refined along actively deforming regions during each remesh procedure with a node spacing of 250 m in the fault zones.The model has a free surface (Andrés-Martínez et al., 2015) with z = 0 km corresponding to a water depth of 4,650 m and constant mantle inflow along the bottom boundary that balances the lateral extension of 7 mm yr 1 on both sides.Spreading rate and water depth represent the amagmatic segments of the SWIR at 62°E-65°E (Cannat et al., 2006(Cannat et al., , 2008)).To stabilize the initial model phase, we trigger the first fault by imposing a linear profile of accumulated strain (Figure 1).
The initial temperature field is based on Bickert et al. (2020), however, we assume a slightly thinner axial brittle lithosphere of 15 km in agreement with recently reported data (Chen et al., 2023).The geometry of the initial temperature field has been adjusted accordingly, see Figure 1.At the domain top and bottom, isothermal boundary conditions are set to 4 and 1,400°C respectively.The lateral domain boundaries are insulating.
For our parameter study, we use one configuration without magmatic input and five configurations with intrusion geometries and periodicities ranging from 0.14 × 1.4 km sized sills every 40 Kyr to 0.5 × 5 km sized sills every 100 Kyr (Fan et al., 2021).The ratio of intrusion volume flux (average over time) to brittle lithosphere flux (15 km × 14 km/Myr) serves as a proxy for the fraction of lithosphere extension that is potentially compensated by magma input through dikes, commonly referred to as the M-factor (Buck et al., 2005).Assuming that the melt budget provided by the intrusions could go into diking, the above intrusion scenarios correspond to equivalent hypothetical values of M between 0.02 and 0.12.Note, however, that spreading in our models is always amagmatic and the actual M-factor in all calculations is zero because neither the volume flux of intrusions nor diking are considered.To investigate the effect of hydrothermal cooling, we vary the magnitude of enhanced hydrothermal fault zone cooling using six values of F max ranging from 0 to 2.5.In each setup, Nu 0 in Equation 21b is adjusted to F max and M and take values between 6.4 and 10.2 to obtain a similar average axial brittle lithosphere thickness for all setups.The exact model parameters for each setup can be found in Figure S1 and Table S1 of the Supporting Information S1.For all experiments, we set the maximum simulation time to 40 Myr.

Results
Lithosphere deformation is controlled by several feedback mechanisms between thermal and rheological subprocesses.In our model calculations we observe various faulting modes that are illustrated in   2021)); vertical exaggeration 2:1; B-Breakaway, E-Emergence, #1 is the currently active detachment, C5-magnetic anomaly #5.Modified after Cannat, Sauter, et al. (2019).
Geochemistry, Geophysics, Geosystems 10.1029/2023GC011331 In Section 3.2 we present the fully coupled models with a dynamically evolving temperature field.We first show one calculation without any hydrothermal cooling (Section 3.2.1) to emphasize the importance of this additional cooling mechanism for reproducing the observed lithosphere thickness.Afterward we present the different faulting modes observed in the parameter study.

Model With Imposed Temperature Structure
Tests with this simplified thermal model show particular sensitivity to weakening parameters and the prescribed thermal structure, similar to the findings of Bickert et al. (2020).We observe a very stable and uniform sequence of flip-flop faults (Figures 3 and 4).The shown faulting sequence starts at 11.1 Myr, after nine flip-flop detachment faults have already formed in the same manner.Each detachment starts as a high angle normal fault (Figure 3a).With increasing vertical displacement, the footwall rolls back under its own weight.Thereby it experiences large flexure, indicated by the formation of a region with increased strain rate pointing from the central footwall to the emergence of detachment D10 (Figure 3b), consistent with the dynamics suggested by Reston and McDermott (2011).Since the temperature field evolves with the seafloor relief, the upward movement of the footwall results in a shallow temperature maximum in the central footwall close to the center of the incipient solid-block rotation.The resulting rheological weakening focuses flexural strain (Figure 3c) and acts as a trigger for the new fault D11 with a polarity opposite to its predecessor D10, in agreement with the assumption that a new fault cuts the hinge of its predecessors footwall (Reston, 2018).Eventually, a shear-localizing feedback loop between deformation and grain size reduction/serpentinization is activated.This leads to the formation of a new fault accounting for a major portion of extension and the cycle starts over again (D11 in Figure 3d).
The stability of the flip-flop cycles are shown in Figure 4. We find that successive strain rate increase and focusing occur simultaneously with grain size reduction prior to the "activation" of the new fault zone, which from then on accommodates a major part of the deformation (dashed vertical lines in the left panel of Figure 4).For evaluating each fault's lifetime we take the time span between the activation of one detachment to the next.Thus, the mean fault duration is also obtained by dividing the total duration of a faulting sequence by the number of faults, allowing for comparison with the SWIR data by Cannat, Sauter, et al. (2019).The results of this simulation are in perfect agreement with the data from the SWIR.
The strain rate of each fault peaks at around 3 × 10 13 s 1 , while fault zone width and minimum grain size converge to 1 km and 1.5 μm, respectively, in agreement with seismic and petrological data (Bickert et al., 2020(Bickert et al., , 2021;;Momoh et al., 2020).
Besides the similarities between our simulation, SWIR observations and the classical flip-flop mechanism (Reston, 2018;Sauter et al., 2013), we also see some differences.From the theoretical model, we would expect the new fault to cut its predecessors exhumed shear plane a few kilometers away from its emergence.Also, with the formation of the new fault, the predecessor should successively fade out as it reaches the colder and thicker off-axis lithosphere.Instead, new faults in our simulations often cut the shear plane around the emergence, and the preceding fault remains active and continues to exhume mantle material to the seafloor.The two simultaneously active faults enclose a central horst with predominantly vertical motion, until it tilts toward the older fault and initiates the mechanism discussed above.Slip on the old detachment ends with the formation of the second next fault with same polarity (D9 becomes inactive when D11 forms; see Figure 3d).Interestingly, this faulting sequence can also explain the segments of abnormal seafloor age pattern at the SWIR (Cannat, Sauter, et al., 2019).To observe more classical flip-flop detachments with more pronounced roll-back, a stronger reorientation of the fault zone from initially straight to curved would be required.Tests have shown that a lower yield stress, that is, lower friction angle and/or cohesion, or an increase of η vp in the visco-plastic regularization lead to faster yielding and wider fault zones, both allowing for more efficient fault zone reorientation and roll-back.
Our model with imposed temperature structure generally confirms the conclusion of Bickert et al. (2020)despite the differences in modeling techniques and parameterizations-that grain size reduction allows for strain focusing in the deep lithosphere at the root of the incipient detachment fault, which in turn facilitates flip-flop faulting.An important prerequisite for triggering a new fault of opposite polarity is, however, a higher temperature in the footwall relative to the colder active fault zone(s) (cf.Figures 2a and 3b).In the next sections we investigate how such a thermal structure can dynamically evolve through the interplay of tectonic, magmatic and hydrothermal processes.

Models With Dynamic Thermal Evolution
We conducted 36 simulations with a dynamically evolving temperature field to test the effects of hydrothermal cooling-widespread plus enhanced to a certain level within and around active fault zones (see Equation 21b)and magmatic input in form of sill intrusions that can be converted into a hypothetical M-factor.In all model simulations, the mean axial brittle lithosphere thickness, defined as the minimum depth of the 750°C isotherm below the seafloor, is in the range of 15.1 ± 1.0 km, which matches observations and facilitates a good comparison of the different model calculations.These simulations are labeled using an index of 1-36 to refer to them later in the summarizing parameter space plot (Figure 11).Additionally, we show simulation 0 in the next section, which has no magmatic input and no hydrothermal cooling at all.Over the following sections we will present the different faulting modes observed in our simulations.

Miniature Flip-Flop (No Hydrothermal Cooling)
Without the additional cooling effect of hydrothermal circulation (Nu = 1), the upwelling hot mantle reduces brittle lithosphere thickness (defined by the 750°C isotherm) to less than 4 km (Figure 5).This result is incompatible with micro-seismicity data indicating a thickness of around 15 km (Chen et al., 2023;Grevemeyer et al., 2019).Also, the seafloor topography is much smoother than observed at the SWIR.We show this calculation to emphasize the importance of hydrothermal cooling not only at predominantly magmatic fast-and intermediate spreading ridges but also at amagmatic slow-and ultraslow ridges.Interestingly, the faulting pattern is very similar to that of the model with imposed temperature structure and could be considered a miniature version of flip-flop faulting.However, such a thin lithosphere is typical for much faster spreading ridges or for sections of slow spreading ridges with a significantly increased melt budget.We therefore consider this setup to be unrealistic for an amagmatic slow spreading ridge.

Long-Lived Detachment Mode
Each of the simulations of the parameter study shows several faulting modes (cf. Figure 2).There is no calculation that maintains a single faulting mode throughout the 40 Myr simulation time.However, certain faulting modes clearly dominate depending on the intensity of fault zone cooling and magmatic input.Note that all calculations include a wide-spread background hydrothermal cooling to achieve a realistic lithosphere thickness.
We first focus on simulation 1 (Figures 6a and 6b), which has neither enhanced cooling of fault zones nor magmatic input (M = 0, F max = 0, i.e.K FZ = 1 in Equation 21b).With increasing displacement on the fault plane, hot mantle material is pulled up with the footwall.The advective velocity and thus the pull-up is largest close to the active fault zone resulting in a temperature maximum at this position (Figure 6a)-similar to the findings of Behn and Ito (2008).With ongoing fault activity, the temperature maximum then moves together with the fault zone toward the hanging wall side (Figure 6b).Thereby, the fault remains in a hot, weak environment instead of migrating into thicker, stronger lithosphere.This generally allows the fault to stay active significantly longer than a typical flipflop detachment and leads to a long-lived detachment that may be accompanied by accommodating footwall faults.
For most simulations with low hydrothermal fault zone cooling and low magmatic input, we observe a behavior similar to that of simulation 1. Long-lived detachment faults and their accommodating faults are the dominant faulting mode.We observe three variations of long-lived detachments: (a) rolling-hinge detachments, (b) detachments with synthetic accommodating footwall faults, and (c) detachments with antithetic accommodating footwall faults.The different modes are illustrated by sketches in Figure 2a and by snapshots from two different simulations in Figure 6.
Rolling-hinge detachments (Figures 6a and 6b) are characterized by a rotation of the footwall that transitions into a horizontal plate motion further off-axis (Buck, 1988;Lavier et al., 1999).This goes along with stretching and shearing in the footwall, but the emerging shear bands do not become focused enough to form a new fault zone (Figure 6a).At the seafloor, a very smooth and kilometer-long shear plane is exhumed.
Synthetic accommodating faults are secondary faults of the same orientation as the detachment but opposite curvature in the upper part (Figures 6c and 6d).Similar faults have been observed by Lavier et al. (2000) and Sandiford et al. (2021).In-between this fault pair, a closed solid-block rotation with very little shear deformation  Geochemistry, Geophysics, Geosystems 10.1029/2023GC011331 begins.The accommodating fault either fades out after some time or realigns its curvature and may become the major active detachment fault itself.
Antithetic accommodating faults (Figures 6e and 6f) cut the footwall of the active detachment with opposite polarity.The general deformation in the footwall remains dominated by the curved long-lived detachment.After producing significant relief, but before experiencing major roll-back, the accommodating fault fades out instead of becoming a full-grown detachment.We find that one long-lived detachment fault can cause and survive multiple antithetic accommodating faults.Interestingly, the resulting fault zone pattern and seafloor topography closely resembles that of flip-flop faulting.

Detachments With Opposite Polarity
One characteristic of flip-flop faulting is the alternating polarity of successive detachments, each of which accommodates a significant amount of extension.A sequence of such detachments with alternating polarity is thus required for considering a simulation to exhibit flip-flop faulting.
Figures 7 and 8 show a temporal evolution from simulation 28 (M = 0.1, F max = 1.5), which is an example for such a sequence and shares many characteristics discussed above for the model with imposed temperature structure.With increasing displacement on the detachment fault, the footwall experiences roll-back (Figures 7a and 7b left column).At the same time, the increased cooling of the active fault zone(s) shifts the isotherms around the detachment to slightly greater depth (Figures 7a and 7b right column).This cooling effect increases the temperature difference between active fault zone and footwall and thereby pushes the temperature maximum further into the footwall.The emplacement of intrusions (size: 4,000 × 400 m every 80 Kyr) near the temperature maximum further weakens the footwall, which eventually allows for enough focusing of flexural strain to initiate the feedback loop between grain size reduction and strain rate increase (Figures 7b and 7c).After the new detachment fault of opposite polarity has formed (Figure 7c), heat is extracted from its root by (in our case parametrized) hydrothermal circulation.The heat extraction causes the temperature maximum again to move into the footwall of the new detachment (Figure 7d) and the next cycle begins.
In the model with a dynamic thermal evolution, footwall roll-back is more pronounced and takes longer compared to that with an imposed temperature structure.This is reflected in a smoother seafloor relief and longer fault life times.The stronger roll-back is possibly related to slightly wider fault zones.Compared to the model with imposed thermal structure (Figure 4), the strain rate of individual faults (Figure 8) does not drop as significantly at the activation of the next detachment, indicating a different partitioning of the total deformation between the active detachments.Furthermore, within the dynamic thermal model flip-flop sequences are shorter and faulting is generally more diverse including cross-cutting faults and more variable fault dips and life times.
While single detachments of opposite polarity are initiated at lower hydrothermal fault zone cooling, it requires F max ≥ 1.5 to achieve stable sequences of alternating polarity in the absence of intrusions (simulation 4).A similar trend is observed in simulations without enhanced cooling of fault zones, where a magmatic input equivalent to M = 0.095 (simulation 25) or more is required to achieve such stable sequences.Combinations of both, moderately enhanced fault zone cooling and moderate magmatic input, also reveal stable sequences of alternating polarity.This indicates that a combination of hydrothermal activity in fault zones and magmatic sill intrusions below the footwall may favor flip-flop faulting.

Deep-Cutting Faults: Criss-Cross and Spider Modes
Another type of faults that we observe in many simulations are deep-cutting faults.The name refers to the case where a new fault of opposite polarity cuts its predecessor several kilometers below the seafloor.The two faults then form an X-shape that encloses a large central block with significant surface relief (Figure 9).After the formation of such a deep-cutting fault pair, we observe two different faulting pattern evolutions.
The first mode is what Bickert et al. (2020) have named spider mode, where both deep-cutting faults remain active almost symmetrically.For geometrical reasons, this configuration requires an alternating activity of the faults below the crossing point, which temporarily induces strain and weak regions in both footwalls.These regions of accumulated strains migrate upwards and form the "spider legs" faintly visible in the strain rate plots (Figures 9b  and 9c).While magmatic intrusions are not necessarily required to activate spider mode, they support and stabilize this mode by weakening the region below the intersection of the two faults.Spider mode typically ends after some million years when one fault becomes more dominant due to an emerging asymmetry.
Criss-cross mode is a new mode identified in our simulations.It represents the second option for the evolution of deep-cutting faults, in which the two faults migrate through the central block by slicing it in several discrete steps.
Each step is a short-lived normal fault that offsets a part of the central block (Figure 9e).The name criss-cross refers to the alternating, almost perpendicular normal faults that cut-off sections of the central block.The number of dissecting normal faults during the criss-cross phase varies, and so does the impact on the enclosed surface relief.Criss-cross mode ends when the intersection point of the two deep-cutting faults reaches the seafloor allowing them to separate and migrate off-axis in opposite directions, similar to "detachments of opposite polarity".
We find criss-cross faults to become more frequent with stronger fault zone cooling.At increased F max (Equation 21b), new faults tend to form at a higher frequency.Thus, there is less time for shifting the temperature maximum away from the root of a fault into its footwall, and the next fault is more likely to form as a deep-cutting fault close to the root of its predecessor.The more spontaneous fault initiation also facilitates the formation of the sequence of dissecting normal faults.

Chaotic Faulting Patterns
At the largest values of F max that we tested, fault formation becomes mostly chaotic.New faults form at high frequency, can be scattered around the axial region and possess variable fault dip and curvature.In this mode, more than two faults may be active simultaneously and accommodate comparable portions of total extension.Since each fault is repeatedly dissected by new faults, the surface relief is continuously reworked, which renders an interpretation of seafloor topography rather complicated.An exemplary sequence from simulation 24 (M = 0.07, F max = 2.5) displaying this chaotic mode is shown in Figure 10.Simulations featuring this mode do not resemble any observations from amagmatic ridges and we therefore consider this setup to be unrealistic.

Faulting Mode Summary
The faulting modes resulting from varying fault zone cooling and magmatic input are summarized in Figure 11 and we observe several trends.The occurrence of long-lived detachments and their accommodating faults is mainly limited to simulations with relatively low fault zone cooling (F max ≤ 1), while criss-cross mode evolving from deep-cutting faults mainly occurs for F max ≥ 1. Chaotic mode is limited to the largest values of F max ≥ 2.0.Spider mode is observed almost equally over the whole range of F max that we tested, however it occurs more frequently with increasing M. Consequently, the number of observed long-lived detachments and criss-cross faults slightly decreases with increasing M. Detachments of opposite polarity are the most abundant faulting mode in our simulations with F max ≤ 2.0.
Over the next sections, we will discuss the mechanisms that produce flip-flop detachment faulting in our model and evaluate different fault sequences from our results in the scope of observations from the SWIR.

Thermal Structure of the Footwall
Our simulations demonstrate that a crucial prerequisite for sequences of alternating detachments that constitute flip-flop faulting is a temperature maximum in the central footwall.In the model with an imposed temperature structure this automatically evolves as the footwall is uplifted.In the fully coupled model with a dynamic thermal evolution, footwall uplift together with hydrothermal fault zone cooling yields a similar thermal structure that is reinforced by magmatic intrusions.To assess the influence of fault zone cooling and magmatic sill intrusions on the temperature field and thus faulting pattern more systematically, we reduce the number of non-linear feedbacks in the fully coupled model.This is done by taking the deformation field from simulation 28 and combining it with different parameter combinations of F max and M to compute alternative temperature evolutions.In Figure 12 we compare the temperature field of simulation 28 with those evolving with the same magmatic input but no enhanced fault zone cooling (setup of simulation 25), same fault zone cooling but zero magmatic input (cf.simulation 4) and without either of the two processes being active (cf.simulation 1).We will analyze the temperature changes relative to the temperature field evolving without intrusions and enhanced fault zone cooling.Figure 12 shows how the parametrized hydrothermal activity efficiently cools the currently active D7 and, to a lesser extent, the preceding fault zone D6 (Figure 12a).This shifts the temperature maximum at the brittle-ductile transition away from the active fault zone up 3 km into the central footwall (Figure 12b).The emplacement of intrusions near the temperature maximum reinforces (Figure 12a) and narrows this high-temperature region by about 6 km (Figure 12c).Note that the cooler regions below the active faults result from stronger background cooling of the brittle lithosphere (i.e., larger Nu 0 ), which is adjusted to the extra heat input from the intrusions to reproduce the observed lithosphere thickness.
The effects of these temperature variations on the faulting mode can be best understood by translating the temperature field in the deep brittle lithosphere and below to a modified viscosity and thus strength structure.Reduced temperature in the deeper fault zone corresponds to a viscosity increase, making the fault zone less efficient after a period of slip and cooling.Magmatic intrusions on the other side cause increased temperatures in the footwall below the brittle lithosphere, serving as a weak spot to focus flexural deformation.This combination of reducing the active fault zone efficiency and providing a weak spot in a geometrically favorable position is what triggers the next detachment of opposite polarity.The repetition of this cycle leads to the flip-flop faulting sequence (Reston, 2018;Sauter et al., 2013).

Hydrothermal Cooling of Active Fault Zones
As a first step toward a hydro-thermo-mechanical model that resolves porous flow within a visco-elasto-plastic deforming lithosphere, we have parametrized the heat transport by hydrothermal fluids taking into account the deformation state of the rock.To relate the magnitude of hydrothermal cooling predicted by our model to natural systems, we convert the enhanced thermal conductivity and resulting temperature gradients to an equivalent hydrothermal power output at a hypothetical fault-related vent site.To facilitate a comparison between all model calculations shown in Figure 11, we aim for a relation that only depends on the factor F max , by which we enhance the cooling effect of fault zones.Doing so allows to place a second axis of hypothetical hydrothermal heat fluxes on the right side of Figure 11.From our simulations we evaluate representative values for a fault-related seafloor temperature gradient of 0.05°C m 1 , and a fault zone width, over which an enhanced heat flux from the lithosphere into the ocean occurs, of 2 km.As fault zone thermal conductivity we take the conductivity of serpentinized mantle at the seafloor reduced by the background cooling.Using Equations 21b and 22 we obtain a fault zone thermal conductivity of κ FZ = κ 0,serp (Nu 0 1) exp (β(1 The value of Nu 0 for each value of F max is taken from the simulation at M = 0. Seafloor hydrostatic pressure P sf is evaluated in the axial valley.Multiplication of temperature gradient, fault zone width, and fault zone conductivity κ FZ provides an estimate for the total power output of the hypothetical vent fields.Although the resulting estimates should be viewed with caution because of the uncertainties involved, they allow at least for an order-of-magnitude comparison with vent field observations.Fault related heat fluxes predicted by our model (axis on the right-hand-side of Figure 11) span a plausible range from the low-temperature Lost City Hydrothermal Field located on a core complex off-axis the Mid-Atlantic Ridge (Lowell, 2017) to the hightemperature Longqi-1 Vent Field located on the hanging wall of a detachment fault in a magmatic section of the SWIR (Tao et al., 2020).
The hydrothermal site that is most relevant for the flip-flop dominated section of the SWIR investigated here is the recently discovered Old City hydrothermal field (Cannat, Agrinier, et al., 2019) located on the exposed fault surface of the currently active detachment in the amagmatic segment.Total heat flux estimates or venting temperature data are still missing, however, the vent field has a larger spatial extent than Lost City (Cannat, Agrinier, et al., 2019) and an estimated comparable venting temperature of <100°C (Fan et al., 2021).We therefore suppose that values of F max much larger than 1.5 seem unlikely, because they would correspond to heat fluxes of more than 20 times that of Lost City.However, venting temperatures of >335°C of the inactive Tianzuo hydrothermal field inferred from geochemical data (Ding et al., 2021) could point toward higher heat fluxes.More data on those and other, presumably existing (Baker et al., 2004) yet undiscovered vent fields in similar tectonic settings, is needed to constrain representative numerical models.Note that our model cannot account for the limited life time of individual vent fields (Humphris & Cann, 2000) but parametrizes the time-averaged cooling effect of all hydrothermal activity.
An open question is, to which depth hydrothermal circulation can efficiently cool the lithosphere.The maximum depth of hydrothermal fluid flow is mainly controlled by the subsurface permeability structure, which is a difficult to assess quantity and suggested values for oceanic lithosphere span orders of magnitude (e.g., Andersen et al., 2015;Barreyre et al., 2018;Hasenclever et al., 2014;Marjanović et al., 2019).Models and data generally suggest a decrease of permeability with depth (e.g., Kuang & Jiao, 2014, and references therein), which would inhibit hydrothermal heat extraction at depths greater than a few kilometers.At temperatures above 600-800°C, pore space is assumed to be thermally closed to vigorous fluid flow (Hirth et al., 1998).Fan et al. (2021) argue that under certain circumstances thermal cracking can extend the reach of hydrothermal convection down to 10 km depth.Tao et al. (2020) model fluid circulation in the brittle lithosphere over a thickness of 13 km if permeability increase in the fractured fault zone is sufficient.Observation of hydrothermally altered minerals at this and comparable tectonic settings indicate that at least a limited amount of fluids reaches down to the brittle-ductile transition, even in the presence of an extremely thick lithosphere (Bickert et al., 2023;Patterson et al., 2021;Vieira Duarte et al., 2020).However, these observations may represent the limit of maximum fluid penetration depth, while active hydrothermal convection and cooling may be efficient only at shallower depth.Our model accounts for the discussed effects by having a pressure-dependent scaling of thermal conductivity leading to a decrease with depth in agreement with previous studies (e.g., Gregg et al., 2009), however possibly underestimated compared to scaling lengths given by Kuang and Jiao (2014), a fading influence across the brittleductile transition and a significant increase of the hydrothermal cooling effect around actively deforming fault zones.Note that the cooling effect incorporated by the Nusselt number parametrization spatially integrates the effect of fluid flow and may thus reach deeper than active hydrothermal circulation.To answer more in-depth questions about the role of active fault zones on hydrothermal circulation pathways-including potential feedbacks with magmatism and faulting-will require to resolve porous flow of hydrothermal fluids in the model instead of parametrizing its cooling effect.
While certainly affected by our model limitations, including its 2-D nature (discussed below), our results nonetheless indicate that enhanced cooling of active fault zones is a key process controlling the tectonic mode of a magma-poor ridge section.Simulations in which the tectonic state has no influence on the parametrized hydrothermal cooling (F max = 0) produce only limited flip-flop faulting, while simulations with a too strong enhancement of fault zone cooling do not produce flip-flop faulting at all.

Magmatic Sill Intrusions-Model Versus Data
In this study we have varied sill sizes as well as emplacement frequencies such that the magmatic input fraction corresponds to a reasonable range of M-values lower than 0.12.This can be considered almost amagmatic, and limited amounts of melt are in agreement with observations (Paquet et al., 2016;Sauter et al., 2013).Intrusions in our model facilitate flip-flop faulting through direct (rheological) and indirect (thermal) weakening.We hypothesize that further mechanical aspects which we have not considered in our model, such as additional stresses resulting from volumetric effects of intrusion emplacement and cooling or melt migration along faults could increase the influence of intrusions and affect the faulting types summarized in Figure 11.Moreover, we cannot speculate on the dominant faulting modes at greater M-values, where volumetric effects become critical and diking is essential to explain the transition to the more magmatic faulting mode of corrugated detachments.
When designing our model our intention was to put as few constrains on the dynamic evolution as possible.We therefore decided to not prescribe the location of intrusions but to choose realistic criteria for where ascending melts might gather to form a sill.We assume that melts migrate up to the shallowest position of the basaltic solidus temperature, which most often is located in the footwall below the active fault.While variable emplacement depths could affect our results, our chosen method is consistent with the predictions of the numerical study by Chen et al. (2022) on the thermal regime of magma-poor segments.After emplacement, cooling intrusions are transported upward where they are cut and sheared by subsequent faults.This cycle is in agreement with recent seismic data (Momoh et al., 2020): Stacked sills in the deep lithosphere can potentially explain a cluster of subhorizontal seismic reflectors between 9.3 and 13 km, while sheared intrusions could intensify the seismic signature of the detachment faults, which are identified by a series of dipping reflectors.

3-D Effects of the Ridge Axis
Both amagmatic corridors of the 62°-65°E segment of the SWIR are enclosed by more magmatic sections characterized by axial volcanism and symmetric normal faulting.This segmentation presumably results from strong along-axis focusing of melts (Cannat et al., 2003).The more magmatic sections do most likely serve as thermal and mechanical anchors for the evolution of their magma-poor neighbor-segments.We presume this to be the main factor missing in our 2-D model, which cannot account for along-ridge variations in magma supply or faulting patterns.Coupling between an amagmatic segment sandwiched between magmatic sections could prevent long-lived detachment faults from shifting tectonic activity several tens of kilometers from the initial position at the ridge axis.Fault life times would consequently be shorter and faulting patterns more stable.
However, being able to reproduce a stable flip-flop faulting mode in a dynamic thermo-mechanical 2-D model without prescribing the ridge axis position lets us speculate that the interplay of hydrothermal cooling (enhanced in fault zones) and the heat input plus weakening effect of sill intrusions are key mechanisms that facilitate flipflop faulting.

Tectonic Significance of Observed Faulting Modes
We performed a three-step analysis on each of the simulations.First, we picked the onset time of every fault zone and categorized it (Figure 11).In the second step, we identified fault sequences of at least 4 consecutive faults that fall in one of the two categories "apparent flip-flop" and "true flip-flop" (Figures 13a and 13b).As "apparent flipflop" we categorize faults that produce off-axis dipping fault zones and relatively regular spaced ridges, since these are the main features used to identify flip-flop detachments from observational data.For the category "true flip-flop" we additionally evaluate if the faulting mechanism resembles the flip-flop mechanism described in the previous sections.Faulting sequences are mostly terminated by a long-lived detachment or a deep-cutting fault/ spider mode phase.As a last step of our analyses, we calculate the mean fault duration during the faulting sequences (Figure 13c).

Fault Life Time
Before we discuss the relevance of the different fault sequences, we will first have a look at the observed fault life times during these sequences.Figure 13c shows that the mean fault duration of simulations with F max ≤ 1.5, including the presented flip-flop sequence from simulation 28, is about 1.6 Myr thus ∼0.5 Myr longer than estimated for the SWIR.This is presumably due to the missing 3-D effects of the ridge axis providing external triggers for the formation of new on-axis faults.Another potential reason is the chosen rheology with relatively high yield stress promoting well defined and stable fault zones.
Fault duration slightly increases with an increase of M. Before fault zone cooling can significantly shift the temperature maximum and thus intrusion emplacement into the footwall, intrusions are emplaced close to the root of the new fault zone, stabilizing and prolonging its initial phase.This mechanism is evident in Figure 12b, where in the first phase of D7 the temperature maximum is shifted toward the fault zone instead of into the footwall, if intrusions are considered (red and yellow lines).
Increasing F max to 2.0 results in a clear drop in the mean flip-flop fault duration.Simulation 5 for example, displays sequences of rather uniform flip-flop faulting with a mean fault duration of 1.0 Myr (see Movie S8).However, we cannot verify the role of a footwall temperature maximum as a trigger for the new faults in these simulations, as it does not coincide as clearly with the evolving shear band that develops into the next detachment of opposite polarity.We rather believe that instead the feedback loop between rheology, deformation and cooling that allows for a more spontaneous fault initiation explains the increased fault frequency.These feedbacks, however, occur very localized in the fault zone and to be able to properly address it, we would require to resolve hydrothermal flow inside the fault zone and to include effects of fluid-rock-interaction on rheology.Consequently, we see F max = 2.0 as the point, where the concept of our parametrization breaks down and we therefore hatched the corresponding region in the parameter space plots in Figure 13.

True Flip-Flop Faulting
From Figure 13b it can be inferred that the activation of flip-flop faulting mechanism is not simply caused by the absence of melt input but seems to result from a sensitive interplay of thermal, tectonic and magmatic processes.
The relevant parameter range that best reproduces the conditions necessary for a stable flip-flop faulting mode in line with the commonly accepted interpretation of SWIR observations is identified by the green band in Figure 13b.
A feature that we observe in many flip-flop producing simulations are antithetic faults in the hanging wall of the active detachment (Figures 7a, 7b, and 7d).These are highly consistent with microseismicity data (Chen et al., 2023).In our simulations, they often help to focus the footwall deformation resulting in the next fault of opposite polarity (Figures 7b and 7c).This could indicate that new faults indeed form around the emergence of their predecessor as predicted by our simulations.While we know of no evidence for still ongoing deformation along the preceding fault zone, that is, D2 in Figure 2b, we could imagine that activity alternates between the two faults on a time scale not resolved by existing observations instead of both being continuously active.
Flip-flop faulting has originally been proposed for the magma-poor rifted west Iberia margin (Reston & McDermott, 2011) and is also observed at another segment of the SWIR (51.5-53.5°E;Liu et al., 2020).However, at virtually amagmatic segments along other ultraslow-spreading ridges like the Gakkel or the Knipovich Ridge this distinct faulting pattern has not been observed yet (Meier et al., 2021;Michael et al., 2003).If flip-flop faulting is indeed absent at these settings, this mode may represent a special case that depends on the specific regional conditions-for example, thermal structure, stress field or the characteristics of neighboring magmatic segments-rather than being the normal state for ridges at this spreading rate and melt supply.This is further supported by the high sensitivity of stable flip-flop faulting to the imposed temperature field and rheology that we observed in tests for the simulation presented in Section 3.1.It could explain, why smooth seafloor areas are not observed more frequently.
Note that the parameter range yielding realistic results possibly goes beyond the simulations displaying the longest and most uniform flip-flop sequences, which is supported by both models and data.Reston and McDermott (2011) suggest alternating polarities to be dominant, but do not preclude successive detachments of the same polarity.Flip-flop sequences observed at mid-ocean ridges are of limited duration: eastern corridor of the 62°-65°E segment: 11 Myr (Sauter et al., 2013); 51.5°-53.5°Esegment: ∼7 Myr (Figure 7 in Liu et al. (2020)).They also show considerable variability: The western magma-poor corridor of the 62°-65°E segment of the SWIR appears to have a lower magmatic budget and more symmetric ridges than the eastern corridor, indicating at least a variable periodicity of the faulting sequence (Sauter et al., 2013).More symmetrical ridges in our simulations often reflect horst-like flip-flop detachments (especially in simulations with imposed thermal structure, Figure 3).Ridges with steeper outward-than inward-facing slopes as observed in the eastern corridor, are indicative of more intensive footwall rollback, or of another faulting mode (see Section 4.3.3).For a more indepth analyses of the surface relief, surface processes such as erosion and mass wasting would need to be considered.Without, our experiments generally produce stronger relief than observed.

"Apparent Flip-Flop" Modes
We now consider sequences featuring other modes besides true flip-flop faulting that are also compatible with some of the most important observations.These modes are antithetic accommodating faults in the wake of longlived detachments (Figures 6e and 6f) and criss-cross faults succeeding a deep-cutting fault.The latter ones are only classified as apparent flip-flop if they offset the enclosed ridge only slightly before they separate at the surface and continue just as regular flip-flop detachments (e.g., Figures 9d-9f).There are some arguments in favor of considering these modes as potential alternative interpretation of the SWIR observations, especially since they occur over a wide and realistic parameter range as Figure 13a shows.Both modes produce a series of off-axis, outward dipping inactive faults and pronounced, more or less regularly spaced ridges.Furthermore, they are generally able to produce abnormal seafloor age pattern (see arrows indicating seafloor age in Figure 2): antithetic accommodating faults as they are passively transported toward their hanging wall side after becoming inactive, and deep-cutting faults like regular flip-flop faults, depending on where they finally cut their predecessors exhumed shear plane.However, this is mainly a conceptual notion, because seafloor age with distance from the SWIR axis is, to our knowledge, mostly inferred from sampling breakaway and fault emergence locations and interpreting these within the conceptual kinematics of flip-flop faulting.
Even though high-resolution magnetic profiles exist for this region (Bronner et al., 2014), only one magnetic anomaly can be confidently traced to the amagmatic corridors (C5 in Figure 3d) due to the extensive serpentinization and very low contribution of magmatic rocks.These data indicate long-term symmetric spreading compatible with flip-flop mode as well as antithetic accommodating (cf. a L-r-r-R-l-l sequence, where uppercase letters stand for long-lived detachments and lowercase letters for accommodating faults) and deep-cutting faulting modes.We are not aware of other data such as high-resolution petrological profiles that would rule out one or the other faulting mode.Moreover, some of the interpreted fault zones are still under debate.For example, D3' interpreted by Corbalán et al. (2021) from seismic velocities has significantly smaller relief than the other detachments and could well be attributed to one of our alternative two modes.Also, both modes can produce asymmetric ridges in agreement with the bathymetry from the eastern magma-poor corridor of the 62°-65°E segment of the SWIR.
Including these possible alternative interpretations provides a much broader fit of our simulations to observations than constraining the results to perfectly uniform flip-flop detachment faulting.We thus hypothesize that while pure flip-flop mode may indeed represent a special case, a mix of the three modes discussed here provides a more robust explanation for magma-poor spreading in a very thick brittle lithosphere.This sheds new light on existing data and emphasizes the need for more densely sampled profile data, such as seafloor age measurements and hypocenter locations of microseismicity, also at other settings potentially hosting flip-flop detachment faulting.

Conclusions
Numerical experiments with an imposed temperature structure have demonstrated the key role of a temperature maximum in the central footwall of a detachment for focusing strain to trigger a new on-axis fault of opposite polarity.Using a dynamic thermo-mechanical model that accounts for the key mechanisms grain size evolution, serpentinization, enhanced hydrothermal cooling of fault zones, and influence of periodic sill intrusions at the base of the lithosphere, we have further explored the evolving faulting modes.Results of experiments in which we systematically varied magnitudes of magmatic input and hydrothermal cooling show that: • Fault zone cooling increases the temperature difference between the active fault and its surroundings.It thereby shifts the temperature maximum caused by the upwelling mantle material into the footwall.• Intrusion emplacement reinforces this temperature maximum and assists strain focusing through thermal and rheological weakening.• Stable flip-flop faulting that is controlled by the large-scale temperature field is limited to simulations with a fault-related long-term average heat flux about 20 times that of the Lost City Hydrothermal Field.Generally, flip-flop faulting is observed at reduced magnitudes of hydrothermal fault zone cooling and magmatic sill intrusions if both processes act together.• Two other faulting modes, antithetic accommodating footwall faults and criss-cross faults, offer an alternative interpretation of the current data base.
Our results indicate that the combination of weak magmatism at the base of the lithosphere and fault zone related hydrothermal activity can explain the thermal structure of the lithosphere and consequently the faulting mode of magma-poor mid-ocean ridges like the 62°E to 65°E segment of the SWIR, which might be more diverse than previously thought.

Figure 1 .
Figure 1.Initial temperature field, finite element mesh and mechanical boundary conditions.The geometry of the temperature field is controlled by the axial depth to the 750 and 1,260°C isotherms (brittle-ductile transition at 15 km and lithosphere-asthenosphere boundary 35 km, respectively), the off-axis depth of the 1,260°C isotherm of 60 km, and the width of the axial thermal structure of 55 km.
Figure2.We distinguish between the formation of individual faults (e.g., a detachment of opposite polarity or an accommodating footwall fault) and the fault sequence (e.g., successive parallel faults or flip-flop).To facilitate a comparison between our model calculations and the ones by Bickert et al. (2020)-despite using different numerical techniques and parametrizations of key processes-we first investigate a setup that is very similar to theirs.To do so we analyze the faulting dynamics of a model with imposed temperature structure, where the initial temperature field is modified only by the vertical movement of the domain top representing seafloor relief.We analyze the faulting dynamics of this model in Section 3.1.

Figure 2 .
Figure 2. (a) Sketches illustrating the different faulting modes observed in our simulations.The numbering indicates the order of the sequences.Colored arrows on top indicate direction of increasing age of the respective seafloor section.Abnormal age gradients (increasing age toward ridge axis) are indicated by an asterisk.In the chaotic mode, seafloor structures are frequently dismantled by the underlying faults, making a classification practically impossible.(b) Interpretation of fault zones and seafloor relief of the eastern amagmatic corridor of the SWIR by Cannat, Sauter, et al. (2019) (D3' from Corbalán et al. (2021)); vertical exaggeration 2:1; B-Breakaway, E-Emergence, #1 is the currently active detachment, C5-magnetic anomaly #5.Modified afterCannat, Sauter, et al. (2019).

Figure 3 .
Figure 3. Central model region displaying flip-flop dynamics in the model with imposed temperature structure.Left: strain rate and white isotherms, vector colors show advection velocity.Fault numbering as in Figure 4. Right: grain size in μm (blue to white) overlain by serpentinization degree (green to yellow) and black isotherms.The same graphical framework is used for all related figures.The full sequence can be seen in Movie S2.

Figure 4 .
Figure 4. Fault zone analysis of the model with imposed temperature structure showing the evolution of strain rate, fault width and grain size.Strain rate and fault width are evaluated in the brittle lithosphere, grain size at the detachment root around the brittle-ductile transition between 500 and 1000°C.Left panel: Time series of subsequent faults.Dashed vertical lines indicate the activation of a new detachment.Right panel: Same data for all detachments in the simulation but as a function of their life time, that is, shifted so that their point of activation (vertical lines in the left panel) match.The second vertical line indicates the onset of the following detachment, giving a mean fault period duration of 1.1 ± 0.1 Myr (SWIR: 1.1 ± 0.3 Myr; Cannat, Sauter, et al., 2019).

Figure 5 .
Figure 5. Sequence showing a model without any hydrothermal cooling (Nu = 1) and no magmatic input.The full sequence can be seen in Movie S3.

Figure 6 .
Figure 6.Sequences illustrating the different types of long-lived detachments observed in the simulations.Labels: B-breakaway of the long-lived detachment; ldlong-lived detachment; sa-synthetic accommodating fault; aa-antithetic accommodating fault.The full sequences can be seen in Movies S4 (a,b,e,f), and S5 (c,d).

Figure 7 .
Figure 7. Flip-flop sequence from simulation 28 including hydrothermal fault zone cooling and magmatic intrusions (M = 0.095, F max = 1.5).First and second column: see Figure 3. Intrusions are marked in red in the central panels.Third column: Nusselt number on a logarithmic scale, maximum value Nu max ≈ 200.Intrusions are marked in red in the central panels and by temperature in the right panels.The full sequence can be seen in Movie S6.

Figure 8 .
Figure 8. Fault zone analysis of a flip-flop sequence from simulation 28 including hydrothermal fault zone cooling and magmatic intrusions (M = 0.095, F max = 1.5).Analysis strategy as described for Figure 4.

Figure 9 .
Figure 9. Exemplary sequences from simulation 28 showing two different modes of deep-cutting faults.(a-c) spider mode: both deep-cutting faults remain symmetrically active.The "spider legs" are visible in (c) in the low strain rate regions on the side of the central block.They are even better visible in the grain size evolution in Movie S6. (d-f) Criss-cross mode: deep-cutting faults repeatedly cut the enclosed block and offset the surface relief.The full sequence can be seen in Movie S6.

Figure 11 .
Figure 11.Parameter space plot showing all 36 simulations of the parameter study as a function of magmatic input fraction (controlled by intrusion size and periodicity) and magnitude of fault zone cooling.Additionally, simulation 0 without any hydrothermal cooling and magmatic intrusions is shown.Pie diagrams show the different faulting modes identified in each simulation.Figure S2 in Supporting Information S1 breaks down the pie diagrams to the percentages of each mode during the 40 Myr model time.Since faults can form as one type and evolve to a different type, for example, synthetic accommodating faults to long-lived detachments, the total number of types per pie diagram generally adds up to more than the total number of faults per simulation.Simulation 32 is blurred because faulting takes place close to the domain boundary, that is, outside of the central well-resolved mesh region, and the simulation is therefore disregarded.Descriptions of the faulting modes are found in the text as is the conversion from F max to hydrothermal heat flux values at the fault zone.Literature values for the total heat output of Lost City (40 kW) and Longqi-1 (250 MW) vent fields are taken from Lowell (2017) and Tao et al. (2020), respectively, and are divided by their fault zone parallel extent.

Figure 10 .
Figure 10.Exemplary sequences from simulation 24 showing the chaotic faulting mode.The full sequence can be seen in Movie S7.

Figure 12 .
Figure 12. Analysis of the influence of different parameter combinations on the temperature field using the deformation from simulation 28.(a) Temperature difference between simulation 28 and its evolution when setting F max = 0 and M = 0 (same snapshot and fault numbers as in Figure 7b).White isotherms are from simulation 28, black ones are for F max = 0 and M = 0.The horizontal bars illustrate the meaning of the parameters plotted in panels (b) and (c).(b) Horizontal shift of the shallowest point of the 750°C-isotherm compared to the setup with F max = 0 and M = 0. Dashed vertical black lines indicate the formation of a new fault.Note that the temperature field is reset at 30.6 Myr when the new fault forms to facilitate this comparison.(c) Width of the 750°C-isotherm, calculated as the geometric mean between five slices at 1 and 5 km below the shallowest point of the isotherm.

Figure 13 .
Figure 13.Parameter space plot showing the analysis of faulting sequences.All three plots use the framework described for Figure 11.Triangular outlined region marks the parameter range producing stable true flip-flop faulting.Gray hatched region at F max = 2.0 is excluded from further analyses, see Section 4.3.2.(a) Proportion of faults within identified fault sequences that reproduce bathymetric and subsurface features without displaying classical flip-flop mechanics.(b) Proportion of faults within identified fault sequences that, additionally to the conditions for a, reproduce the flip-flop mechanism.(c) Mean fault duration of faults within identified flipflop sequences (apparent and true).The black hatched region gives mean general fault duration for simulations without SWIR/flip-flop resembling faulting sequences.

Table 1
Description and Values of Model Parameters and Variables Used for This Study