On the formation of oceanic detachment faults and their influence on intra-oceanic subduction initiation: 3D thermomechanical modeling

Abstract Extensional detachment faults, which have been widely documented in slow-spreading and ultraslow-spreading ridges on Earth, can effectively localize deformation due to their weakness. After the onset of oceanic closure, these weak oceanic detachments may directly control the nucleation of a subduction zone parallel to the former mid-ocean ridge, as is suggested for the Neotethys in Middle Jurassic times. So far, this hypothesis has only been tested by 2D numerical models, whereas the geometry of detachment faults is intrinsically three-dimensional. Here, we conducted a series of 3D numerical thermomechanical experiments in order to investigate the formation of detachment faults in slow oceanic spreading systems and their subsequent response upon inversion from oceanic spreading to convergence. Numerical results show that during the oceanic spreading stage, the formation of detachment faults strongly depends on the magnitude of the healing rate of faulted rocks in the oceanic lithosphere, that reflects the stability of hydrated minerals along fractured rocks. The detachment faults formed in our 3D numerical models deviate from the “rolling hinge model” of oceanic detachment faulting where fault footwalls are rotated and oceanic core complexes are thereby formed. Our results accentuate that the controlling physical parameters for the development of oceanic core complexes and detachment faults can differ, and that their coupled development in nature remains a key target for future research. Upon modeled transition to compression, previously formed asymmetric spreading patterns are prone to asymmetric inversion, where one oceanic plate thrusts under the other. Our results suggest that detachment faults accommodate significant amounts of shortening during the initiation of oceanic closure, but, in contrast to the previously proposed simple conceptual model, no direct inversion of a single detachment fault into an incipient subduction zone is observed. Instead, a widespread interaction of multiple detachment faults occurs after the onset of convergence. Ultimately, the nascent subduction zone cuts through the base of several pre-existing detachment faults, thereby forming an initial accretionary wedge in the incipient fore-arc.


