CERTIFER Jean--Louis Boulanger
email: jean--louis.boulanger@certifer.eu

The new CENELEC EN 50128 and the used of formal method

Keywords: Embedded, CENELEC 50128:2011, Formal Method, Model checking, Proof, Safety, Software, SSIL, Tool, Qualification, V&V

The standard CENELEC 50128 [CEN 01, 11] identifies a complete process for software development of railway application. The new version 2011 introduced many new needs and a complete new structure. In this paper we present the new version of the CENELEC EN 50128 and describes how we can instantiate it. This new version introduce some new activities such the tools qualification, the software deployment, etc. and develop some activities such the data preparation and the software maintenance.

In railway, we used from many year some formal methods for the specification, the conception and the verification by proof or model checking. This paper present the new CENELEC 50128 ([CEN 11]) and describe how the formal method can be used and what is impact on the recommended activities.

CENELEC 50128

Introduction

The CENELEC EN 50128 [CEN 01a] standard is specifically dedicated to the development aspects of software for the rail sector. Regarding software, the SSIL (Software SIL) makes it possible to define different levels of criticality: from 0 (no danger) to 4 (critical). The CENELEC EN 50128 [CEN 01,CEN 11] standard specifies the procedures and the technical prescriptions applicable for the development of programmable electronic systems used in applications for rail command and protection. The CENELEC EN 50128 standard is thus not normally applicable to all software applications of the rail sector.

Figure 1 introduced the structure of the new standard CENELEC 50128:2011. This standard is based on a new notion the software assurance. The software assurance is composed with the quality management, the verification and the validation, the assessment and the tools qualification. For the clause 7 that described the development of generic software, is possible to use the formal method at different level: specification, architecture and design but also for verification. The clause 8 concern the application data and introduce a process for manage the data preparation process.

Application

The specification of a software application is thus at the very least, a set of requirements. A first analysis of the specifications provided by the client must make it possible to identify functional needs. The difficulty resides then in the definition of the concept of a requirement. There are several studies that attempt to identify what a requirement is and how to account for it. [HUL 05] presents one of the most complete syntheses. In parallel with the identification of the functional needs, it is possible to begin analyses linked to dependability, the objective of which is to define the nonfunctional requirements: safety requirements but also availability, reliability, or other requirements. It is, however, necessary to know a bit more about the needs of the software application. Within the application, it is thus necessary to introduce:

-interfaces with the environment;

-the notion of state in establishing a partition between safe functioning, decline, and dangerous states;

-the notion of correct behavior and of dangerous behavior;

-the notion of requirement linked to safety; this type of requirement must allow for characterization of dangerous behavior;

-the integrity level of the software (written SSIL/DAL).

The environment of the software application is composed of interfaces with the hardware resources (memory, specific address, input/output, watchdog, etc.), with other software applications (base software, related application, etc.) and/or with the operating system.

Figure 2 reveals a software application environment which is composed of three entries (Ei), two exits (Sj) and three interfaces (Ik) with the hardware resources (example access to a specific memory address).

Standards recommend that the specification of a software application be composed from a textual description of the need (the requirements) and all the notations necessary for facilitating understanding of the need. So, for the EN 50128 [CEN 01, CEN 11] standard, there is table A.2.

Classically, designers of a software application go directly from the textual requirements to the code without having been able to verify coherence of the requirements and without always mastering the unity of requirements. -a formalization of requirements phase (see figure X.10). This formalization phase can rely on structured methods, a model or formal methods (controller, Petri net, Grafcet, B-method, SCADE, etc.); -an architecture phase.

which can be translated into set-theoretic logic of the first order in ∀ t1,t2 ∈ T , hence Dt ∩Dt =φ , if t1 ≠ t2 , with Dt which is the domain of train speed i and [T] which is the whole set of trains. In the railway domain, structured method (SADT and SART for example), semi--formal method and formal method (B method, SCADE, etc.) are used to formalize the need (set of requirement). More generally, we used some model.

Model provide the capability :

-to manage the need; -to formalize the need and for verify the coherency and the completeness; -to help for tests case selection; -4 Formal method

Definition

Section B.30 of the CEI/IEC 61508 standard indicates that a formal method (HOL, CSP, LOTOS, CCS, linear temporal logic (LTL), VDM 1 [JON 90], and Z 2 [SPI 89]) enables an unambiguous and coherent description of a system at a stage of development (specification, architecture, and/or design) to be produced. The description takes a mathematical form and can be subjected to a mathematical analysis. This mathematical analysis may be performed automatically.

A formal method generally offers a notation, a technique for processing a description in this notation, and a verification process for controlling correction of the requirements. 1 In the IEC 61508 standard, the references to VDM include VDM++ [DUR 92], which is a real--time and object--oriented extension of VDM. To find out more, visit: www.vdmportal.org 2 In the IEC 61508 standard, the B--method [ABR 96] is seen as a method associated with Z.

…" Req_11:"…" Req_12:"…" Req_13:"…" …"

The CEI/IEC 61508 standard indicates that it is possible to make transformations right down to "a logic circuit design" 3 . Petri nets and state machines (mentioned in the outline of semi--formal methods) can be considered as formal methods, according to the degree of conformity to a rigorous mathematical basis of their uses.

Application

Two types of approaches were presented in different chapters:

• the first consists of starting from a specification to create a formal model (see Figure 5) and to carry out verifications on the model; • the second consists of carrying out formal analyses on a code carried out traditionally (in C, ADA, or C++ for example) starting from a specification (see Figure 6). …" Req_11:"…" Req_12:"…" Req_13:"…" …" int"xx;" main"()" {" …" }"

Analyse"

Impact

Since the first version of the standard CENELEC EN 50128 in 2001, formal methods were introduced and highly recommended. In the new version (2011) of the standard, the V--cycle is recommended and testing is the basic verification technique. But from specification to design it is possible to use formal methods.

The standard introduced the possibility to use the formal techniques (proof, model--checking, etc.) in place of tests for the verification (see the table A.5 of the CENELEC EN 50128).

For verification, you can choose a combination of techniques, the standard proposes some best choices (1 and 4, 3 and 4 and the last 4,6,7), and if you introduce a new combination you need to explain why it's a good choice and why this set of techniques have the same efficacy.

Each time, formal methods and/or formal techniques are used; you need to explain why the efficacy (for some error class) is similar.

In railway, formal methods are used more and more and at different levels: specification, architecture, design, in replacement of unit tests and of software/software integration test, in data preparation, etc. Other works, such as [MON 00, HAD 06], provide a panorama of the formal methods in a more scientific approach.

Return of experience

In the first book "Static Analysis of Software -the abstract interpretation" [BOU 12a], we present some formal tools based on the static analysis and abstract interpretation of code.

In the second book [BOU 12b], we presented some approach based on the modelization (SCADE, B Method, Mathlab/Simulink, and ControlBuild). In third book, several environments and formal methods were presented: Spark Ada, Matelo, AltaRica, Polyspace, Escher Verification Studio Perfect Developer, SCADE,B method.

Actually, in all railway projects we used some modelisations based on SCADE, CONTROL--BUILD, B--Method at different levels (system, software and complete software or for some specific function). For some project, the proof of properties is applied with a reduction of some tests activities.

One of the big used of proof in replacement of the test activities, is related to the data validation. The railway software manages many data call configuration data and the main efficient technics to validate these data are the proof ([BOU 07]). One of the topics discussed is related to the acceptance by authority (what is the evidence, how formalize some proofs, etc.) and the impact on the certification.

Conclusion and Future Work

Figure 1 :

 1 Figure 1 : CENELEC EN 50128.

Figure 2 :

 2 Figure 2 : Software application environment.

Figure 3 .

 3 Figure 3. From the requirements to the code

Figure 4 .

 4 Figure 4. Formalization of requirements

Figure 5 :

 5 Figure 5 : Requirement and model.

Figure 6 :

 6 Figure 6 : Formal analysis of a code.

  This paper presented the new standard CENELEC 50128:2011 and introduced some discussions on the impact of the used of formal method in development and in the verification of the generic software and on the data preparation. Some examples of real used are presented and explain what are the difficulties (what is a metric, what is the code coverage for graphic view, what is the set of document to produce and why, etc.).