Introduction
Subduction occurring at multiple locations on Earth is considered to be a key driving process of the modern plate tectonic regime (e.g., Gurnis et al., 2004;Stern and Gerya, 2017). Despite its importance, many aspects of subduction remain poorly understood. One of the most unresolved and controversial challenges in current geodynamic research involves the initiation of subduction zones (e.g., Gurnis et al., 2004;Stern, 2004;Stern and Gerya, 2017 and references therein). Based on theoretical considerations and natural data, Stern (2004) proposed two major modes of subduction initiation: spontaneous (i.e., caused by forces originating at the subduction initiation site and not elsewhere) and induced (i.e., caused by ongoing plate motion, or changes in plate motion caused by changes in force balance away from the subduction initiation site). In the first mode, cold, old lithospheric plates sink spontaneously in the mantle under their own weight due to their inherent gravitational instability (Vlaar and Wortel, 1976;Stern, 2004). The dynamics of the latter subduction initiation mode are governed by the interaction of forces that drive and resist subduction. In order to initiate an induced subduction zone, the externally applied compressional forces must overcome the resistive forces, i.e. the plates' elastic resistance to bending and the frictional shear resistance at pre-existing fault zones (McKenzie, 1977;Matsumoto and Tomoda, 1983;Gurnis et al., 2004). This constraint has been resolved in some cases by former numerical studies, such as the induced initiation of intra-oceanic subduction zones along transform faults and fracture zones in rift systems (e.g. Matsumoto and Tomoda, 1983;Gurnis et al., 2004;Stern and Gerya, 2017 and references therein).
The majority of recently initiated subduction zones are intra-oceanic, stressing the importance of a proper understanding of the various processes leading to modern intra-oceanic subduction initiation (Gurnis et al., 2004;Stern and Gerya, 2017). One particular intra-oceanic subduction initiation mechanism that has received attention recently concerns the inversion of spreading ridges to trenches (e.g., Agard et Accepted manuscript, Gülcher et al. (2019), EPSL, 506:195-208, https://doi.org/10.1016EPSL, 506:195-208, https://doi.org/10. /j.epsl.2018 a,l,. 2015; Maffione et al., 2015;van Hinsbergen et al., 2015;Duretz et al., 2016). Ophiolites overlying metamorphic soles are commonly thought to be generated after subduction has started near a ridge (Wakabayashi and Dilek, 2003), yet intra-oceanic subduction initiation parallel and close to oceanic spreading ridges remains controversial. Although young and thin lithosphere adjacent to ridges could potentially be easily deformed, a young lithosphere might not be negatively buoyant enough to favor subduction (Cloos, 1993) and a ridge push force will oppose compression. It is therefore not clear how easy it is to initiate subduction in this geodynamic setting (e.g. Stern and Gerya., 2017 and references therein).
Based on geochemical, geological and paleomagnetic data, the Mirdita ophiolite (Albania) is interpreted as a supra-subduction zone ophiolite that formed in the Middle Jurassic in the western Neotethys parallel to and near a spreading center in a forearc during subduction infancy (Stern and Bloomer, 1992;Maffione et al., 2015). Together with geological and paleomagnetic observations, the geochemical signatures of the Mirdita ophiolite led to the interpretation that the subduction zone initiated along an oceanic detachment fault parallel to the spreading center (Bortolotti et al., 2013;van Hinsbergen et al., 2015). These observations were supported by 2D numerical experiments undertaken by Maffione et al. (2015), suggesting that subduction can initiate along oceanic detachment faults formed during seafloor spreading.
Oceanic detachment faults are long-lived normal faults exposed at the seafloor and have been documented in many slow and ultraslow-spreading systems (e.g. Dick et al., 2003;MacLeod et al., 2002MacLeod et al., , 2009Escartin et al., 2008). They are often associated with compositionally and tectonically complex oceanic lithosphere. One of the most striking features are massifs in which lower crustal and upper mantle rocks are exhumed at the seafloor, better referred to as oceanic core complexes (e.g. Cann et al., 1997;Escartin et al., 2003;Smith et al., 2008). Detachment faults can cut the lithosphere up to a depth of 7-8 km , which allows seawater to infiltrate to greater depths. Mantle lithosphere can rise up to shallower levels, resulting in greater hydrothermal alteration of the surrounding lower crust and mantle rock, modifying its texture and forming secondary hydrous minerals -such as serpentinite and talc (Boschi et al., 2006). As a result of the weakening effects of these hydrous minerals, detachment faults are associated with a substantially decrease of the overall strength of the oceanic lithosphere (McLeod et al., 2002;Boschi et al., 2006), which may encourage the localization of plastic strain upon compression.
Detachment faults are therefore a likely candidate for the initiation of subduction when compressional far-field forces are applied .
The hypothesis of intra-oceanic subduction initiation controlled by detachment faults requires further testing. So far, this subduction initiation mechanism has only been investigated with 2D numerical models with a prescribed, rather than spontaneously formed, detachment fault .
Furthermore, the geometry of detachment faults in oceanic spreading systems is intrinsically threedimensional, demanding a 3D modeling approach (e.g., Tucholke and Lin 1994;Escartin et al., 2003Escartin et al., , 2008, whereas most previous explorations of faulting in extensional systems were done in 2D (e.g. Lavier et al., 1999;Beaumont, 2002, 2003;Buck and Poliakov, 1998;Buck et al., 2005). Detachments represent relatively small structures in oceanic basins and it remains unclear how discrete detachments can link upon compression to form a continuous, wide subduction zone . Thus, the formation and geometry of detachment faults in a three-dimensional slow-spreading oceanic basin as well as how detachments can initiate a self-sustaining subduction zone parallel to the ridge axis, are left unresolved.
In this paper, we conducted a series of 3D numerical thermomechanical experiments in order to investigate the formation of detachment faults in slow oceanic spreading systems and their subsequent response upon inversion from oceanic spreading to convergence. The research is divided into two main parts: (i) Improve our understanding of the formation of detachment faults in (ultra)slow oceanic spreading systems in three-dimensions. What physical parameters control their 3D geometry and spatial distribution?
(ii) Investigate the response of a slow-spreading rift system with detachments upon compression.
Under which conditions will subduction be initiated and how is the incipient subduction suture affected by the detachments?

Numerical approach
Numerical experiments were performed using the 3D thermomechanical I3ELVIS code (Gerya and Yuen, 2007;Gerya, 2013), based on a combination of a Finite-Difference method applied on a fully Material properties are advected on the moving Lagrangian markers.

Model set-up
The initial model set-up is shown in Figure 1  boundary condition is approximated on top of the crust by defining a 5 km sticky water layer (Schmeling et al., 2008). A high thermal conductivity is used for the upper sticky water layers to ensure an efficient heat transfer at the upper surface of the plate. Spreading, and subsequent convergence, is assumed to be caused by far-field forces and apply on all rock types. Boundary conditions are constant symmetrical spreading rates in x-direction where "#$%&'()* = ,-%./ + ,$(*1/ and ,-%./ = ,$(*1/ (see Fig. 2).
Compensating vertical influx velocities which are constant along the whole upper and lower boundaries are derived from , ( ) to ensure conservation of model volume as well as constant average thickness of the sticky water layer. Because detachment faults form in ultraslow to slow-spreading oceanic ridges (e.g. Dick et al., 2003;Smith et al., 2006), full spreading rates were taken in a limited range of 1.5 to 2.5 cm/yr, characteristic of ultraslow to slow-spreading systems (Dick et al., 2003). Inversion from divergence to convergence is implemented by a linear interpolation of the prescribed velocity boundary condition from Inversion is set to happen after 5 million years of oceanic spreading, such that a mature spreading pattern can first be formed at the ridge. The transition time from extension to compression was set to 3-5 Myr to allow for gradual changes in deformation.

Numerical implementations of oceanic spreading processes
In order to accurately model plate break-up and oceanic spreading, an implementation of the successive melting/accretion-related processes is required. The following key processes were suggested to be critical for this purpose: (i) thermal accretion of the oceanic mantle lithosphere resulting in plate thickening; (ii) partial melting of the asthenospheric mantle, melt extraction and percolation towards the ridge resulting in crustal growth; (iii) crystallization of the new oceanic crust beneath the ridge and (iv) hydrothermal circulation at the axis of the ridge resulting in excess cooling of the crust (e.g. Buck et al., 2005;Gregg et al., 2009;Katz, 2010). These processes are included in our numerical model in a simplified manner as described by Gerya (2013). It is important to note that the adopted magmatic model is designed to reproduce a gross-scale distribution of melt production processes below a ridge but cannot replicate the exact volcanic/magmatic activity, which is commonly observed around slow-spreading ridges in nature (Gerya, 2013).

Rheology
A constant, low viscosity of 10 18 Pa s is implied for the partially molten crustal and mantle rocks, whereas a non-Newtonian visco-plastic rheology is implied for the lithospheric plates. Upper and lower cut-off viscosity limits for the lithosphere were set to 10 24 and 10 18 Pa s. The viscous rheology contains both stress-independent diffusion creep and stress/strain-rate dependent dislocation creep (Ranalli, 1995). Brittle-plastic rheology of the plates assumes fracture-related strain weakening (e.g. Lavier et al., 2000;Huismans and Beaumont, 2002;Hieronymus, 2004;Gerya, 2013). This is implemented by using a Drucker-Prager yield criterion defining the visco-plastic transition: where AA is the second stress invariant (Pa), is the dynamic pressure on solids (Pa), . is the hydrostatic fluid pressure ( . = . with . = 1000 kg/m 3 , g=9.81 m/s 2 and =vertical coordinate), is the rock strength at − . = 0 and is the internal friction angle coefficient. Strain weakening is obtained by strain-dependent linear interpolation between parameterized minimum and maximum values of and where Q and T are the initial and final strength values for the fracture-related weakening, respectively.
≥ 0 is the integrated plastic strain (with Q = 1 as the upper strain limit for fracture related weakening) and is calculated via: here, ( Z is the plastic strain rate tensor and 1 %&-()* is the fracture healing rate. Fracture-related weakening is therefore implemented as a linear interpolation from Q to T and from φ=0.6 to φ=0 for ≤ Q . This assumed strain weakening model is adopted from previous studies of mid-ocean ridges and can correspond to multiple mechanism of weakening such as formation of weak hydrated phases by fluid percolation (Escartin et al., 2003) and/or fault gouge formation (Huismans and Beaumont, 2002). The amount of strain weakening is mainly investigated in this study through the ratio of Q to T , hereafter termed as the cohesion ratio. An increase of cohesion ratio may represent a decrease in coherency of brittle material due to initial pressure-independent fracturing, and vice versa (Eq. (5); Allken et al., 2012;), EPSL, 506:195-208, https://doi.org/10.1016/j.epsl.2018.042 Gerya, 2013). Partial strain healing of deactivated shear zones is applied by imposing a constant plastic healing rate 1 %&-()* that reduces the accumulated plastic strain with time. When the second strain rate invariant T W (( Z ) W is greater than the plastic healing rate for 1 %&-()* , strain weakening occurs, whereas strain healing of the faults occurs when the second invariant is smaller than 1 %&-()* . Following former numerical investigations of spreading systems in 3D, magnitudes of the healing rate were taken in the range of 0 -10 -13 s -1 in the numerical experiments (from Gerya, 2013;Puthe and Gerya, 2014). In cases where 1 %&-()* is set to zero, maximum weakening of fractured rocks occurs, and these models thereby represent spreading systems in which hydrated minerals are able to permanently weaken fractured rocks.
The upper limit of 1 %&-()* = 10 -13 s -1 reflects a fast healing of deactivated fractures, on the scale of geodynamic processes, which undergo no chemical change to a weaker phase (Lyakhovsky and Ben-Zion, 2008).

Numerical results
The conditions and results of the performed numerical experiments investigating spreading patterns and inversion are shown in table 2. In these experiments, different initial and boundary conditions were adopted: (i) spreading/inversion rates; (ii) inversion timing (iii) healing rate 1 %&-()* and (iv) the cohesion ratio C0/C1 of solid rocks. The results of modeling of the 3D (ultra)slow oceanic spreading patterns and the responses of these spreading systems upon inversion will be discussed in the sections below.
3.1 Development of 3D oceanic spreading patterns

Reference model development
The evolution of the reference model (diaa) during 5 million years of oceanic spreading is presented in Figure 2. In this reference model, the healing rate was set to zero and the C0/C1 cohesion ratio to 10/3.
The initial straight ridge produces two conjugate normal faults bounding half-grabens in the brittle lithosphere along which deformation localizes. Asymmetric accretion occurs as more motion is accommodated by one of the two faults, which cuts deeper into the lithosphere and is more pronounced than the other (Fig. 2a). This detachment fault localizes a significant amount of extensional slip in its initial stage. Another inherent feature is that asymmetric accretion spontaneously varies along the ridge axis.
The dip direction of the dominant detachment fault is abruptly reversed at z=75 km (Fig. 2a). This 50° angle to the ridge, which is visible in the bathymetric profiles (Fig. 2d). We refer to this bathymetric pattern as a 'Christmas tree pattern'.
Even though many detachment faults are formed, mantle rock is never exhumed up to the surface, and thus no oceanic core complexes are formed throughout the evolution of the model. This is because new detachment faults are formed closer to the ridge before the faults accumulate sufficient offset to exhume mantle rocks.  Increasing the C0/C1 cohesion ratio produced models with fewer detachment faults, which have a larger variability in shape and length and are separated by a greater distance of 10 to 15 km (Figs. 3a-c).

Influence of model parameters
Irregularity of the ridge is also promoted by larger cohesion ratios, and segmentation of the ridge by ridgeperpendicular transform offsets occurs in the models combining large cohesion ratios with low healing rates (Figs. 3c,f).

Reference model development
The evolution of the reference model (diaa) during inversion from spreading to convergence is presented in Figure 4. Inversion took place from 5-8 Myr. The mature spreading pattern after 5.0 million years of spreading is shown in Figure 4a. Many detachment faults up to 80 km long run parallel to the spreading ridge and localize extensional deformation. As the spreading rate decreases during the initial stage of inversion, the temperature decreases near the ridge, leading to lithospheric thickening (Fig. 4b).
Just after plate motion is inverted to convergence (Fig. 4c) faults, which interact with one another within a nascent accretionary wedge. The relatively narrow highresolution models cannot be used to further investigate the incipient subduction as the retreating trench starts to interact with one of the lateral boundaries (Fig. 4e).

Influence of model parameters
The inversion results at time t=9.0 Myr for models with different 1 %&-()* values are shown in Figure   5. Symmetrical responses upon compression are encouraged by a large healing rate of 10 -13 s -1 : upon inversion, the model diac with an initially symmetrical, smooth spreading pattern shows two active conjugate thrusts formed on both sides of the spreading ridge. The ridge swell forms atop the original spreading ridge and evolves in a symmetrical pattern. In the three other models with lower fracture ), EPSL, 506:195-208, https://doi.org/10.1016/j.epsl.2018 healing rates, an asymmetrical pattern with underthrusting forms. As follows from the comparison of these three models, an increase in fracture healing rate does not change the mechanism of formation but influences the geometry of the accretionary wedge-like structure and the morphology of the thrust. The pre-existence of regularly spaced deep detachment faults (models diaa and diad, Figs. 5a,b) encourages a stable horizontal underthrust associated with an accretionary wedge in which the detachments faults interact with one another to form back-thrusts. Without the presence of such detachment faults (model diab, Fig. 5c), the downgoing plate plunges directly into the mantle at a steeper angle and no evident accretionary wedge forms.

Detachment faults in slow-spreading systems
Detachment faults and off-axis strain localization have been documented in many of the conducted experiments. In our 3D models, oceanic spreading geometries have shown to have a significant dependence on the value of the healing rate 1 %&-()* of fractured rocks and cohesion ratio Q / T , as illustrated in Figure 3. The numerical results indicate that the formation and preservation of detachment faults in oceanic slow-spreading systems are mainly promoted by low to zero healing rates (1 %&-()* ≤ 10 -15 s -1 ). The low healing rate would correspond to hydration of the fault zone through fluid infiltration, which is unlikely to heal (dehydrate) during seafloor spreading and subduction initiation (Escartin et al., 2003). Detachment faults formed in our models cut the lithosphere to depths of 7-8, which would allow for hydration along the fault zone  and thereby the results are consistent with this suggested relationship between low healing rates and hydrothermal alteration of fractured rocks. The pressure-independent fracture weakening (i.e. cohesion ratio), mainly affects the geometry and spacing of the detachment faults but not their presence or absence in spreading patterns (Fig. 3).
Oceanic spreading models with ample healing of fractured rocks are inherently symmetric  whereas increased strain weakening breaks this symmetry and promotes asymmetric accretion. (Figs. 3a-i). The positive feedback between strength reduction and increased strain, leading to further weakening of fractured rocks, and its direct correlation with the development of asymmetric extensional systems was acknowledged in former 2D numerical studies (e.g. Beaumont, 2002, 2003) and is endorsed by our 3D results. The lifetime of faults is also greatly affected by weakening of fractured rocks: increasing the healing of fractured rocks decreases the average lifetime of faults and vice versa (Puthe and Gerya, 2014), which is encouraged by our models in which the lifetime of detachment faults decreases from 0.6 Myr in the reference model (diaa,Figs. 2;3a) to zero when the healing rate is too large to maintained long-lived weakening of the lithosphere (models diad, diap, diat; Figs. 4j-l). The influence of the coherency of the rock, i.e. the cohesion ratio, in our 3D numerical models are consistent to former 2D numerical studies of extensional faulting in two ways: (i) the correlation of a larger cohesion drop, i.e. cohesion ratio, and greater normal fault offsets (Buck et al., 1998;Lavier et al., 1999;Huismans and Beaumont, 2003), and (ii) larger spacing between normal faults for models with greater cohesion ratios (Huismans and Beaumont, 2003). The third dimension of our models reveals that numerical experiments with an increased cohesion ratio are more prone to transformal offsets, in addition to wider spacing of detachment faults. The evolutionary sequence of detachment fault formation and migration away from the ridge also agrees with previous 3D numerical studies of (ultra)slow oceanic spreading systems (Puthe and Gerya, 2014). Yet, the displacement accommodated along the detachment faults is typically moderate, and these structures do not exhume mantle rock to the surface. This faulting geometry deviates from the "rolling hinge" model of detachment faulting (Lavier et al., 1999), in which detachment faults initiate as steeply dipping normal faults but then rotate to shallower dip as a result of flexure in the footwall. It is suggested that the development of long-lived, low-angle detachment faults and accompanying oceanic core complexes require that a small part (~30-50%) of the total extension is accommodated by magmatic dyke accretion while the rest is accommodated by lithosphere thinning (e.g. Buck et al., 2005;Behn and Ito, 2008;Tucholke et al., 2008). As our model does not account for any dyke- ), EPSL, 506:195-208, https://doi.org/10.1016/j.epsl.2018 related processes, the entire oceanic crust is accreted by the cooling of the roof and the walls of a large shallow magma reservoir that forms spontaneously under the ridge as the result of melt extraction from the mantle (Gerya, 2013). Our numerical set-up thereby deviates notably from the previous 2D models, which could explain the observed difference in detachment fault geometries in numerical results.
The three-dimensionality of our models does provide a fascinating new aspect to slow-spreading systems. The "Christmas tree" bathymetric pattern, generated consistently in many of the numerical experiments (Fig. 3), and represent variation along the ridge in asymmetric accretion, expressed as a reversal in the dip direction of active detachment faults. The patterns require asymmetrical accretion along a bended ridge, and are thereby created in models with low healing rate values (1 %&-()* ≤ 10 -15 s -1 ) and relative low cohesion ratios (C0/C1 = 10/3).

Comparison of oceanic spreading models with natural data
Several common features of natural oceanic spreading systems have been reproduced numerically in this study. Orthogonal ridge-transform patterns are common in slow oceanic spreading ridges on Earth (e.g. Gerya, 2013 and references therein) and are formed spontaneously in several of our experiments (models diam, diaq, dian, diar; Figs. 3b-f). The "Christmas tree" bathymetric pattern (oblique patterns ~50° to the spreading ridge offsetting detachment faults) is systematically generated in the remaining part of the numerical experiments that accrete in an asymmetrical way (Fig. 3). Although the exact physics behind its development remains to be explored, the characteristic pattern also has analogous observations on spreading systems on Earth. Similarities between natural observations and several numerical models are shown in Figure 7. Arabian Sea (Kattoju et al., 2015).
Most of the experiments show rough topographies with highs and lows running parallel to the ridgeaxis. These structures are very similar to abyssal hills, one of the most prominent features of slowspreading mid-ocean ridge flanks (Buck et al., 2005). The abyssal hill topographic pattern is interpreted as a product of extensional faulting at slow mid-ocean ridges and has been mapped in many places on Earth.
Ridge and the Southeast Indian Ridge. From these bathymetric profiles, an average separation between abyssal hills of 5-7 km was measured. This wavelength agrees with the average distance between the detachment faults generated in the numerical models with abundant weakening of rocks at low cohesion ratios (Figs. 3a and d).

Ridge inversion
Numerical results shown in Figure 5 indicate that (ultra)slow-spreading ridges with asymmetrical faulting patterns manifested by the presence of deep detachment faults are prone to develop subductionlike underthrusting when ridge-perpendicular compressional forces are applied. Inversion of the spreading models show that underthrusting never localizes at the spreading ridge, but instead at one of the pre-existing detachment faults, and to that extend endorse the findings by numerical explorations by Maffione et al. (2015). Yet, none of our 3D experiments showed the direct inversion of a single detachment fault into an incipient subduction zone (as proposed by Maffione et al., 2015). The suboptimal orientation of the reactivated detachment fault (i.e. steep dip angle and limited lateral extent along the ridge) prevents direct plunging of the down-going plate along this single detachment fault. Instead, the thrust plane choses a lower angle geometry and thereby cuts through the base of several former detachment faults which interact with one another to form an accretionary wedge-like structure (Fig. 4e).
Also, in the absence of exhumed mantle rock that forms an oceanic core complex in the footwall, the gravity anomaly of the underthrusted lithosphere is insufficient to induce a bending of the oceanic lithosphere. This explains why the down-going lithosphere is underplated below overriding plate and stagnates around 20 km depth.
Although the numerical models do not allow for the development of a full-fledge subduction zone, they can be partially correlated with the geological record of ophiolites in Albania and Greece (Bortolotti et al., 2013). For simplicity we will only discuss correlations with the Mirdita ophiolite as the mechanism of nearridge subduction initiation through detachment fault inversion was based on this ophiolite only (Maffione et al., 2015 and references therein). A major argument in favor of the validity of the models in this study is that they reproduce the high peak pressure-temperature conditions recorded in the metamorphic sole at the time of subduction initiation (624°C to 769°C with pressure <0.7 GPa; Gaggero, 2009) as well as the inverted metamorphic gradient (assuming that the base of the hanging wall in the model corresponds to the metamorphic sole; Fig. 4). In many of the numerical experiments, the downgoing plate is very fragmented during oceanic closure, thereby forming an accretionary wedge atop the main thrust mainly composed of oceanic lithospheric fragments (Figs. 5a-b). This fragmented prism cannot be observed in the Mirdita ophiolite or in the Melange below, which is dominated by sedimentary fragments from the Adria continental margin (Bortolotti et al., 2013), but it is likely that these fragments have disappeared through the process of subduction erosion by the time of obduction and that only a few remnants are preserved in the mélange below the ophiolite. Lastly, there are significant differences in the way the detachment faults are reactivated in the numerical models and the geological record of the Mirdita ophiolite. In our results, it is the detachment fault closest to the ridge that is reactivated during compression (Fig. 4). The Mirdita ophiolite contains a fossil detachment fault that is interpreted as closer to the ridge but not reactivated during subduction initiation (Maffione et al., 2013). This discrepancy might be because, unlike in nature (Maffione et al., 2013), the detachment faults in our model are short-lived and do not form large-slip detachments that expose oceanic core complexes (e.g. Lavier et al., 1999).
Large-slip detachment faults might be weaker (e.g. more strain weakening) and reach deeper in the lithosphere and therefore will be more likely to reactivate in inversion (Buiter et al., 2009) than short and steep detachment faults closer to the ridge. More studies would be necessary to properly elude this question. In particular, taking into account many similarities between naturally observed and numerically reproduced oceanic spreading patterns, our simulated ridge inversion patterns with an accretionary prism composed of multiple seafloor fragments should be of rather broad significance and call for further numerical studies and systematic comparison with other known ophiolite complexes worldwide (e.g., Vaughan and Scarrow, 2003).

Limitations and outlook
It is important to note that our numerical experiments were performed on a relatively small spatial scale with a high resolution and a simple geometry, aiming to form an initial basis for large-scale formation and inversion of complex spreading systems in 3D. Due to the model dimensions, we are limited in our ability to correlate our relatively small-scale and short-term ridge inversion models to geological observation in ophiolites. Also, the amount of induced plate convergence (<50 km) in the numerical models does not span the entire period of subduction initiation (e.g. Stern and Gerya, 2017 and references therein). Therefore, further high-resolution 3D numerical experiments investigating larger segments of an oceanic spreading system and the underlying mantle are still needed to develop a better understanding of detachment fault-controlled subduction initiation.
Our results emphasize the influence of the healing and weakening mechanisms occurring in faulted rocks on oceanic spreading patterns. It must be noted that the employed simplified strain weakening/healing model does not directly reflect the physics involved in the various weakening/healing processes, such as fluid percolation, hydration reactions and grain damage or growth. As fluid-rock interactions are greatly associated with detachment tectonics (e.g. Escartin et al., 2008), further theoretical development and numerical implementation of a more sophisticated and realistic lithospheric weakening/healing model might bring a better understanding of the relation between detachment faults, oceanic core complexes and subduction initiation.

Conclusions
We performed high-resolution three-dimensional numerical modeling of (ultra)slow oceanic spreading systems with detachment faults and their response upon subsequent ridge-perpendicular convergence.
Our numerical models suggest that asymmetric spreading and detachment fault development is mainly controlled by the low fracture healing rates of brittle/plastic lithospheric rocks. This low fracture healing rate reflects ample strain weakening along faults, possibly driven by hydrothermal alteration of the fractured rocks. The modeled spreading patterns share key features with former 2D and 3D numerical studies as well as natural systems, such as asymmetrical rift development, ridge-parallel abyssal hills and transformal offsets. Numerical results deviate from the simple conceptual "rolling hinge" model of oceanic detachment faulting, where a long-lived main detachment fault is intrinsically associated with the development of an oceanic core complex (Lavier et al., 1999). Instead, detachment faults are regularly spaced and have moderate offsets. The melt supply by dykes in the spreading systems has not been assessed in this study, and we propose that in the future studies the magmatic vs. tectonic extension ratio should be coupled with the fracture weakening/healing parameters, to improve our understanding of extensional faulting behavior in slow-spreading systems.
Upon ridge-perpendicular convergence, the pre-existing detachment faults have been shown to localize extensive deformation. Our results indicate that asymmetric spreading patterns with multiple detachment faults are prone to asymmetric inversion, where underthrusting of one oceanic plate under the other occurs. Instead of a single detachment faults directly inverting into an incipient subduction zone, multiple ridge-parallel detachment faults interact with one another upon compression. The main thrust plane cuts through the base of several detachment faults, thereby forming an incipient accretionary wedge in the incipient fore-arc.
The formation of a continuous, self-sustaining subduction zone parallel to a pre-existing mid-ocean ridge remains an outstanding large-scale three-dimensional problem and future numerical modeling research is needed to account for larger spatial and temporal scales. Key targets for future threedimensional modeling is (i) understanding the formation of and coupling between detachment faults and oceanic core complexes in slow-spreading systems, possibly by assessing both the melt supply as the healing/weakening parameters; and (ii) the effect of detachment faults on the initiation of a selfsustaining, laterally continuous subduction zone.