It should be noted that initially SCADE was a development environment based on the language LUSTRE and that since Version 6, SCADE became a language of its own (the code generator for Version 6 works well inputting a SCADE


Blandine Gardonio
email: gardonio@geologie.ens.fr

Romain Jolivet
email: romain.jolivet@ens.fr

E Calais

H Leclère

M w 6.5 Botswana Earthquake: An Intraplate Event Triggered by Deep Fluids

Our InSAR-based estimate of the source parameters indicates a deep event with two probable families of models.

• We detect a precursory phase using template matching at teleseismic distances to the mainshock.

• A transient pulse of fluids is likely required to trigger this event, maybe a common stress perturbation at the origin of intraplate events

While earthquakes primarily occur along plate boundaries, where most of the tectonic strain is released, large events also strike stable continental interiors, although much less frequently [e.g. [START_REF] Johnston | The Seismicity of 'Stable Continental Interiors[END_REF][START_REF] Calais | A new paradigm for large earthquakes in stable continental plate interiors[END_REF]. The New Madrid region of the Central U.S., with four M>7 events between December 1811 and February 1812, is a type locale for large earthquakes in such settings [e.g. [START_REF] Nuttli | The Mississippi Valley earthquakes of 1811 and 1812: Intensities, ground motion and magnitudes[END_REF][START_REF] Johnston | Seismic moment assessment of earthquakes in stable continental regions-II. Historical seismicity[END_REF][START_REF] Hough | On the Modified Mercalli intensities and magnitudes of the 1811-1812 New Madrid earthquakes[END_REF], and examples can be found in all continents. South Africa was struck in 1969 by the M w 6.3 Ceres earthquake [START_REF] Kruger | The 29 September 1969, Ceres, South Africa, Earthquake: Full Waveform Moment Tensor Inversion for Point Source and Kinematic Source Parameters[END_REF], and in 1998, the Tennant Creek sequence of M s 6.3, 6.4, and 6.7 events shook central Australia [START_REF] Bowman | The 1988 Tennant Creek, northern territory, earthquakes: A synthesis[END_REF]. Careful studies of such rare events are important as the mechanism leading to stress release on faults in regions of very low tectonic deformation remains poorly understood, leading to large uncertainties in hazard assessment in populated continental interiors [e.g. [START_REF] Allman | Increasing Loss Potential in Earthquake Risk-A Reinsurance Perspective[END_REF][START_REF] Liu | Mid-continental earthquakes: Spatiotemporal occurrences, causes, and hazards[END_REF]. We focus on the largest stable continental interior earthquake since the 1998 Tennant Creek events, a M w 6.5 normal-faulting earthquake that ruptured a blind fault in Botswana on April 3rd, 2017, within the Kalahari cratonic domain, far from any identified active fault (Fig. 1).

The M w 6.5 Bostwana earthquake occurred close to the junction between the Archean Kaapval craton and the Paleoproterozoic Mahalapye granite of the Limpopo belt (Fig. 1).

It may have reactivated one of the high angle faults that mark the boundary between these old and seismically quiet geological units [START_REF] Kolawole | Aeromagnetic, gravity, and Differential Interferometric Synthetic Aperture Radar analyses reveal the causative fault of the 3 april[END_REF]. Seismic activity in its epicentral region (-25 • to -21 • south, 23 • to 27 • east) is more elevated than typical stable continental interiors, with 23 earthquakes of magnitude larger than 4 detected since January 1st, 2000. The mainshock was followed by about 30 aftershocks of magnitude greater than 3, as reported by the International Seismological Center (ISC). Aftershocks appear to first extend over a large region before focusing closer to the epicenter, though this may be biased by the location capability of small events at teleseismic distances. One of the objectives of this study is to expand this aftershock catalog and, more importantly, to search for potential foreshocks, using for the first time a template matching technique at teleseismic distances from the mainshock.

Although southern Africa experiences diffuse, low magnitude seismicity [e.g. [START_REF] Pastier | Is the Okavango Delta the terminus of the East African Rift System? Towards a new geodynamic model: Geodetic study and geophysical review[END_REF], no active structure has been identified in the epicentral region of the Botswana earthquake. The nearest active faults are 350 km to the north bounding the isolated Okavango tectonic basin, where a M6.7 earthquake occurred in 1952 [START_REF] Modisi | Rift kinematics during the incipient stages of continental extension: Evidence from the nascent Okavango rift basin, northwest Botswana[END_REF][START_REF] Reeves | Rifting in the Kalahari?[END_REF]. This basin has been interpreted by some as an incipient rift at the western end of a belt of weak seismicity that extends from the Okavango basin through the Kariba graben to the east and connects to the Rukwa region of the East African Rift [START_REF] Scholz | Evidence for incipient rifting in southern africa[END_REF][START_REF] Modisi | Rift kinematics during the incipient stages of continental extension: Evidence from the nascent Okavango rift basin, northwest Botswana[END_REF], although that interpretation has recently been challenged [START_REF] Pastier | Is the Okavango Delta the terminus of the East African Rift System? Towards a new geodynamic model: Geodetic study and geophysical review[END_REF]. The 2017 Bostwana earthquake is located well away from these structures. Its GCMT mechanism and fault plane solutions derived from Differential In-SAR data show a fault plane striking perpendicular to this possible direction of rift propagation [START_REF] Kolawole | Aeromagnetic, gravity, and Differential Interferometric Synthetic Aperture Radar analyses reveal the causative fault of the 3 april[END_REF][START_REF] Albano | Did anthropogenic activities trigger the 3 april 2017 mw 6.5 botswana earthquake?[END_REF]. However, these published InSARderived fault models did not explore the whole range of possible fault planes and another objective of this study is to re-evaluate the source mechanism and depth of this earthquake.

In the following, we show that the Bostwana earthquake (1) occurred at a depth of 29 km, in the lower crust, on a 73 • or a 17 • -dipping fault, and (2) was preceded by two foreshock sequences. We argue that this earthquake sequence was triggered by a local and transient pulse of elevated pore fluid pressure in a lower crust where viscous deformation should otherwise prevail.

InSAR analysis and Bayesian source parameters estimation

We compute an interferogram from Sentinel-1 acquisitions in interferometric wideswath mode on March 30th and April 11th 2017 along ascending orbits (see Supp. Mat.).

We identify surface displacements resulting from this normal faulting event as the ovalshaped set of fringes located in the epicentral region (Fig. 2). We measure about 4 cm of coseismic surface displacement in the satellite Line-Of-Sight (LOS). Given the incidence angle of the Radar signal in the epicentral region and assuming negligible horizontal motion, this corresponds to about 6 cm of subsidence.

We use a Bayesian formulation to invert for the source parameters of the mainshock. The goal is to find wether it is possible to discriminate between the 2 possible fault planes provided by the GCMT solution. We select a subset of the interferogram and downsample the interferometric phase using a curvature-based quadtree algorithm [Fig 2 and Supp. Mat. Fig. S1;[START_REF] Simons | Coseismic deformation from the 1999 mw 7.1 hector mine, california, earthquake as inferred from insar and gps observations[END_REF]. Our model set up includes a single fault plane with constant slip embedded in an elastic homogeneous medium. We solve for the fault centroid position (longitude, latitude and depth), its orientation (strike and dip), and the amount of slip along dip. We consider its size fixed (10 km along dip and 30 km along strike, chosen based on orders of magnitude for M w 6 earthquakes) as these size parameters trade-off with slip. In addition, we solve for a constant offset to add to the InSAR data.

Following Bayes' theorem, we write the posterior Probability Density Function (PDF)

of the model, P(m|d), as proportional to the product of the prior PDF, P(m), and the likelihood, P(d|m),such as

P(m|d) ∝ P(m)P(d|m), ( 1
)
where m is the vector of model parameters and d the vector of data to invert (i.e.

here, the downsampled surface displacements from InSAR). We use uniform distributions for the prior PDF, imposing positivity on slip (i.e., the fault has to be a normal fault). We use a Gaussian formulation for the likelihood, P(d|m), that writes

P(d|m) ∝ e -1 2 (Gm-d) T C -1 d (Gm-d) , (2)
where, G is the matrix of Green's functions relating source parameters to surface displacements and C d is the data covariance matrix. We build the Green's functions using Okada's formulation of the surface displacements produced by slip on a rectangular dislocation [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF]. We use the covariance function of the interferogram to build the data covariance matrix, describing the influence of turbulent tropospheric noise [see Supp.

Mat.; [START_REF] Jolivet | Aseismic slip and seismogenic coupling along the central San Andreas Fault[END_REF].

In the present case, the posterior PDF must be multimodal as both southwest-and northeast-dipping fault planes should be able to fit surface displacements. Using a classic sampling approach, for instance based on the long-used Metropolis algorithm, will unlikely resolve directly the complete shape of the PDF as the probability for one Markov chain to jump between isolated, high probability regions of the model space is very low. This is likely why a previously published Bayesian model does not show multiple modes for the posterior PDF, hence multiple possible dip angles [START_REF] Albano | Did anthropogenic activities trigger the 3 april 2017 mw 6.5 botswana earthquake?[END_REF]. Even if this probability is not null and sampling all regions of the model space should theoretically be possible with a classic Metropolis algorithm [START_REF] Xu | Identifying active faults by improving earthquake locations with insar data and bayesian estimation: The 2004 Tabuk (Saudi Arabia) earthquake sequence[END_REF], a prohibitive number of steps would be required. Furthermore, it would be difficult to assert the respec-tive importance of each mode. We use instead the AlTar sampler, specifically designed for high dimensional problems and complex PDF sampling (see Supp. Mat.).

Our posterior PDF shows four family of models, including two equally likely and two equally less likely that we separate using a K-mean clustering algorithm (Supp. Fig. S3). The two most likely families of models have an average depth of 29 ± 4 km, an average magnitude of 6.54 ± 0.05 and their centroid location is consistent with that from GCMT. Members from these 2 families only differ by their dip angle, one being on average 17 • ± 4 • and the other at 73 • ± 4 • , hence the two possible families of strike angle at 180 • from eachother. These values are similar to those of the conjugate planes of the GCMT solution, with slightly steeper or shallower dip angles for the steep -and a shallow -angle planes, respectively. The InSAR derived centroid location and magnitude are consistent within uncertainties with those determined by GCMT. We will not further consider the less likely families as their fault strike and magnitude, centroid location, and depth are not consistent with the seismologically-derived ones. Finally, our results indicate that it is not possible to discriminate between a steep and a shallow angle normal faulting event as both families of models are equiprobable.

Aftershocks and foreshocks detection

We then seek to detect aftershocks and possible foreshocks of the Botswana earthquake. Since no data from local seismic networks were available at the time of this work, we apply template matching to continuous signals recorded at teleseismic distances (1200 to 2000 km) to the mainshock from November 2016 to April 2017. This technique has been used to detect low frequency earthquakes (LFE) within tremor signals [START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF] or to recover missing events in aftershock [START_REF] Lengliné | Inferring the coseismic and postseismic stress changes caused by the 2004 mw=6 parkfield earthquake from variations of recur-rence times of microearthquakes[END_REF][START_REF] Peng | Migration of early aftershocks following the 2004 parkfield earthquake[END_REF] and foreshock sequences [START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF][START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Lengliné | Decay and expansion of the early aftershock activity following the 2011, Mw9.0 Tohoku earthquake[END_REF][START_REF] Kato | Multiple slow-slip events during a foreshock sequence of the 2014 Iquique, Chile Mw 8.1 earthquake[END_REF][START_REF] Gardonio | Seismic and aseismic slip: the japanese subduction zone[END_REF]. It is however more challenging at teleseismic distances because of the much lower signal-to-noise ratio of the seismic records compared to near-field observations.

We select as templates 18 M>3 aftershocks that occurred in the week following the main event to compute their coherence with continuous seismic records at five teleseismic stations (Fig. 1). We obtain a continuous record of the coherence per template and per day. We define a coherence threshold above which we consider that an event has been de-tected (see Supp. Mat.). Figure 3 shows the cumulative number of events detected at least at 2 teleseismic stations using coherence thresholds of 0.89, 0.9, and 0.91 after removing auto-detections, i.e. detections of a template by itself. As expected, the number of detections increases with decreasing threshold. All events are seismic in origin but some of the lower-threshold ones are located up to 35 km away from the mainshock epicenter, suggesting a fairly broad seismically active region both before and after the mainshock. We also observe, for all detection thresholds, an increase in seismicity between December 4th and 30th, 2016 -four months before the main event -while only one M4.1 earthquake was reported by the ISC within 120 km of the main shock during that same time interval.

At the 0.91 coherence threshold, we detect a total of 20 new events (see Supp. Mat.).

We manually checked these detections for each station. This number is expectedly much smaller than near-field template matching studies [START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Lengliné | Decay and expansion of the early aftershock activity following the 2011, Mw9.0 Tohoku earthquake[END_REF][START_REF] Kato | Multiple slow-slip events during a foreshock sequence of the 2014 Iquique, Chile Mw 8.1 earthquake[END_REF] but nevertheless shows a remarkable temporal distribution, with two sequences of foreshocks. The first one occurs between 4 and 3 months before the mainshock, with 9 events detected by templates south-west of the mainshock epicenter (Supp. Mat. Fig. S4). Its is followed by a seismically quiet time interval until early

March. The second foreshock sequence occurs during the 2 weeks that precede the main event, with 5 events detected by three templates north-west of the mainshock.

A precursory geodetic signal?

By reference to other preparatory phases identified before large earthquakes, the foreshocks detected here may be embedded in a broader aseismic event [START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 iquique Mw 8.1 earthquake[END_REF][START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF]. In order to detect such an event, we processed all available Sentinel-1 acquisitions between April 2015 and August 2017 (Supp. Mat.). The resulting time series of relative surface displacements between the epicentral region (< 10 km from the epicenter) and the far-field (> 60 km from the epicenter) shows no significant signal within the precursory phase during which small earthquakes are detected (Fig. 3). The same holds for the full time series (Supp. Fig. S8 to S11). The coseismic signal is well recovered, with a 2-3 cm offset at the time of the main shock. Averaging over the last 11 acquisitions of the time series suggests that the coseismic offset is followed by up to 1 cm of post-seismic displacements.

4 Present-day regional strain accumulation?

The location of the Bostwana earthquake far from areas of concentrated seismic activity, in a region with no morphological evidence of recent tectonic activity, raises the question of present-day strain accumulation in southern Africa. The low strain rates expected in such an intraplate setting may be difficult to measure geodetically, especially since the distribution of GPS stations is sparse and uneven. Nevertheless, it is useful to try place an upper-bound on regional strain accumulation using the existing permanent GPS stations. We therefore updated the analysis of data from openly available, continuouslyrecording GPS stations in southern Africa to derive a continental-scale velocity field in order to search for region-wide deformation (Supp. Mat. and Fig. 1).

We search for deviations from a purely rigid behavior by estimating a single rigid rotation for the whole region considered here and examining residual velocities. We find no spatial pattern in residual velocities with respect to a rigid plate (Fig. 1) with an RMS misfit of 0.25 mm/yr (maximum residual of 0.95 mm/yr), and a χ 2 of 1.1. The comparison of the distribution of residual velocities, normalized by their uncertainty, with that of a two-dimensional, unit variance, normal distribution shows that the velocities are well described by a random process (Fig. 1). Although the uneven geographic distribution of GPS stations in southern Africa is not optimal for this type of study, our results rule out the possibility of broad-scale deformation at a level of 0.25 mm/yr, consistent with observations in several other plate interiors [START_REF] Nocquet | Present-day kinematics of the Mediterranean: A comprehensive overview of GPS results[END_REF][START_REF] Tregoning | A decade of horizontal deformation from great earthquakes[END_REF][START_REF] Craig | Strain accumulation in the New Madrid and Wabash Valley seismic zones from 14 years of continuous gps observation[END_REF], as well as previous results in the same region [START_REF] Hackl | Estimation of velocity uncertainties from GPS time series: Examples from the analysis of the South African TrigNet network[END_REF][START_REF] Saria | Present-day kinematics of the East African Rift[END_REF].

What caused the Botswana earthquake?

The 2017, M w 6.5, Botswana earthquake occurred in an area with no previous evidence of similar magnitude events, low and diffuse background seismicity, and no topographic features indicative of repeated recent faulting, characteristics that are shared by most large earthquakes in stable continental regions [START_REF] Calais | A new paradigm for large earthquakes in stable continental plate interiors[END_REF]. That the earthquake focal mechanism shows purely normal faulting is consistent with the occurrence of other normal faulting earthquakes within stable southern Africa [START_REF] Heidbach | World stress map database release 2016[END_REF] and with stress models derived from horizontal gradients of gravitational potential energy [START_REF] Coblentz | Tectonic stresses in the African plate: Constraints on the ambient lithospheric stress state[END_REF][START_REF] Stamps | Current kinematics and dynamics of Africa and the East African Rift System[END_REF]. Both indicate an extensional stress state with a subvertical maximum principal compressive stress (σ 1 ) and a subhorizontal least principal compressive stress (σ 3 ).

However, a striking feature of the event is its depth, which all authors consistently find between 25 and 30 km, well into the lower part of a crust that is around 35 km-thick in the epicentral area [START_REF] Nguuri | Crustal structure beneath southern Africa and its implications for the formation and evolution of the Kaapvaal and Zimbabwe cratons[END_REF][START_REF] Tedla | A crustal thickness map of Africa derived from a global gravity field model using Euler deconvolution[END_REF][START_REF] Youssof | Moho depth and crustal composition in Southern Africa[END_REF]. The occurence of earthquakes in the lower crust is often interpreted as evidence for a mafic composition [START_REF] Shudofsky | Unusually deep earthquakes in East Africa: Constraints on the thermo-mechanical structure of a continental rift system[END_REF][START_REF] Nyblade | East African earthquakes below 20 km depth and their implications for crustal structure[END_REF][START_REF] Albaric | Crustal rheology and depth distribution of earthquakes: Insights from the central and southern East African Rift System[END_REF][START_REF] Craig | Earthquake distribution patterns in Africa: their relationship to variations in lithospheric and geological structure, and their rheological implications[END_REF]. Seismic data in the Kaapval craton however show a low Poisson ratio of 0.25 for the whole crust and a 2,860 kg/m 3 lowermost crust density, indicative of a felsic composition [START_REF] James | Crustal structure of the Kaapvaal craton and its significance for early crustal evolution[END_REF], consistent with the lack of mafic granulite xenoliths from the lower Kaapvaal crust [START_REF] Schmitz | Ultrahigh-temperature metamorphism in the lower crust during Neoarchean Ventersdorp rifting and magmatism, Kaapvaal Craton, southern Africa[END_REF].

Given a crustal geotherm derived from local surface heat flow measurements, crustal yield strength enveloppes for the region (Fig. 4a) show that brittle failure should not happen at such depth, under hydrostatic pore-fluid pressure, as a felsic lower crust flows at low differential stress. In the absence of a mafic lower crust, (1) the maximum differential stress that can be maintained at hypocentral depth is about 50-100 MPa and (2) deformation at such depth is controlled by viscous flow. However, if pore-fluid pressure becomes sub-lithostatic, brittle failure on a 73 • -dip fault is allowed at a differential stress lower than 50 MPa (Fig. 4c). Alternate mechanisms are also possible, such as thermal shear runaway, rupture of a brittle asperity, or dehydration reactions [e.g. [START_REF] Green | The Mechanics of Deep Earthquakes[END_REF][START_REF] Prieto | Seismic evidence for thermal runaway during intermediate-depth earthquake rupture[END_REF]. The former two require significant, on-going, shear motion, which is unlikely in this stable cratonic environment. The latter implies a phase transition, which would require a recent rejuvenation of the cratonic crust that is not documented. The stability of the cratonic crust requires an external forcing to trigger this event.

Therefore, a likely explanation for the Botswana earthquake is that it was triggered by elevated, sub-lithostatic, pore fluid pressure that enabled failure at the low differential stress that prevails in the viscous lower crust [START_REF] Gold | Fluid ascent through the solid lithosphere and its relation to earthquakes[END_REF]. The observed foreshock swarm-like sequences may be the signature of the initiation of a pulse of high pore fluid pressure [START_REF] Hainzl | Indications for a successively triggered rupture growth underlying the 2000 earthquake swarm in Vogtland/NW Bohemia[END_REF][START_REF] Reyners | The role of fluids in lower-crustal earthquakes near continental rifts[END_REF][START_REF] Balfour | Localization of intraplate deformation through fluid-assisted faulting in the lower-crust: The Flinders Ranges[END_REF].

Field observations show numerous evidence of fluid-assisted embrittlement in the viscous regime of deformation [START_REF] Handy | Continental fault structure and rheology from the frictional-to-viscous transition downward[END_REF][START_REF] Wehrens | Deformation at the frictional-viscous transition: Evidence for cycles of fluid-assisted embrittlement and ductile deformation in the granitoid crust[END_REF]. Lower-crustal earthquakes in the northern Alpine foreland [START_REF] Deichmann | Structural and rheological implications of lower-crustal earthquakes below northern Switzerland[END_REF] and beneath the Flinders Ranges of South Australia [START_REF] Balfour | Localization of intraplate deformation through fluid-assisted faulting in the lower-crust: The Flinders Ranges[END_REF] have been interpreted as the result of a decrease of effective stress on pre-existing faults by fluids at near-lithostatic pore pressure, allowing a switch from viscous to brittle deformation. In the later case, the authors argued for a deep fluid source from a remnant hydrated mantle on the basis of elevated 3 He/ 4 He ratios in springs, as also observed in Vogtland, Bohemia, and Eger Rift intraplate seismicity areas of Central Europe [START_REF] Weise | Transport of mantle volatiles through the crust traced by seismically released fluids: a natural experiment in the earthquake swarm area Vogtland/NW Bohemia[END_REF][START_REF] Bräuer | Earthquake swarms in non-volcanic regions: What fluids have to say[END_REF]. In the Taupo active rift (New Zealand), lower-crustal earthquakes in the viscous regime are interpreted as triggered by fluids migrating upward from the hydrated Hikurangi subduction mantle wedge [START_REF] Reyners | The role of fluids in lower-crustal earthquakes near continental rifts[END_REF].

Southern Africa is a largely cratonic province characterized by widespread kimberlite outcrops, which take their source in carbonate-rich matrices and parental magma [START_REF] Kamenetsky | Towards a new model for kimberlite petrogenesis: Evidence from unaltered kimberlites and mantle minerals[END_REF]. Rapid kimberlite melt ascent through the crust is assumed to be driven by the exsolution of a H 2 O-and CO 2 -rich fluid phases at mantle depths [START_REF] Russell | Kimberlite ascent by assimilation-fuelled buoyancy[END_REF]. This requires the presence of volatiles in the mantle, which can be hosted in the lower lithosphere until remobilization under large-scale tensional tectonic stresses as shown in the Virunga volcanic field of the East African Rift [START_REF] Hudgins | Melt inclusion evidence for CO2-rich melts beneath the western branch of the East African Rift: implications for long-term storage of volatiles in the deep lithospheric mantle[END_REF]. In the Okavango basin, 350 km north of the Botswana earthquake, a thermal anomaly measured in the absence of surface magmatism and of a thinned or altered lithosphere is interpreted as the signature of fluids advected from a metasomatized lithospheric mantle [START_REF] Leseane | Thermal perturbations beneath the incipient Okavango Rift Zone, northwest Botswana[END_REF]. Therefore, several lines of evidence point to the presence of significant amounts of fluids in the upper mantle underneath southern and eastern Africa. The upward migration of these deep and buoyant fluids could perhaps explain the locally elevated pore-fluid pressure necessary to trigger seismicity in the otherwise ductile lower crust of the Kaapval craton.

Conclusion

The occurrence of the April 3rd, 2017, Botswana earthquake in a felsic lower crust where viscous deformation should prevail indicates that pore-fluid pressure elevated to sub-lithostatic played a key role in triggering the rupture. The two swarm-like sequences of earthquakes that preceded the main Botswana event in December and March 2017 may be further evidence for fluid movement in a critically loaded fault network, that eventually led to a large event. Finally, the damage caused by the mainshock potentially led to a de- Normalized Detection 2 0 1 6 -1 1 -0 1 2 0 1 6 -1 2 -0 1 2 0 1 7 -0 1 -0 1 2 0 1 7 -0 2 -0 1 2 0 1 7 -0 3 -0 1 2 0 1 7 -0 4 -0 1

Figure 1 .Figure 2 .

 12 Figure 1. Regional context -Map of the southern tip of the African continent. The NEIC epicenter of the April 3, 2017, Bostwana earthquake is shown with a yellow circle, with its GCMT focal mechanism. Black triangles indicate the location of seismological stations used for template matching. Grey dots are earthquake epicenters from the NEIC catalog covering the 1974-2018 time interval. Dotted lines show the boundaries of the Kaapval and Zimbabwe archean cratons. Solid black lines show major faults associated with recent topographic features. Red arrows show residual GPS velocities with respect a rigid plate assumption with their 95% confidence ellipse. Uncertainties include site-dependent, time-correlated colored noise. Red labels indicate the code name of some of the GPS stations. Lower right inset shows the normalized cumulative distribution of GPS velocity residuals (both north and east components). The red line shows the theoretical χ 2 distribution expected if residuals are normally distributed in two dimensions with a unit variance.

  Figure 1.

  Figure 4.

Table 1 .

 1 Catalog of the newly detected events. The time of occurrence corresponds to the time at which coherency between the template and the continuous signal is the highest. The locations of the detected events are the same as the location of the template they matched, as indicated in the ISC catalog. Magnitudes are estimated by calculating the amplitude ratio (see Supp. Mat.).

	491
	492
	493
	494
	Year Month Day Hour Minute Second Latitude Longitude Magnitude
	2016	12	11	17	9	32.05	-22.5367	24.9746	3.3
	2016	12	11	19	28	14.05	-22.6610	25.0010	3.5
	2016	12	13	21	46	46.05	-22.5367	24.9746	3.3
	2016	12	14	17	22	52.05	-22.6610	25.0010	3.6
	2016	12	27	12	1	44.05	-22.6913	25.1021	3.6
	2016	12	27	16	16	18.05	-22.6610	25.0010	3.8
	2016	12	28	10	56	54.05	-22.9860	25.1260	3.7
	2016	12	28	15	12	22.05	-22.6610	25.0010	3.4
	2016	12	29	7	39	24.05	-22.9870	24.9980	3.7
	2017	03	14	18	38	20.05	-22.5646	25.0868	3.1
	2017	03	15	15	12	2.050	-22.6610	25.0010	3.9
	2017	03	23	17	43	12.05	-22.5367	24.9746	3.5
	2017	03	27	10	18	4.050	-22.6610	25.0010	3.9
	2017	04	02	8	30	40.05	-22.5367	24.9746	3.5
	2017	04	05	12	29	10.05	-22.3206	25.4211	3.4
	2017	04	11	11	57	32.05	-22.5367	24.9746	3.5
	2017	04	12	19	3	28.05	-22.8180	24.9340	3.8
	2017	04	21	11	52	4.050	-22.5367	24.9746	3.5
	2017	04	21	21	36	48.05	-22.6610	25.0010	3.6
	2017	04	27	5	36	56.05	-22.6784	25.1558	3.3

Acknowledgments

We thank IRIS for granting access to seismic data (http://ds.iris.edu/SeismiQuery). Sentinel 1A and B acquisitions were downloaded through the Plateforme d'Exploitation des Produits Sentinel (PEPS, https://peps.cnes.fr) of the Centre National des Études Spatiales (CNES). The seismic catalog was downloaded through the ISC website (http://www.isc.ac.uk/).

Global Precipitation

Measurement data are available for download at https://pmm.nasa.gov/GPM. BG received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (REALISM project, grant agreement 681346). RJ received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Geo-4D project, grant agreement 758210). EC acknowledges support from the Institut Universitaire de France and from the INSU/CNRS Tellus-Rift program. We thank K. Materna and an anonymous reviewer for the time they took at evaluating our work and their constructive comments.

Confidential manuscript submitted to GRL crease in pore-fluid pressure locally, turning off the activity of this swarm-like sequence, hence the detection of a classic Omori decay of aftershock productivity.

The Botswana earthquake therefore did not require localized, present-day, stress or strain accumulation, contrary to plate boundary events resulting from the near-fault accrual of stress imposed by plate and block motions [START_REF] Kanamori | The physics of earthquakes[END_REF]. In a pre-stressed crust able to store reversible strain on long timescales [START_REF] Feldl | Great Himalayan earthquakes and the Tibetan plateau[END_REF][START_REF] Craig | Evidence for the release of long-term tectonic strain stored in continental interiors through intraplate earthquakes[END_REF] with faults at or close to failure [START_REF] Townend | How faulting keeps the crust strong[END_REF], short-term fault strength transients, such as those triggered by fluids leaks from the upper mantle, may be all it takes to trigger large events. 2 0 1 6 -1 1 -0 1 2 0 1 6 -1 2 -0 1 2 0 1 7 -0 1 -0 1 2 0 1 7 -0 2 -0 1 2 0 1 7 -0 3 -0 1 2 0 1 7 -0 4 -0 [START_REF] Rutter | Experimental intracrystalline plastic flow in hotpressed synthetic quartzite prepared from Brazilian quartz crystals[END_REF] and for dry quartzite from [START_REF] Hirth | An evaluation of quartzite flow laws based on comparisons between experimentally and naturally deformed rocks[END_REF]. Strain rate is 10 × 10 -18 s-1. We use a surface heat flow of 44 mW/m 2 , the average of four close-by measurements [START_REF] Ballard | Terrestrial heat flow in Botswana and Namibia[END_REF], and computed the corresponding crustal geotherm [START_REF] Russel | A steady state conductive geotherm for the north central Slave, Canada: Inversion of petrological data from the Jericho Kimberlite pipe[END_REF]. Blue lines show friction law for hydrostatic (λ V = 0.37) and sub-lithostatic (λ V = 0.85) pore fluid pressure [START_REF] Byerlee | Friction of rocks[END_REF]. B Geotherms derived from a surface heat flow of 44 ±10% mW/m 2 calculated following [START_REF] Russel | A steady state conductive geotherm for the north central Slave, Canada: Inversion of petrological data from the Jericho Kimberlite pipe[END_REF]. C Differential stress (σ 1 -σ 3 ) required for frictional reactivation of cohesionless normal faults at 30 km depth as a function of their dip angle for different values of the pore fluid factor λ V , following [START_REF] Sibson | High-angle reverse faulting in northern New Brunswick, Canada, and its implications for fluid pressure levels[END_REF]. Lithostatic conditions correspond to λ V = 1. Histograms are the probability density functions of the earthquake source parameters estimated from InSAR data, including centroid depth and dip angle.


Pierre Romanet
email: pierreromanet@geologie.ens.fr

Harsha S Bhat

Romain Jolivet

Raúl Madariaga

Fast and slow slip events emerge due to fault geometrical complexity

Fault geometry can be a natural source of slip complexity in earthquake cycle modeling, resulting in slow slip events (SSE) and earthquakes.

• A simple two overlapping fault model produces different observed scaling laws for earthquakes and for slow slip events.

• All observed complexities emerge with uniform loading and rate weakening friction properties on the fault.

Introduction

Since their discovery in the late nineties, Slow-Slip Events (SSE) have been widely observed along various subduction zones (Central Ecuador [START_REF] Vallee | Intense interface seismicity triggered by a shallow slow slip event in the central ecuador subduction zone[END_REF], Southwest Japan [START_REF] Hirose | A slow thrust slip event following the two 1996 hyuganada earthquakes beneath the bungo channel, southwest japan[END_REF], Guerrero [START_REF] Lowry | Transient fault slip in guerrero, southern mexico[END_REF], Cascadia [START_REF] Dragert | A silent slip event on the deeper cascadia subduction interface[END_REF][START_REF] Rogers | Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip[END_REF], Hikurangi [START_REF] Douglas | Slow slip on the northern hikurangi subduction interface, new zealand[END_REF], Northern Chile [START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 iquique mw 8.1 earthquake[END_REF] and others). The discovery of SSEs mainly came from the development and the installation of networks of permanent GPS stations around subduction zones. Although GPS is still nowadays the main SSE detection tool, new observations now allow for the detection of slow-slip, like InSAR [START_REF] Rousset | An aseismic slip transient on the north anatolian fault[END_REF][START_REF] Jolivet | Spatio-temporal evolution of aseismic slip along the haiyuan fault, china: Implications for fault frictional properties[END_REF], networks of sea-bottom pressure gauge [START_REF] Ito | Episodic slow slip events in the japan subduction zone before the 2011 tohoku-oki earthquake[END_REF][START_REF] Wallace | Slow slip near the trench at the hikurangi subduction zone, new zealand[END_REF] or, indirectly, via the migration of microseismicity, repeating earthquakes and tremors [START_REF] Igarashi | Repeating earthquakes and interplate aseismic slip in the northeastern japan subduction zone[END_REF][START_REF] Kato | Propagation of Slow Slip Leading Up to the 2011 Mw 9.0[END_REF], thus increasing significantly the probability of their detection. SSEs, like earthquakes, correspond to an accelerating slip front propagating along a fault. However, unlike earthquakes, SSEs themselves do not radiate any detectable seismic waves and are hence sometimes nicknamed "silent events". Until the discovery of SSEs, it was thought that only earthquakes release the accumulated strain energy along a fault. Since SSEs also contribute to this release of energy, they should play an important role in the estimation of seismic hazard along subduction zones [START_REF] Obara | Connecting slow earthquakes to huge earthquakes[END_REF]. In addition, SSEs exhibit very specific characteristics. Their propagation speed -2-Confidential manuscript submitted to Geophysical Research Letters along the fault (about 0.5 km/h in Cascadia [START_REF] Dragert | Geodetic and seismic signatures of episodic tremor and slip in the northern cascadia subduction zone[END_REF] to about 1 km/day in Mexico [START_REF] Franco | Propagation of the 2001-2002 silent earthquake and interplate coupling in the oaxaca subduction zone, mexico[END_REF]) contrasts with the rupture propagation speed of earthquakes (at about 3 km/s). The slip velocity of SSEs (from about 1mm/yr in the Bungo Channel, Japan to about 1 m/year in Cascadia) is around one or two orders of magnitude greater than plate convergence rates but orders of magnitude smaller than earthquakes slip rates (of the order of 1m/s) [START_REF] Schwartz | Slow slip events and seismic tremor at circum-pacific subduction zones[END_REF].

Although the exact influence of SSEs in the seismic cycle is not yet fully understood, they seem closely related to earthquakes. Several seismic and geodetic observations suggest that SSEs may have happened just before and in regions overlapping with earthquakes. The 2011 M w 9.0 Tohoku-Oki event and the 2014 M w 8.1 Iquique event are two examples in subduction zones where a SSE apparently occurred just before the earthquake, within a region overlapping with the area where seismic slip nucleated [START_REF] Kato | Propagation of Slow Slip Leading Up to the 2011 Mw 9.0[END_REF][START_REF] Brodsky | Recognizing Foreshocks from the 1 April 2014 Chile Earthquake[END_REF][START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 iquique mw 8.1 earthquake[END_REF][START_REF] Mavrommatis | Long-term acceleration of aseismic slip preceding the Mw 9 Tohoku-oki earthquake: Constraints from repeating earthquakes[END_REF]. More recently, geodetic evidence of a large SSE triggering an earthquake was pointed out in the Guerrero subduction zone [START_REF] Radiguet | Triggering of the 2014 mw7.3 papanoa earthquake by a slow slip event in guerrero, mexico[END_REF]. There are also suggestions that SSEs may be triggered by earthquakes either by stress-waves and/or static stress transfer [START_REF] Itaba | A slow slip event triggered by teleseismic surface waves[END_REF]Zigone et al., 2012;[START_REF] Kato | Detection of a hidden boso slow slip event immediately after the 2011 mw 9.0 tohoku-oki earthquake, japan[END_REF][START_REF] Wallace | Large-scale dynamic triggering of shallow slow slip enhanced by overlying sedimentary wedge[END_REF]. On the other hand some SSEs occur without an accompanying large earthquake as in the Cascadia subduction zone, where SSEs occur periodically [START_REF] Rogers | Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip[END_REF], or in the Hikurangi subduction zone [START_REF] Wallace | Slow slip near the trench at the hikurangi subduction zone, new zealand[END_REF]. From the above examples, it seems that there may or may not be a connection between slow slip events and fast earthquakes. Some authors [Obara and Kato, 2016, for e.g.] have suggested that slow slip events, because of their sensitivity to very small stress perturbations, can act as a stress meter of the current stress in the crust. However, this still needs to be confirmed. Also, the exact role of SSE's in hazard assessment remains largely unknown. All SSEs have the same direction of slip as earthquakes, i.e. opposite to the plate convergence direction, and are accompanied by a positive stress drop which corresponds to a reduction in the accumulated strain energy.

In the absence of external forcing mechanism, this necessitates SSEs to occur in a slip, or slip rate, weakening region which is also prone to rupture as a fast dynamic event. These observations, put together, raise the first question. What physical mechanism explains slowslip and fast, dynamic earthquakes occurring under similar frictional boundary conditions along active faults? Our key finding is that fault geometrical complexity gives rise to the variety of modes of slip along an active fault without any other complex mechanism involved.

Furthermore, earthquakes and SSEs seem to follow different scaling laws [START_REF] Ide | A scaling law for slow earthquakes[END_REF], which remain out of reach of numerical models until now [START_REF] Ide | Modeling fast and slow earthquakes at various scales[END_REF]. The seismic moment of earthquakes scales with the cube of their duration (M ∝ T 3 ) whereas the corresponding moment of SSEs is proportional to their duration (M ∝ T), raising the second question. Is such different scaling a general feature of earthquakes and SSEs, highlighting different physical mechanisms [START_REF] Ide | Bridging the gap between seismically and geodetically detected slow earthquakes[END_REF][START_REF] Peng | An integrated perspective of the continuum between earthquakes and slow-slip phenomena[END_REF][START_REF] Ide | Modeling fast and slow earthquakes at various scales[END_REF]? We address the above questions using physics-based numerical modeling of active faults governed by rate-and-state friction [START_REF] Dieterich | Time-dependent friction and the mechanics of stick-slip[END_REF] and develop a unified framework that addresses all the observations about earthquakes and SSEs mentioned above.

Modeling slow, aseismic slip

SSEs were discovered to emerge spontaneously from numerical models in the rateand-state framework for the modeling of subduction zones [Liu andRice, 2005, 2007].

In this framework, fault areas with weakening properties will preferentially host seismic slip (i.e. earthquakes) while strengthening regions will host stable continuous creep or post-seismic slip. Numerical experiments and theoretical works have shown that the main physical control on the emergence of SSEs in models is how the characteristic length of a weakening patch [START_REF] Ruina | Slip instability and state variable friction laws[END_REF][START_REF] Rice | Constitutive relations for fault slip and earthquake instabilities[END_REF][START_REF] Dieterich | Earthquake nucleation on faults with rate-and state-dependent strength[END_REF][START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF] compares to the specific nucleation length scale [START_REF] Liu | Aseismic slip transients emerge spontaneously in 3d rate and state modeling of subduction earthquake sequences[END_REF][START_REF] Rubin | Episodic slow slip events and rate-and-state friction[END_REF]. If the length of a fault patch is large compared to the nucleation length scale, earthquakes have enough room to grow and become dynamic, so this fault patch will generate only dynamic, seismic events. If the length of the fault is small compared to this length scale, earthquakes can never grow large enough to become dynamic or no events will occur at all (i.e. permanent creep). It is therefore necessary, under this framework, to tune for the right fault length compared to the nucleation length scale to allow modeling of both slow and fast ruptures. Given the observed spatial size over which some SSEs propagate i.e. on the order of tens of kilometers, this would lead to unrealistically large nucleation sizes, preventing the occurrence of any earthquakes. A possible explanation for such large nucleation lengths could be the presence of high-pressure pore fluids released during metamorphic dehydration reactions. However it has been shown recently that regions of high fluid pressure and slow slip events do not always overlap along all the subduction zones [Saffer -4-Confidential manuscript submitted to Geophysical Research Letters [START_REF] Saffer | The frictional, hydrologic, metamorphic and thermal habitat of shallow slow earthquakes[END_REF]. One solution to overcome this issue is to appeal to other competing frictional mechanisms like dilatant-strengthening [START_REF] Segall | Dilatancy, compaction, and slip instability of a fluidinfiltrated fault[END_REF][START_REF] Rubin | Episodic slow slip events and rate-and-state friction[END_REF][START_REF] Segall | Dilatant strengthening as a mechanism for slow slip events[END_REF] with or without thermal-pressurization [START_REF] Segall | Slow-slip evolves into megathrust earthquakes in 2d numerical simulations[END_REF]. Although we do not include these additional frictional mechanisms in our modeling below, we acknowledge that it would broaden the range over which we are able to observe slowslip.

As the above models suggest, a set of competing mechanisms are required for slowslip and earthquakes to coexist. However, there is one ubiquitous feature that is often ignored for computational reasons: the geometric complexity of active faults. Indeed, faults are rarely planar over length scales of tens of kilometers and in fact, fault segmentation and geometric complexity are visible at multiple scales [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. Subduction zones also show geometrical complexities like subducting seamounts [START_REF] Das | Effect of subducting seafloor topography on the rupture characteristics of great subduction zone earthquakes[END_REF]. It is also known that subduction zones have large normal faults that connect the main slab and can sometimes be reactivated during seismic events [START_REF] Hicks | Seismic slip on an upper-plate normal fault during a large subduction megathrust rupture[END_REF][START_REF] Hubbard | Coseismic slip on shallow décollement megathrusts: implications for seismic and tsunami hazard[END_REF].

This non-planarity of faults should introduce a natural stress based interaction between faults. Several lines of evidence suggest that geometric complexity should be considered in conjunction with various observed slip dynamics. Aseismic slip has been observed with earthquake swarms in the northern Apennines (Italy) along splay faults [START_REF] Gualandi | Aseismic deformation associated with an earthquake swarm in the northern apennines (italy)[END_REF]. It has been detected along the Haiyuan fault (China) [START_REF] Jolivet | Spatio-temporal evolution of aseismic slip along the haiyuan fault, china: Implications for fault frictional properties[END_REF], the North Anatolian Fault [START_REF] Rousset | An aseismic slip transient on the north anatolian fault[END_REF][START_REF] Bilham | Surface creep on the north anatolian fault at ismetpasa, turkey, 1944-2016[END_REF] and, in earlier publications, along the San Andreas Fault [START_REF] Murray | Spatiotemporal evolution of a transient slip event on the san andreas fault near parkfield, california[END_REF]]. SSE's have been observed in the very shallow part of subduction zones, such as in Hikurangi [START_REF] Wallace | Slow slip near the trench at the hikurangi subduction zone, new zealand[END_REF] and Nankai [START_REF] Araki | Recurring and triggered slow-slip events near the trench at the nankai trough subduction megathrust[END_REF], among others. The only known common ingredient of all of these different seismotectonic settings is the geometrical complexity of faults across scales.

In this work, we have restricted ourself to only one type of geometric complexity i. e. two overlapping faults. Of course, this geometry cannot be interpreted directly as a subduction zone or any other natural setting. However, we suggest that if this simple geometry can give rise to a complex slip behaviour in the seismic cycle then a more realistic description of fault zones with multiple slip surfaces should not be ignored.

-5-Confidential manuscript submitted to Geophysical Research Letters

Model set-up

Our aim is to test the influence of fault geometry on the behavior of slip along a fault. We build a conceptual model in which fault slip is controlled by an unstable frictional rheology (rate weakening) without any lateral variation. Doing so, we introduce no a priori complexity in initial and boundary conditions. We load the faults with constant stress loading rate and observe the variety of modes of slip.

In our conceptual model, we consider two overlapping faults of the same length L (see geometry in Fig. 1). This geometry is chosen to illustrate the effect of complex stress interactions between neighboring faults or fault segments and is in no way supposed to be interpreted as the only geometrical configuration of faults in a fault network. Friction on both faults is controlled by rate-and-state friction with aging state evolution. Frictional resistance decreases with increasing slip rate and is spatially uniform, i.e. the fault is rateweakening. Loading is imposed using a constant rate of shear stress increase on the fault.

We model elastic interactions using out-of-plane static stress interactions with a radiation damping approximation [START_REF] Rice | Spatio-temporal complexity of slip on a fault[END_REF]. The computation of static stress interactions is accelerated using the Fast Multipole Method, allowing us to compute all stages of the earthquake cycle in a tractable computational time [START_REF] Greengard | A fast algorithm for particle simulations[END_REF][START_REF] Carrier | A fast adaptive multipole algorithm for particle simulations[END_REF] (See Methods section for more details).

To better understand the role of multi-fault interactions, we explore the influence of the distance between faults, D, the length of the faults, L, and the ratio of the constitutive frictional parameters, a/b. For rate-weakening faults, a/b ranges between 0 and 1.

Because of the importance of the nucleation length scale L nuc in this problem, all geometrical parameter are non-dimensionalized by L nuc ,

L nuc = 2 π µD c σ n b(1 -a/b) 2 ; a/b → 1 (1)
where, a and b are rate-and-state constitutive friction parameters, D c is the characteristic slip distance, µ is the shear modulus of the medium and σ n the normal stress acting on the fault [START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF][START_REF] Viesca | Stable and unstable development of an interfacial sliding instability[END_REF]. This formulation provides good insights on the nucleation phase of earthquakes along a fault that is mildly rate-weakening (a/b → 1).

For computational reasons, we restrict our experiments to fault lengths L/L nuc ∈ {1, 2, 3, 4}. Our parameter space includes also distances between faults D/L nuc ∈ {0.1, 0.

Results

For each of the parameters identified above, we initiate the model, and compute slip velocity over time (Fig. 1). We observe cycles of quiescence and earthquakes as expected for a rate-weakening rheology but, unlike in a model with a single, flat fault with no geometrical complexity, we also observe episodes during which slip is slow. In our conceptual model, we see regular earthquakes with a clear nucleation, dynamic and afterslip phases and these events happen without any evident periodicity. We observe what would be considered in nature as the slow nucleation of earthquakes, the slow phase of recovery following an earthquake, earthquakes of variable slip duration and velocity and slow slip events. It appears then, that the sole introduction of a simple geometrical complexity leads to the emergence of the complete range of modes of slip, even with a uniform rate-weakening rheology. Slow-slip events emerge spontaneously without prescribing the necessary conditions for slow slip. In our model, a fault that slipped seismically can also potentially host slow slip, as in the region of overlap of co-and post-seismic slip or along the shallow portion of a creeping fault [START_REF] Wallace | Slow slip near the trench at the hikurangi subduction zone, new zealand[END_REF][START_REF] Rousset | An aseismic slip transient on the north anatolian fault[END_REF]. Once again, without the introduction of a second fault, and its associated stress perturbations, the fault behaves like a simple spring-slider system with weakening properties, with similar earthquakes happening periodically (see Figure S3 in Supp. Mat.).

We believe the choice of such geometry brings realistic perturbations in stress along the fault and these perturbations lead to the emergence of the observed variety of modes of slip. Fig. 1 illustrates the complexity that emerges by only appealing to stress perturbations from a neighboring fault and/or non planarity of the fault. Now considering that faults are geometrically complex at all scales, it appears natural to extend this conclusion and consider that the whole range of modes of slip observed in nature may result, among other mechanisms, from these geometrically-induced stress complexities. In addition, it may be safe to think that models that do not include such complexities will require ad-hoc -7-Confidential manuscript submitted to Geophysical Research Letters tuning, which might not be necessary, to reproduce observations. We have not yet identified the precise conditions leading to an earthquake or a slow slip event, but clues should be found in the analysis of the evolution of stresses and state variable along the fault.

A phase diagram of slip

We allow our model to undergo multiple earthquake cycles before measuring slip and rupture velocity of each slow and dynamic event. We identify SSEs and earthquakes based on their slip and rupture velocity. SSEs are events with a slip velocity V in the range of 1µm/s to 1 mm/s and a rupture velocity V rup lower than 0.001c s , where c s is the shear wave speed. Earthquakes are events with a slip velocity greater than 1 mm/s and a rupture velocity greater than 0.001c s . We also define nucleation as the moment before an earthquake, where slip velocity is higher than 1µm/s until it reached 1 mm/s. We purposefully chose a relatively small threshold value for rupture velocity, because quasi-dynamic simulations lead to much slower rupture velocity than dynamic simulations [START_REF] Thomas | Tremor-tide correlations and near-lithostatic pore pressure on the deep San Andreas fault[END_REF]. As our faults are one dimensional, we define the equivalent moment for a seismic or aseismic event as M = µ DL rup × 1km, where L rup is the total length of the fault that slipped during an event (SSE or earthquake) and D is the slip averaged over the length L rup . For earthquakes, we compute separately the seismic moment during the nucleation phase and the dynamic phase. For SSEs, moment accounts for the entire duration when the slip velocity exceeds 1 µm/s. We obtained about 3000 individual earthquakes and about 500 SSEs in our calculations when the faults hosted both earthquakes and SSEs.

We identify five different domains of fault slip behavior (Fig. 2). For small faults (L << L nuc ), there is a damped domain in which the fault experiences no events at all as the fault length is too small for any type of instability to grow. For long faults (L >> L nuc ) with strongly rate-weakening properties (a/b < 0.5), we observe periodic earthquakes, similar as in a case with no geometric complexity. This is perfectly normal as both our faults are flat and the longer they are, the larger the portion that is left unaffected by the geometrical complexity (i.e. if the faults are long, their edges are independent and dominate the general behavior of slip, reducing this setting to a case with no geometrical complexity). For mildly rate-weakening faults (1 > a/b > 0.6) and whatever the length of the fault, we observe a complex behavior with a mixture of slow and rapid slip for fault sizes between 1 and 4 times the nucleation length and only complex earthquakes (partial ruptures, aperiodic events, variable after slip) for longer faults. That -8-Confidential manuscript submitted to Geophysical Research Letters is, although the length over which we observe slow slip events is increased compared to the case where there is no additional fault, we are still limited by the nucleation length scale. Therefore like in other studies, we will require another mechanism. This can just be low effective normal stress, additional frictional mechanisms like dilatant strengthening or even stronger geometrical complexities. The domain where both slow and fast earthquake coexist, shrinks when the distance between the faults is increased. All this put together confirms our intuition that stress perturbations from one fault to another help modulate the mode of slip along faults.

Scaling

Geodetic and seismological observations in nature suggest two different scaling relationships for moment of slow slip on one side and rapid, dynamic slip events on the other side [START_REF] Ide | A scaling law for slow earthquakes[END_REF][START_REF] Peng | An integrated perspective of the continuum between earthquakes and slow-slip phenomena[END_REF]. Considering the statistics of slip events produced by our model, we also find that the moment of both seismic and aseismic events modeled by rate and state friction law follows two different scaling laws as observed in nature (Fig. 3). Because we conducted our calculations in 2D, the moment of a dynamic slip event should scale with its duration squared: M ∝ T 2 . This scaling emerges naturally from our conceptual model without imposing any complexity in the spatial variation of frictional properties. If we do not include any geometrical complexity, periodic, identical earthquakes are observed impeding our ability to observe any potential scaling. Although we do not preclude the possibility that other models, that have produced SSE's and earthquakes, also reproduce such scaling laws, geometrical complexities give rise to a wide range of modes of slip and the resulting events obey similar scaling laws as in nature.

We note the moment of our simulated events clearly depends on the ratio of constitutive parameters a/b. Since the nucleation length L nuc increases with a/b and since we compare models with non-dimensionalised fault length, the real length of the fault, L, also increases when a/b→1, leading to bigger moment release and longer duration for events.

To verify the robustness of this scaling law, we changed the maximum slip velocity criteria used to distinguish SSEs and earthquakes by one order of magnitude. This does not change the observed scaling.

Another interesting feature that emerges from our calculations is that the moment of the nucleation phase of earthquakes also follows the same linear scaling with duration -9-Confidential manuscript submitted to Geophysical Research Letters as slow-slip events. However, we cannot argue that this similarity in scaling may be preserved in 3D. We finally notice that by adding the nucleation and after-slip moment of earthquakes, the clear scaling distinction between earthquakes and SSEs starts vanishing (see Figure S1 in Supp. Mat.). This observation is in favor of a continuum of modes of slip ranging from slow to rapid, dynamic slip.

We can find some physical intuition about this relative scaling between SSEs and earthquakes in the temporal evolution of rupture length and slip for each event (Fig. 4).

For earthquakes, the average growth of both rupture length and slip are linear with event duration, independent of a/b, hence independent of the actual length of the fault as we non-dimensionalised length scales by L nuc . As a consequence, seismic moment grows quadratically with event duration. In other words, earthquakes propagate as an expanding crack: slip and rupture length are proportional to each other.

For SSEs, however, the temporal evolution of slip and rupture length shows a clear dependence on the fault length. For a given a/b, final rupture length is constant i.e. it is independent of event duration. However, slip grows linearly with duration. If we now increase the fault length (i.e. increase a/b), the accumulated slip decreases (compared to the low a/b case) while the final rupture length increases. These two effects exactly counterbalance each other, such that the final moment scales linearly with duration and is independent of fault length (i.e. for different a/b). This highlights an interesting fact that SSEs are not necessarily self-similar in our calculations.

Finally, we observe that the moment of the nucleation phase scales linearly with its duration. The evolution of slip and rupture length for the nucleation phase is scale independent contrary to SSEs. Slip and final rupture length for nucleation phases evolve, individually, with the square root of the event duration, which might point to a significant difference between these processes.

Stress drop

Interestingly, static stress drops of both slow and rapid slip events in our model are comparable (see Figure S4 in Supp. Mat.). We evaluate this parameter in three different ways following [START_REF] Noda | Comparison of average stress drop measures for ruptures with heterogeneous stress change and implications for earthquake physics[END_REF] (see Supp. Mat. for more details). Regardless of the method, stress drops of SSEs and earthquakes are of similar order of magnitude. Earthquake stress drops are on an average about twice as large as those for SSEs. This is not -10-Confidential manuscript submitted to Geophysical Research Letters completely in agreement with observations where SSEs stress drop is generally 1 or 2 orders of magnitude smaller than for earthquakes [START_REF] Gao | Scaling relationships of source parameters for slow slip events[END_REF]. However, it has also been shown that earthquake stress drops can vary by several orders of magnitude [START_REF] Goebel | Stress-drop heterogeneity within tectonically complex regions: a case study of san gorgonio pass, southern california[END_REF]. Finally, and as expected, the stress drop scales with the moment of individual earthquakes and SSEs. Such observation emphasises the relative importance of slow events in the stress/energy budget of active faults.

Conclusion

We have shown that one simple geometrical complexity (two overlapping faults) can naturally result in a complex seismic cycle (with SSEs, earthquakes, partial ruptures etc.), without appealing to complex friction rheology on the fault. We believe that geometry of fault systems, that have been shown to control the dynamics of ordinary earthquakes [START_REF] Lay | An asperity model of large earthquake sequences[END_REF], are also a primary cause of the source of complexity in the seismic cycle.

In recent years, many models have attempted to explain the nearly ubiquitous presence of slow-slip events in subduction zone. Current models using rate and state friction can only produce slow and fast dynamics in a very narrow range of parameters. Extension of this range required considering additional competing frictional mechanisms. Our work here suggests that complex stress interaction due to geometric complexity of faults could also act as a complementary mechanism to enhance the presence of slow slip in models. This work is an exploratory work on the role of fault geometric complexities in an earthquake cycle. We think that the role of fault geometry in earthquake cycle models has been under-emphasised compared to the role of friction laws in earthquake cycle modelling probably because of the inherent computational limitation of modelling on non-planar geometries. We argue that a unified model that would explain all observations needs to account for geometric segmentation and/or the non-planar nature of active faults as this is a first-order and well documented feature that results in a spatiotemporally inhomogeneous stress accumulation rate [START_REF] Mitsui | Slow slip events controlled by the slab dip and its lateral change along a trench[END_REF][START_REF] Matsuzawa | Comprehensive model of short-and long-term slow slip events in the shikoku region of japan, incorporating a realistic plate configuration[END_REF][START_REF] Li | Spatiotemporal evolution of slow slip events in a nonplanar fault model for northern cascadia subduction zone[END_REF]. As this work shows, the simplest of geometrical complexity can lead to very complex modes of slip on a fault network. Example of a calculation that gives rise to complex slip behaviour on faults. Here L/L nuc = 2, D/L nuc = 0.1 and a/b = 0.9. To avoid any artefact from initial conditions, the first 10 events of the simulation were removed. Left panel shows the maximum slip velocity for fault 1 (blue) and fault 2 (red). Right panel represents the space-time evolution of slip velocity on the faults. The highlighted duration of events corresponds respectively for earthquakes and slow events to the time when the slip velocity exceeds 1mm/s or 1µm/s for the first time to the time when it decelerates below 1mm/s or 1µm/s. Bottom panel gives the geometry used for this example. Events 2,3 and 6 are slow-slip events. Events 1, 4, 5, 7 and 8 are earthquakes.

Event 5 and 7 are small earthquakes that did not rupture the entire fault. Event 1 and 7 clearly show afterslip contrary to events 4 and 8. The table lists the seismological (∆σ M ), spatially averaged (∆σ A ) and slip averaged (∆σ E ) stress drops for the events. et al., 2010;[START_REF] Gao | Scaling relationships of source parameters for slow slip events[END_REF][START_REF] Gomberg | Reconsidering earthquake scaling[END_REF] (top panel) and from our all our calculations (bottom panel). We only used the seismic moment of the dynamic part of an earthquake. The original scaling [START_REF] Ide | A scaling law for slow earthquakes[END_REF] also included data from tremors, very low frequency earthquakes and low frequency earthquake. However because we are not reproducing any of these events, we cut the data to show only slow slip events.

-

  Zigone, D., D. Rivet, M. Radiguet, M. Campillo, C. Voisin, N. Cotte, A. Walpersdorf, N. M. Shapiro, G. Cougoulat, P. Roux, et al. (2012), Triggering of tremors and slow slip event in guerrero, mexico, by the 2010 mw 8.8 maule, chile, earthquake, J. Geophys. Res., 117, B09,304, doi:10.1029/2012JB009160.

Figure 2 .Figure 3 .

 23 Figure 2. Phase diagram showing the evolution of mode of slip along the 2 fault system given the distance between the faults. This figure includes a broader set of simulations in comparison to the paper. Dampeddomain is a domain within which the fault experiences no event at all. SSE & EQ is the domain of coexistence of both slow events and earthquakes. Complex EQ is a domain within which we get only earthquakes but with spatio-temporal complexities. Periodic EQ is a domain within which earthquakes are periodically rupturing the entire fault. And finally, Slip Bursts is a domain within which the entire fault is destabilized at the same time, there is no propagation of the rupture. This corresponds for small faults compared to the nucleation lenghscale and small a/b. This domain is called the no-healing regime[START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF].

Figure 4 .

 4 Figure 4. Final moment, slip and rupture length with time for slow-slip events, earthquakes and nucleation phase of earthquakes.

  Figure1.

L

  Figure2.

  ∈ {0.7, 0.8, 0.85, 0.90, 0.95}. For illustrative purposes we provide a table of dimensional values of L and D in the supplementary section. The smallest faults are 200 m long separated by distance of 21 m. The largest faults are about 20 km long separated by a distance of about 2 km. In fact, it is possible to distinguish between different domains of behavior, that mainly depend on a/b, L/L nuc and the scaled distance between the faults D/L nuc .

5, 1, 2, 3, 4}, and constitutive parameters a/b

  Event ∆σ M (MPa) ∆σ A (MPa) ∆σ E (MPa)

			2.3s
	12.8s
		2h51m
			27.7s
	5.8s
		1h22m	6h07m
	7.3s
	1	2.22	1.99	2.43
	2	0.76	0.94	1.27
	3	0.79	0.83	2.75
	4	2.61	2.57	2.66
	5	2.15	1.94	2.30
	6	0.04	0.08	0.73
	7	0.59	0.75	0.86
	8	2.42	2.21	2.63
	Figure 1.

  Event ∆σ M (MPa) ∆σ A (MPa) ∆σ E (MPa)

			2.3s
	12.8s
		2h51m
			27.7s
	5.8s
		1h22m
			6h07m
	7.3s
	1	2.22	1.99	2.43
	2	0.76	0.94	1.27
	3	0.79	0.83	2.75
	4	2.61	2.57	2.66
	5	2.15	1.94	2.30
	6	0.04	0.08	0.73
	7	0.59	0.75	0.86
	8	2.42	2.21	2.63

Acknowledgments

Numerical computations were performed on the S-CAPAD platform, IPGP, France. P. R. and H. S. B. are grateful to Leslie Greengard and Zydrunas Gimbutas for the FMMLIB2D library. This article benefited from discussions with Robert Viesca and Pierre Dublanchet. P. R. acknowledges the GPX program, funded by the French National Research Agency (ANR), CGG, TOTAL and Schlumberger, for his PhD fellowship. This work received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Geo-4D project, grant agreement 758210).

The software developed for the paper and all the relevant data is available permanently at http://www.geologie.ens.fr/~bhat/romanetGRL2018/.

Confidential manuscript submitted to Geophysical Research Letters

Confidential manuscript submitted to Geophysical Research Letters