Docteur

Rosemary Fayjaloun
email: rose-marie.fayjaloun@univ-grenoble-alpes.fr

Thèse Dirigée Par

CR Christophe Voisin

M Pascal

IPGP, Rapporteur Bernard Physicien

Mme Françoise Courboulex

Pierre-Yves Bard

M Salah

M Christophe Voisin

Mme Cécile Cornou

Laurent Stehly

Pierre Boué

Mathieu Causse

Cecile Cornou

Fabrice Cotton

THÈSE Pour obtenir le grade de

Keywords: pulse, period, magnitude, rupture length, rupture speed

Je remercie les Professeurs Pascal Bernard et Françoise Courboulex qui se sont acquittés de la délicate tâche de rapporteur sur cette thèse. Toute ma gratitude s'adresse aussi à Messieurs Pierre-Yves Bard et Salah Sadek qui m'ont fait l'honneur d'examiner ce manuscrit et de prendre part au jury de soutenance de cette thèse.

Je tiens à remercier mes directeurs de thèse pour l'orientation, pour le temps qu'ils m'ont consacré, pour leurs conseils et pour les pistes pour poursuivre ma recherche ! Merci Cécile d'avoir été à l'initiative de cette aventure. Merci pour ta patience et pour ta confiance. Merci Christophe d'avoir encouragé mes initiatives au travers de la grande liberté d'actions que tu m'as autorisée. Merci Mathieu d'avoir encadré cette thèse et pour ta présence permanente. Merci pour ta relecture attentive de ce manuscrit, pour ta patience et ton soutien qui m'a été précieux afin de mener mon travail à bon port.

Un grand merci aux professeurs Michel Bouchon et Ralph Archuleta, toujours à disposition pour répondre aux questions de manière très cordiale. Mes remerciements s'étendent également à

To Mom and Dad

This thesis is divided into 2 parts. Part 1 focuses on better understanding the seismic rupture process and its relation with the near-fault ground motion. The mechanisms of peak ground motion generating are investigated for homogeneous as well as for heterogeneous ruptures. A quantitative sensitivity analysis of the ground motion to the source kinematic parameters (mainly the distribution of the final fault slip and the rupture speed and the correlation between these parameters) is presented, for sites located in the vicinity of the fault rupture, as well as far from the rupture. A second chapter is dedicated to a major near-fault source effect: the directivity effect. This phenomenon happens when the rupture propagates towards a site of interest, with a rupture speed close to the shear-wave speed (Vs); the waves propagating towards the site add up constructively and generate a large amplitude wave called the pulse. The features of this pulse are of interest for the earthquake engineering community. In this chapter, a simple equation is presented that relates the period of the pulse to the geometric configuration of the rupture and the site of interest, and to the source parameters.

Part 2 is dedicated to better estimate the seismic hazard in Lebanon by simulating the strong ground motion at sites near the main fault (the Yammouneh fault). Lebanon is located in an active tectonic environment where the seismic hazard is considered moderate to high. Historically, destructive earthquakes occurred in the past, the last one dates back to 1202. However, strong motion has never been recorded in Lebanon till now due to the presently infrequent large-magnitude seismicity, and therefore facing an alarming note of potential new ruptures. The Yammouneh fault is a large strike-slip fault crossing Lebanon, making all its regions located within 25 km from the fault. At first, the crustal structure tomography of Lebanon, in terms of Vs, is performed using the ambient noise, in order to characterize the wave propagation from the rupture to the ground surface. To our knowledge, this is the first study of the 3D Vs tomography in Lebanon. Afterwards, a hybrid approach is presented to simulate broadband near-fault ground motion (0.1 -01 Hz). At lowfrequencies (≤ 1 Hz), potential ruptures of Mw 7 are simulated (capturing the source effect in the near-fault region defined in the previous chapters), and the generated slip rate functions are convolved with the Green's functions computed for the propagation medium defined by the Vs tomography. The ground-motion is complemented by a high-frequency content (up to 10 Hz), using a stochastic model calibrated by nearfault recordings and accounting for the presence of the directivity pulse. The computed peak ground acceleration is compared to the design acceleration in Lebanon.

Key words: Near-fault peak ground motion, Source numerical (kinematic) simulation, Peak ground acceleration, Sensitivity analysis, Source parameters, Directivity effects, Pulse period, Ambient noise, 3D crustal Tomography, Lebanon, Seismic hazard, Near-fault ground motion hybrid model.

Résumé

Les données accumulées sur les mouvements du sol apportent des connaissances très importantes sur les processus de rupture des séismes, les caractéristiques du milieu de propagation, les effets d'amplification du site sur le mouvement du sol et la relation entre le mouvement du sol et les dommages des structures. Cependant, les séismes de faible et moyenne amplitude étant plus fréquents que les grands événements sismiques, les bases de données de mouvements de sol utilisées dans le développement de modèles de prédiction du mouvement du sol ne contiennent pas beaucoup de données de forts séismes. Le point le plus critique concerne les stations proches de la rupture de la faille, dites en champ-proche, pour lesquelles les bases de données restent mal échantillonnées. C'est pourtant là où les secousses sismiques sont les plus fortes. Les pays à sismicité modérée ou élevée pour lesquels des failles majeures peuvent se briser à proximité de ses grandes villes, sont donc confrontés à un risque sismique majeur, mais le manque d'enregistrements du mouvement ne permet pas une bonne prédiction des mouvements fort du sol. Il est donc nécessaire de simuler le mouvement fort. L'un des défis actuels en sismologie est la mise au point de méthodes fiables pour simuler les mouvements du sol à proximité des failles, en tenant notamment compte du manque de connaissance sur les caractéristiques d'une rupture potentielle.

Cette thèse est divisée en 2 parties. La partie 1 se concentre sur une meilleure compréhension de la rupture sismique et de son rapport avec le mouvement du sol proche de la faille. Les mécanismes de génération des valeurs de pics du mouvement du sol sont étudiés pour des ruptures homogènes et hétérogènes. Une analyse quantitative de sensibilité du mouvement du sol aux paramètres cinématiques de la rupture (principalement la distribution du glissement et de la vitesse de rupture, et la corrélation entre ces paramètres) est présentée, pour des sites au voisinage de la rupture ainsi qu'en champ lointain. Un second chapitre est consacré à un effet de source majeur en champ proche: l'effet de directivité. Ce phénomène se produit lorsque la rupture se propage vers un site, avec une vitesse de rupture proche de la vitesse de l'onde de cisaillement Vs; les ondes se propageant vers le site interfèrent de manière constructive et génèrent une onde de grande amplitude appelée « pulse ». Les caractéristiques de ce pulse, notamment sa durée, représentent des paramètres d'intérêt pour le génie parasismique. Dans ce chapitre, une équation simple est présentée pour relier la durée du « pulse » à la configuration géométrique de la rupture et du site d'intérêt et aux paramètres de la source.

La partie 2 est consacrée à une meilleure estimation de l'aléa sismique au Liban en simulant le mouvement fort pour des sites proches de la faille principale : la faille de Yammouneh. Le Liban est situé dans un environnement tectonique actif où le risque sismique est considéré comme modéré à élevé. Historiquement, des tremblements de terre destructifs se sont produits dans le passé, le dernier remontant à 1202. Cependant, en raison de la sismicité de grande ampleur actuellement peu fréquente, aucun mouvement fort n'a jamais été enregistré au Liban à ce jour. La période de retour de ces séismes historiques ayant été atteinte, le Liban est soumis à de nouvelles grandes ruptures potentielles. La faille de Yammouneh est une grande faille en décrochement traversant le Liban du Nord au Sud, situant toutes les villes et infrastructures à moins de 25 km de la faille. Dans un premier temps, une tomographie de la structure de la croûte du Liban, en termes de vitesse des ondes de cisaillement Vs, est réalisée en utilisant le bruit ambiant, pour caractériser la propagation des ondes sismiques de la source à la surface. À notre connaissance, il s'agit de la première étude de la tomographie Vs 3D au Liban. Par la suite, une approche hybride est utilisée pour simuler le mouvement du sol en champ proche sur une large bande de fréquences (0.1 -10 Hz). Aux basses fréquences (≤ 1 Hz), des ruptures potentielles de Mw7 sont simulées (prenant en compte l'effet de la source en champ proche définie dans les chapitres précédents), et les fonctions sources obtenues sont convoluées aux fonctions de Green calculées pour le modèle de propagation des ondes issu de la tomographie Vs afin d'estimer le mouvement du sol à proximité de la faille. Le mouvement du sol est complété par un contenu haute fréquence (jusqu'à 10 Hz), en utilisant un modèle stochastique calibré par des enregistrements en champ proche, et en tenant compte de la phase impulsive due à la directivité de la rupture (« pulse »). L'accélération maximale du sol calculée est comparée à l'accélération réglementaire utilisée pour la conception sismique au Liban.

Mots clés: Mouvement du sol en champ proche, Simulation numérique (cinématique) de la source, accélération maximale du sol, Étude de sensibilité, Paramètres de la source, Effet de directivité, Période du pulse, Bruit ambiant, tomographie 3D de la croûte, Liban, Alea sismique, Modèle hybride.

Table of Contents

Introduction

Throughout the ages, earthquakes have been one of the most destructive natural hazards. Between 2000 and 2015, more than 800,000 people were killed worldwide because of earthquakes, according to the United States Geological Survey (USGS). 2004 was the deadliest year among the last two decades with more than 230,000 people killed by the M9.1 earthquake in Sumatra (Indonesia) and the triggered tsunami. Although these risks cannot be prevented, their effects can be lessened. In order to mitigate risk, one must first view the problem in its entirety portrayed by the relation [START_REF] Ambraseys | Earthquakes in the Mediterranean and Middle East: A Multidisciplinary Study of Seismicity up to 1900[END_REF]):

[Earthquake Risk] = [Structural Vulnerability] * [Earthquake Hazard]
There is a clear distinction between structural vulnerability and earthquake hazard. Structural vulnerability is the degree of structural damage or loss resulting from an earthquake of a given magnitude and is the subject matter of earthquake engineering. It is determined by the physical characteristics of structures. Measures to reduce vulnerability can be thought of either as long-term, e.g. earthquakeresistant design and construction, appropriate physical planning of settlements, or as short-term action in response to the post-earthquake hazard. Earthquake hazard is the probability of a damaging earthquake happening within a specific period of time and given area. It is still beyond human control. Assessment of earthquake hazard is a subject matter of Earth sciences. At the present level of technology, earthquakes cannot be predicted. Since we cannot know what will happen in the future, to estimate earthquake hazards we have to find out what happened in the past and extrapolate from there. Strong ground motion prediction is one of key factors for mitigating disasters for future earthquakes. Ground motion observations are the result of a long history of instrument development and use, allowing strong ground motion prediction for mitigating disasters for future earthquakes. The world's first accelerogram was recorded in California in 1933, and the first to be recorded in Europe was 30 years later in Macedonia in 1967 [START_REF] Bommer | Current Empirical Ground-Motion Prediction Equations for Europe and Their Application to Eurocode 8[END_REF]. Since the late 1960's, a rapid growth in the deployment of special instrumentation to record strong ground motions throughout the seismic regions of the world has resulted in the accumulation of a large data bank of useful records from key sites like California, USA, Italy, Greece, Japan, New Zealand and Turkey. 5.5. In Japan, the National Research Institute for Earth Science and Disaster Prevention (NIED) operates the Strong-Motion Sismograph Networks (K-NET, KiKnet) with 660 strong motion stations. Each station records triaxial accelerations both at the surface and at sufficient depth in rock to understand the physics of earthquake fault rupture and to directly observe linear and nonlinear seismic wave propagation in the shallow crust. These borehole-surface data have provided fundamental new constraints on peak ground motions [START_REF] O'connell | Earthquake Ground Motion Estimation[END_REF], direct observation of nonlinear wave propagation, and new constraints on ground motion variability [START_REF] Rodriguez-Marek | Analysis of Single-Station Standard Deviation Using the KiK-Net Data[END_REF]. Global and national databases were also developed: COSMOS Consortium of Organizations for Strong Motion Observation Systems (San Francisco), Center for Engineering Strong-Motion Data (California), Strong Motion Database of Turkey, ITalian ACcelerometric Archive, RESORCE and European Strong Ground Motion Database (Europe), Euroseistest database (Greece), Unified HEllenic Accelerogram Database (Greece), Swiss National Strong Motion Network, Institute of Engineering Seismology and Earthquake Engineering (Greece), Guerrero Accelerograph Network (Mexico), GeoNet (New Zealand)… Accumulated data of strong ground motions have been providing us very important knowledge about rupture processes of earthquakes, propagation path, site amplification effects, relation between ground motion and damage, and so on.

Researchers developed Ground Motion Prediction Equations (GMPEs) that relate a ground motion parameter (e.g., peak ground velocity and acceleration, PGV and PGA respectively; acceleration response spectra) to a set of explanatory variables describing the earthquake source, wave propagation path and local site conditions (e.g. [START_REF] Douglas | Earthquake Ground Motion Estimation Using Strong-Motion Records: A Review of Equations for the Estimation of Peak Ground Acceleration and Response Spectral Ordinates[END_REF]. In their simplest form, these empirical GMPEs predict peak ground motions based on a limited parametric description of earthquake and site characteristics. These variables include the earthquake magnitude and faulting mechanism, the reduction (attenuation) of ground motion amplitudes with increasing distance from the fault (geometric spreading), and the local site characteristics using either site classification schemes or a range of quantitative measures of shallow to deeper velocity averages or thresholds. Peak ground motion amplitudes generally increase with increasing magnitude up to a threshold magnitude range where peak accelerations saturate, i.e., only slightly increase or stay nearly constant above the threshold magnitude range [START_REF] Campbell | Near-Source Attenuation of Peak Horizontal Acceleration[END_REF]. Similarly, observed peak ground motion amplitudes decrease with increasing distance from the earthquake fault, but saturate at close distances to faults such that the decrease in amplitudes with increasing distance is small within several km of faults. Recent GMPEs are also parameterized with rupture directivity effects, hanging wall effects, non-linear magnitude scaling and seismic wave propagation distinction (e.g. high and low Q zone, Moho reflection and deep basin effects) [START_REF] Abrahamson | Summary of the ASK14 Ground Motion Relation for Active Crustal Regions[END_REF][START_REF] Boore | NGA-West2 Equations for Predicting PGA, PGV, and 5% Damped PSA for Shallow Crustal Earthquakes[END_REF][START_REF] Campbell | NGA-West2 Ground Motion Model for the Average Horizontal Components of PGA, PGV, and 5% Damped Linear Acceleration Response Spectra[END_REF][START_REF] Idriss | An NGA-West2 Empirical Model for Estimating the Horizontal Spectral Values Generated by Shallow Crustal Earthquakes[END_REF])… In the past five decades many hundreds of GMPEs for the prediction of PGA and linear elastic response spectral ordinates have been published [START_REF] Douglas | Ground-Motion Prediction Equations 1964-2010[END_REF]. The range of application for an empirical ground motion prediction model is constrained by the range of the data coverage used in the analysis [START_REF] Arroyo | On the Forecasting of Ground-Motion Parameters for Probabilistic Seismic Hazard Analysis[END_REF].

However, since the occurrence of small and moderate magnitude earthquakes is more frequent than the occurrence of large seismic events, most of ground motion databases used in the development of ground motion prediction models are primarily comprised of accelerograms produced by small and moderate earthquakes. Hence, as magnitude increases, the sets of ground motions become sparse. For instance, in the database used in the Next Generation of Ground Motion Attenuation Models project (NGA-West2), there are only nine earthquakes with M>7 and two of them yielded only two recordings per earthquake (e.g. [START_REF] Chiou | Update of the Chiou and Youngs NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF].

In the beginning, due to lack of enough information, seismologists were not able to predict the location and the rupture zone for a large earthquake, and it was difficult to find the strategic placement for the seismographs (at first very costly and difficult to maintain) in order to capture the best details of the earthquake near the zone of the rupture. As a result, almost all recorded ground motions were of earthquakes far from the rupture area [START_REF] Cox | Characterization of Large Velocity Pulses for Laboratory Testing[END_REF], where the ground shaking has been attenuating as distance from the fault increases. Awareness of the importance of the near-fault effects date back as far as the 1971, after the San Fernando-California earthquake that was one of the first instances of building collapse that was associated with near-fault source effect. Serious concern was raised following the 1994 Northridge-California earthquake and the Kobe earthquake in Japan (1995). Nonetheless, for large magnitude short distance ranges, databases are poorly sampled; the number of strong ground recordings in the near fault is not large enough to develop GMPEs based on recorded data. Moreover, GMPEs take into account the source rupture by the magnitude parameter only. However, this source parameter is not sufficient to describe the source process: the empirical upper bound on peak accelerations are showing only a weak magnitude dependence [START_REF] Anderson | An Overview of the Largest Amplitudes in Recorded Ground Motions[END_REF].

GMPEs make efforts to take into account the directivity effect, one of the major source effects, that is significant in the near-fault. This phenomenon usually happens when the rupture propagates towards a site of interest; the energy reaches the site adds up constructively, and generates a large amplitude wave called the pulse. To represent the directivity pulse-like features, five directivity models were developed as corrections to be added to the current median NGA-West2 GMPEs. Those directivity models modify the response spectra predicted by existing GMPEs: either by amplifying response spectra at all periods (called broadband models, e.g., [START_REF] Somerville | Modification of Empirical Strong Ground Motion Attenuation Relations to Include the Amplitude and Duration Effects of Rupture Directivity[END_REF][START_REF] Spudich | Directivity in NGA Earthquake Ground Motions: Analysis Using Isochrone Theory[END_REF][START_REF] Bayless | Bayless-Somerville Directivity Model[END_REF], or by only amplifying response spectra within a narrow range of periods that depends on the earthquake magnitude (called narrowband models, e.g., [START_REF] Shahi | An Empirically Calibrated Framework for Including the Effects of Near-Fault Directivity in Probabilistic Seismic Hazard Analysis[END_REF]Shahi 2013;[START_REF] Spudich | The Spudich and Chiou NGA-West2 Directivity Model[END_REF]). However, these models need further refinement [START_REF] Spudich | Comparison of NGA-West2 Directivity Models[END_REF][START_REF] Spudich | Final Report of the NGA-West2 Directivity Working Group[END_REF]. Despite continuous expansion of the database of recorded earthquake ground motions, recorded near-fault ground motions exhibiting forward directivity pulses remain scarce to perform a regression analysis for the features of the pulses. Hence, numerical models are a relevant approach for properly modeling the near-fault source effect like the directivity effects.

When databases of strong motion recordings are too sparse to develop empirical attenuation relations from recorded data, strong motion simulations become then very useful. Strong motion simulation methods have the advantage of allowing the incorporation of information about earthquake source, seismic wave propagation, and local site characteristics that are specific to the region and to the site in question. These characteristics may include rupture directivity effects, hanging wall/foot wall effects, Moho bounce effects, and site effects (for instance 3D basin or topographic effects). Besides, on the contrary to GMPEs, numerical simulations provide seismic ground motion time series, which are useful for nonlinear dynamic analyses of structures needed to design earthquake-resistant buildings and critical structures such as bridges, lifelines, and electric power plants.

The strong ground motion, especially near-fault, is highly sensitive to the fault rupture process (directivity effects, supershear ruptures, rupture nucleation hypocenter, rupture roughness…). A major issue in simulating ground motion for future potential earthquake in near-and far-fault remains the definition of these rupture parameters, which are poorly known and then difficult to estimate a priori. Therefore, there is still a real need of in depth investigations of the effects of rupture parameters on ground motion in order to extract rupture parameters that mostly impact ground motion. In addition, ground motion simulation need calibration by comparing with the predicted motions predictions by GMPEs (and/or the observed ground motions if available). It is important to note that both schemes of ground motion simulation and GMPE have no controversial role and complement each other. The goal of this PhD work is to simulate broadband ground motion for moderate to large event in the vicinity of a fault, incorporating the complexity of the rupture process.

PhD thesis flowline

This thesis is divided into 2 parts. The first part of my PhD thesis focuses on understanding the impact of various source parameters and the directivity effects of the source on the ground motion near the fault. The second part is dedicated to extract the shear-wave velocity model at crustal scale in Lebanon and, hereafter, to simulate the ground motion and its variability that would occur at various sites in Lebanon for a M7 earthquake on Yammouneh fault that crosses the country.

Part 1 begins with an overview (chapter 1) about the different approaches for predicting the seismic ground motion. Moreover, it sheds the light on the difficulties encountered when it comes to defining the values of seismic source parameters. An important challenge of the current seismology is the development of reliable methods of simulating ground motion near the fault, taking into account the lack of knowledge of the rupture characteristics. Because source parameters are often not readily available or fully understood, the resultant uncertainties of source characterization can be the dominant contributions to uncertainty in ground motion prediction. Nowadays, researchers are capable of simulating very large set of ground motions in complex 3D structures ( e.g. [START_REF] Moschetti | 3D Ground-Motion Simulations of Mw 7 Earthquakes on the Salt Lake City Segment of the Wasatch Fault Zone: Variability of Long-Period (T≥1 s) Ground Motions and Sensitivity to Kinematic Rupture Parameters[END_REF]). However, it remains very hard to well-describe the many source parameters in order to evaluate the seismic hazard for a given fault. Given the complexity and our lack of knowledge of the rupture process in terms of fault slip, speed of the rupture front and stress drop, among other source parameters, one cannot escape to statistical approaches to describe the source. Chapter 2 (Manuscript Submitted for Publication) has thus the main goal to better understand what source parameter controls the most the ground motion in the near fault. It presents a quantitative sensitivity analysis of the surface ground motion (up to 5 Hz) to several source kinematic parameters (in particular level of source heterogeneity and correlation between rupture parameters), in the vicinity of the fault rupture, as well as far from the rupture. In this chapter, the mechanisms of generating the peak ground motion are illustrated for homogeneous as well as for heterogeneous ruptures.

Near-fault ground motions possess distinct characteristics that can have strong influence on structural response due to forward directivity (e.g. [START_REF] Baker | Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis[END_REF][START_REF] Spudich | Comparison of NGA-West2 Directivity Models[END_REF]): when the rupture propagates from the hypocenter toward a site near the fault, the site is said to be located in the forward directivity region; when the rupture front propagates toward the site and at a velocity almost equal to the shear-wave velocity of the ground, all the seismic energy radiated from the fault rupture arrives at the site in a single, short-duration pulse. Hence, the site may experience a large-amplitude, short duration pulse at the beginning of the velocity time series [START_REF] Somerville | Development of an Improved Representation of near Fault Ground Motions[END_REF]. The forward directivity effect is reflected by a peak in the response spectrum near the period of the directivity pulse [START_REF] Somerville | Conditions That Give Rise to Unusually Large Long Period Ground Motions[END_REF]). [START_REF] Baker | Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis[END_REF][START_REF] Mavroeidis | A Mathematical Representation of Near-Fault Ground Motions[END_REF][START_REF] Shahi | An Empirically Calibrated Framework for Including the Effects of Near-Fault Directivity in Probabilistic Seismic Hazard Analysis[END_REF][START_REF] Bray | Characterization of Forward-Directivity Ground Motions in the near-Fault Region[END_REF][START_REF] Alavi | Consideration of Near-Fault Ground Motion Effects in Seismic Design[END_REF] proposed empirical equations to predict the period of the directivity pulses as function of magnitude. In chapter 3 (BSSA Publication), we present a simple model to predict the pulse period. The parameters for the equations are related to the source rupture process, and depend on the location of the observation points with respect to the rupture, on the rupture speed, and on the rise time. This model, though simple, fairly well explains the spatial variability of the pulse periods observed.

Part 2 focuses on simulating ground motion for a M7 earthquake in Lebanon within the frequency range from 0.1 to 10 Hz. Although located in low to moderate seismicity region, Lebanon faced several destructive earthquakes in the past. Indeed Lebanon is crossed by the Dead Sea Fault that splits into four main fault branches in Lebanon, each branch being capable to produce M7+ earthquake as observed in the past [START_REF] Nemer | Coseismic Growth of Sedimentary Basins along the Yammouneh Strike-Slip Fault (Lebanon)[END_REF][START_REF] Daëron | 12,000-Year-Long Record of 10 to 13 Paleoearthquakes on the Yammoûneh Fault, Levant Fault System, Lebanon[END_REF][START_REF] Elnashai | Earthquake Hazard in Lebanon[END_REF]). The faulting system in Lebanon, which is a small country, makes all its regions located no more than 25 km away from a large fault. Part 2 begins with an overview about the Dead Sea Fault in the Middle East region, and more particularly in Lebanon, its historical and recent seismicity with an overview of its geology (chapter 4). However, the Earth crustal velocity structure in Lebanon is largely unknown. Chapter 5 (Manuscript in Preparation for Submission) focuses on the shear-wave tomography by means of seismic ambient interferometry. Finally, chapter 6 wraps up the main findings from each of the preceding sections on rupture parameters and crustal structures in Lebanon in order to simulate the ground motion in Lebanon for critical, however, realistic scenarios for a M 7 rupture along the northern part of the Yammouneh Fault, a strike-slip fault crossing Lebanon. Source parameters are chosen according to chapter 2, where subshear and supershear ruptures (source rupturing at a speed smaller or larger than the shear wave speed) are considered, and the wave propagation medium is based on chapter 5. A broadband hybrid model to simulate strong ground motion on a broad frequency range (~0.1 -10 Hz) is finally presented, that combines pseudo-dynamic source rupture models reflecting the physics of the rupture and the directivity effects, and stochastic modeling approach calibrated to worldwide recordings of large earthquakes in the near-fault area to cover a broadband frequency range while respecting the characteristics of the low-frequency ground motion. The computed peak ground acceleration is compared to the design acceleration in Lebanon. An earthquake may last only a few seconds, but the processes that cause earthquakes have operated within the Earth for millions of years. Over the course of geological time, earthquakes and other natural events have helped to shape the surface of our planet. For centuries, people wondered what caused the Earth to shake. In the 1960s, scientists settled on the theory of plate tectonics [START_REF] Oreskes | Plate Tectonics: An Insider's History of the Modern Theory of the Earth[END_REF][START_REF] Ohnaka | The Physics of Rock Failure and Earthquakes[END_REF]): Although our Earth feels solid as we walk along its surface, it is really only partly so. The Earth is divided into three main layers: there is a hard outer surface (the crust), a softer middle layer (the mantle), and a central core. The crust and the upper portion of the mantle are referred to together as the lithosphere, with an average depth of 100 km. The outermost layer of the Earth is broken into irregular pieces, called the Earth's lithospheric plates. These pieces are not static but in very slow constant motion. The convection currents are the force that drives the plates. It is a system of heat exchange that forms in the Earth's mantle: beneath the lithosphere, the mantle is semi-molten to a depth of about 260 km. Its plastic-like material rises in response to heat and sinks when the temperature drops. This convective movement acts as a drag on the underside of the lithospheric plates, causing them to move. Plates move in three different ways: colliding with each other (convergent movement), spreading apart (divergent movement), or sliding past one another (strike slip or transform movement). On the boundaries of these plates are faults, which stick together while the rest of the plate keeps moving. As a result of plate motions, the rocks are either squeezed (they are under compressional stress), or being pulled apart (they are under tensional stress). The rocks will behave elastically; that is, they absorb the shear stress by changing their shape, and the change in shape is called strain. When this occurs, the energy that would normally cause the plates to move past one another is stored up, until eventually, the fault strength reaches its limit, the fault becomes unstable and begins to slip: the fracture nucleates at a point -the hypocenter -and propagate with a rupture velocity 𝑉 𝑟 . The rock mass on either side moves abruptly, and the strain energy accumulated between the pieces is suddenly released. This sudden release is what we call elastic rebound [START_REF] Coburn | Earthquake Protection[END_REF]. One measure of the size of the earthquake is the moment magnitude Mw (or M used interchangeably), based on estimated rupture length L and width W, average slip 𝐷 ̅ , and the average rigidity 𝐺 over the rupture area.

Before slippage, the initial value of shear stress is 𝜏 0 . After slip motion has stopped, the shear stress has reached its final value 𝜏 1 . During slippage, work is done against the frictional stress 𝜏 𝑓 . With these notions, the static stress drop is thus ∆𝜏 = 𝜏 0 -𝜏 1 , and the dynamic stress drop is ∆𝜏 𝐷 = 𝜏 0 -𝜏 𝑓 . In this thesis, if the parameter 'stress drop ∆𝜏' is mentioned, it is static stress drop that is being referred to. During this process, the potential energy (strain and gravitational energy) of the system is lowered by ∆𝑊 = Ε 𝑠 + Ε 𝑓 + Ε 𝑔 (Figure 1-1), where Ε 𝑠 is the energy radiated in form of seismic waves, Ε 𝑓 is the frictional energy loss in form of heat and Ε 𝑔 is the energy expended to create new fault surface, also called fracture energy [START_REF] Rivera | Representations of the Radiated Energy in Earthquakes[END_REF]. Earthquakes have caused, and can cause in the future, enormous loss of life, injuries, destruction of properties, and economic and social disorder. We cannot accurately predict when and where an earthquake will occur. Although scientists do create sophisticated models of earthquakes and study the history of quakes along fault lines, no one has enough of an understanding about the state of stress on active faults at depths, and how rupture initiates and stops in a complex medium characterized by a variable stress field acting on complex fault geometries. Likewise, scientists do not have a full understanding of the conditions -the rock materials, minerals, fluids, temperatures and pressures -at the depths where quakes start and grow… to be able to predict them. Our observations of earthquakes are always at distance, viewed indirectly through the lens of seismic waves, surface faulting and ground deformation. To predict earthquakes, we would need to have a good understanding of how they occur, what happens just before and during the start of an earthquake [START_REF] Geller | Earthquakes Cannot Be Predicted[END_REF]). Although identifying the exact time and size of an earthquake is currently impossible, scientists can estimate the probability of an earthquake occurring in a region or on a fault over a span of decades. To do so, we need information about how fast the fault is sliding over the long term -typically a few millimeters to centimeters of slip per year -and how big the earthquakes are likely to be. We calculate how much slip is used up in each earthquake, and thus how often earthquakes must occur (the return period), on average, to keep up with the long-term slip rate. Knowing the date of the last earthquake helps improve forecasting. Note however, that any time could mean tomorrow or 100 years from now.

Although we cannot predict the timing of an earthquake, we can predict the damaging seismic waves generated by a potential earthquake. Determining the earthquake intensity does help developers make good decisions about where to build and what type of forces those buildings should be constructed to withstand. If the buildings are strong, we will be safe no matter when the ground happens to shake, and we can use that knowledge to make our communities and ourselves resilient. The most complete description of seismic ground motion is given by the recorded time series (Figure 1-2) on the surface from which we can extract parameters characterizing the amplitude, duration and frequency content. The records provide the displacement, the velocity or the acceleration time history during the earthquake. The maximum value of the ground motion is an important seismic quantity (e.g., Peak ground acceleration 𝑃𝐺𝐴 is the maximum absolute values of the acceleration waveforms). As well, the frequency content provides valuable information on the earthquake source as well as the Earth structure. For earthquake engineers, the response of the structures towards earthquake shaking is what matters. Buildings and bridges are very sensitive to the frequency at which they are loaded [START_REF] Kramer | Geotechnical Earthquake Engineering[END_REF]. A structure suffers the greatest damage from the ground motion shaking at a frequency close or equal to its natural resonance frequency. And therefore, each structure performance depends on the shape (configuration) of the building, the construction material, the height of the building… As a rule of thumb, the natural frequency of a structure 𝑓 0 is related to the number of stories 𝑁 by statistical relation like 𝑓 0 = [ Buildings are designed to resist certain level of force applied to it, which is proportional to the response acceleration of each structure 𝑆𝑎 (𝑇 0 ) . Seismologists study the acceleration ground motion since the acceleration is related to the loads on the structures. We normally associate acceleration (the second derivative of the displacement) with high frequencies, because, in the frequency domain, the differentiation introduces the high-frequency amplification: spectral components are amplified in direct proportion to their frequencies as shown in equation (1-1):

𝑭𝑻 ( 𝒅 𝜶 𝒇(𝒕) 𝒅𝒕 𝜶 ) = (𝒊𝝎) 𝜶 𝑭(𝝎) (1-1)
Accordingly, the acceleration maximal value (PGA) is mainly controlled by the high frequencies.

And since the damage of a building depends on the earthquake ground acceleration (frequency and amplitude content), and since the acceleration maximal value (PGA) is mainly controlled by the high frequencies, the principal objective of earthquake engineering and engineering seismology is thus to provide quantitative and reliable estimates of expected levels of seismic ground-motion and response of the structures, especially at high frequencies.

Fundamental equations for earthquake ground motion

Earthquakes generate seismic waves that travel from the source to the surface and cause surface ground motions over a wide range of frequencies. Indeed, seismic waves are generated as part of the strain energy released from the rupturing of a fault (source effect). The seismic waves then propagate through the Earth's geological structure (wave propagation or path effect) and approach the surface of the Earth, where they undergo further modifications while propagating through shallow soils (site effect). The ground motion recorded at the surface therefore is the end product of the interaction between source, path and site effects.

The wave equation and the Green's functions

A large part of our understanding of the physics of earthquakes and wave propagation is recapitulated in the elasto-dynamic seismic wave equation, which is basically, the combination of two fundamental theorems: Newton's second law that connects forces in a continuous medium to observable displacements and Hooke's law that relates stress and strain in a linear elastic medium. Assuming a homogeneous (i.e., the elastic properties 𝜌 0 , 𝜆 and 𝜇 are constants) and isotropic medium (i.e., the elastic properties are equal in all directions), the equation of motion is presented in equation (1-2) (Aki and Richards 2002):

𝝆 𝟎 𝒖 ̈= (𝝀 + 𝑮) 𝛁(𝛁. 𝒖) + 𝑮𝛁. 𝛁𝒖 + 𝒇 (1-2)
in terms of the motion variables 𝒖 and 𝒖 ̈ and the seismic source term 𝒇, where 𝜆 and 𝜇 are the Lamé parameters. When the fault ruptures, the energy is released and it radiates outward through the ground in the form of body waves. Body waves that travel through the Earth are either P-(for Primary) waves or S-(for Secondary) waves, where the P-wave velocity, 𝑉 𝑃 , and the S-wave velocity, 𝑉 𝑠 , are respectively given by equations (1-3) and (1-4):

𝑽 𝑷 𝟐 = 𝝀 + 𝟐𝑮 𝝆 𝟎 (1-3) 𝑉 𝑠 2 = 𝐺 𝜌 0 (1-4)
We can use the latter equations to rewrite the elastic wave equation directly in terms of the P and S velocities, presented in equation (1-5):

𝒖 ̈= 𝑽 𝑷 𝟐 𝛁𝛁. 𝒖 -𝑽 𝒔 𝟐 𝛁 × 𝛁 × 𝒖 + 𝒇 (1-5) P-waves travel faster than S-waves. The two types together are called body waves because they travel through the body of the Earth, on the contrary to surface waves, which travel along the surface of the Earth. Such waves do not exist in homogeneous media, but in the case of layered media or at the free surface of the media. When 𝒇 is an impulse double couple of unit forces, the displacement 𝒖 is called Green's function, which represents the response of the medium.

The Representation theorem and the far field approximation

The source term should describe the rupture process over a finite extent fault plane.

It is then necessary to make use of the representation theorem which relates slip on the fault to ground motion at an arbitrary point in the field. It allows computing the ground displacements resulting from an earthquake if the slip on the fault plane is known. It consists of integrating the Green's functions weighted by the displacement discontinuity over the fault surface, as shown in equation (1-6):

𝒖 𝒏 (𝒙, 𝒕) = ∫ 𝒅𝝉 ∬ [𝒖 𝒋 (𝝃, 𝝉)] 𝒄 𝒋𝒌𝒑𝒒 𝑮 𝒊𝒑,𝒒 (𝒙, 𝒕; 𝝃, 𝝉) 𝒗 𝒌 𝒅𝚺 𝚺 +∞ -∞ (𝝃) (1-6)
where 𝑐 𝑗𝑘𝑝𝑞 are the elastic moduli, 𝐺 𝑖𝑝,𝑞 (𝒙, 𝑡; 𝝃, 𝜏) is the derivative of the Green function with respect to 𝜉 𝑞 (refer to (Aki and Richards 2002) for a complete definition of parameters). At distances far from the fault, that is, at distances larger than a few rupture lengths, the fault rupture can be approximated by an equivalent point source with seismic moment 𝑀(𝑡) = 𝐺 ∬ 𝑢(𝝃) 𝑑Σ Σ (𝝃). In addition, the near field terms of the body waves become negligible. Equation (1-6) can then be approximated by equation (1-7). According to this equation, the displacement recorded at the stations looks bell-shaped and identical with the source time function prescribing the slip velocity evolution during rupture propagation along the fault (Figure 1-3). This leads to the concept of the source seen as an effective point source from far away and the seismic waves generated by this point source may be described in terms of a double couple moment tensor (Aki and Richards 2002). Hence the farfield approximation states that in a homogeneous, isotropic, unbounded medium, at distance larger than the source dimension and several wavelengths of the considered signal, the displacement is proportional to the moment rate function 𝑀 ̇(𝑡), i.e. the time derivative of the seismic moment 𝑀 0 . Considering a shear wave, the far-field displacement is then given by equation (1-7):

𝒖 𝑭𝑭 (𝑿, 𝒇) = 𝟏 𝟒𝝅𝝆 𝟎 𝑽 𝒔 𝟑 𝟏 𝐗 𝐑𝐏 𝑴 ̇(𝒇) 𝒆 - 𝟐𝝅𝒊𝒇𝑿 𝑽 𝒔 𝒆 - 𝝅𝒇𝑿 𝑽 𝒔 𝑸 𝒔 (1-7)
Where, 𝑋 is the distance from the rupture, 𝜌 0 is the rock density, 𝑉 𝑠 is the shear wave speed, 𝑅𝑃 is the radiation pattern of the shear wave, 𝑄 𝑠 is the attenuation factor, and 𝑓 is the frequency. The attenuation of ground motion is simply modeled by a 1/𝑋 geometrical attenuation and by the term 𝑒𝑥𝑝 (-𝜋𝑓𝑋 𝑉 𝑠 𝑄 𝑠

) to represent the anelastic attenuation. This simple relation between the waveforms and the fault slip function will be used in this thesis to analyze the sensitivity of the far-field ground motion to some source properties. However, this relation does not include the details of the rupture as a finite-extent process, and therefore a complete determination of the slip function requires observations near the seismic source. In addition, during propagation, the waves are affected by scattering, spreading, focusing, multipath interference and other complex path effects. One way of minimizing the effects of the path complexity is to make observations at short distance from the seismic source, therefore near-fault data (at small distance from the fault) are essential for more complete study of the source mechanism (Aki and Richards 2002).

1.1.3 Earthquake ground motion modeling 1.1.3.1 Ground motion prediction equations (GMPEs) and ground motion variability

Strong motion databanks offer the possibility to derive empirical or semi-empirical ground motion prediction equations (𝐺𝑀𝑃𝐸) to estimate ground-shaking levels for future earthquake. The GMPEs relate predictor variables 𝑌 such as PGA to variables such as earthquake magnitude, source-to-site distance, faulting style, and site class [START_REF] Abrahamson | Summary of the Abrahamson & Silva NGA Ground-Motion Relations[END_REF][START_REF] Boore | NGA-West2 Equations for Predicting PGA, PGV, and 5% Damped PSA for Shallow Crustal Earthquakes[END_REF]Bindi, Spallarossa, and Pacor 2017b), describing the source scaling, the attenuation with distance and the site amplifications. The global world databank (NGA-West2) of strong motion data contains thousands of records thus providing estimates of the ground shaking as a function of magnitude and distance from the rupture. On the other side, the analysis of such databases outline the very high degree of variability of the ground motion parameters even at a given distance for a given earthquake (K Campbell and Bozorgnia 2014). This drives seismologists into better understanding, quantifying and modeling the observed variability of earthquake ground motion in order to correctly estimate expected levels of seismic ground-motion. Ground motions computed using GMPEs are then given in terms of the median of the natural logarithm of 𝑌 (𝜇 ln(𝑌) ) and its standard deviation (𝜎 ln(𝑌) ), typically referred to as ground motion variability. Ground motion variability associated with ground motion prediction results from imperfect modeling, that is, uncertainties on the model parameters (epistemic uncertainty) or partial knowledge about the physical processes driving the ground motion (aleatory variability). The term 𝜎 ln(𝑌) can then be further subdivided into an intraevent (also called within-event W) component 𝜑 ln(𝑌) (i.e., the average variability in site conditions and path effect over all event) and an interevent (also called between-event B ) component 𝜏 ln(𝑌) (i.e., the variability due to the natural source randomness) (L. Al Atik et al. 2010), illustrated in Figure 1-4. The between event variability is the variability in ground motion one would expect at a given station that recorded many different earthquakes of the same magnitude and corrected for path effects (David M Boore and Atkinson 2008; K.W. Campbell and Bozorgnia 2014). As such, 𝜏 ln(𝑌) provides an upper bound of the variability due to source effects only, that is the variability one would observe at a given station for repeating event on the same fault. The total variability is given by equation (1-8):

The reported values of the between-event variability are around 0.3 for the natural logarithmic values of the PGA (𝜏 ln(𝑃𝐺𝐴) =0.3). [START_REF] Causse | Are Stress Drop and Rupture Velocity of Earthquakes Independent? Insight from Observed Ground Motion Variability[END_REF] collected recent observations of the between event variability presented in It is important to note that 𝜎 ln(𝑌) has significant impact on Probabilistic Seismic

Hazard Analysis [START_REF] Bommer | Why Do Modern Probabilistic Seismic-Hazard Analyses Often Lead to Increased Hazard Estimates[END_REF]. Therefore, it is of paramount importance to understand and precisely quantify the sources of ground motion variability to improve ground motion prediction for future earthquakes. Studies have found earthquake-magnitude, distance, azimuth and rupture style dependence of 𝜎 ln(𝑌) [START_REF] Abrahamson | Summary of the Abrahamson & Silva NGA Ground-Motion Relations[END_REF][START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF][START_REF] Boore | NGA-West2 Equations for Predicting PGA, PGV, and 5% Damped PSA for Shallow Crustal Earthquakes[END_REF][START_REF] Imtiaz | Is Ground-Motion Variability Distance Dependent? Insight from Finite-Source Rupture Simulations[END_REF]. However, in GMPEs, source effects are, in general, accounted for only by magnitude and a term characterizing the fault mechanism in a simplified form [START_REF] Causse | Constraining the Roughness Degree of Slip Heterogeneity[END_REF]. Due to the improvement of strong motion instrumentation, and consequently the increasing numbers of reliably recorded data and metadata, current 𝐺𝑀𝑃𝐸s can be considered well constrained in the distance range 20 -30 𝑘𝑚 for moderate (𝑀 𝑤 = 6 -6.5) earthquakes. However, data coverage remains sparse for large earthquakes (𝑀 𝑤 > 7), especially at closer distance (𝑅 < 20 𝑘𝑚) [START_REF] Strasser | Sigma: Issues, Insights, and Challenges[END_REF]. Besides, recording stations available in these databases do not fulfill a complete coverage of chancy rupture faults and ground motions close to the faults, where the strongest shaking does occur. The 1992 Landers earthquake provided one of the first set of strong ground motion recordings near a fault [START_REF] Peyrat | Dynamic Modeling of the 1992 Landers Earthquake[END_REF]. Shortly after, the 1994 Northridge and 1995 Kobe earthquakes and then the Izmit and Kocaeli earthquakes in Turkey and the Chi-Chi earthquake in Taiwan in 1999 added more near-source records [START_REF] Halldórsson | Near-Fault and Far-Field Strong Ground-Motion Simulation for Earthquake Engineering Applications Using the Specific[END_REF]. However the total number of recorded near-source ground motions still remains too limited to enable empirically based near-fault ground motion prediction. Therefore near-source ground motion predictions based GMPEs for large earthquakes are highly uncertain.

Numerical ground motion modeling

To overcome the strong ground motion variability and the scarcity of near-fault recordings for large earthquakes in strong ground motion predictions based on GMPEs, numerical modeling of strong ground motion offers an attractive alternative. Therefore, simulation techniques that include complex rupture processes sometimes coupled with complex wave-propagation effects are more and more used to compute and analyze near-field ground motion (e.g. [START_REF] Moschetti | 3D Ground-Motion Simulations of Mw 7 Earthquakes on the Salt Lake City Segment of the Wasatch Fault Zone: Variability of Long-Period (T≥1 s) Ground Motions and Sensitivity to Kinematic Rupture Parameters[END_REF]. The simulations are then based on the above-mentioned representation theorem. Given the importance of the earthquake source process in the observed ground motion and its variability, this raises the question of how to properly model the rupture on an extended fault. The next section is then devoted to the source process.

1.1.4 The source rupture process 1.1.4.1 Omega squared source spectrum [START_REF] Aki | Scaling Law of Seismic Spectrum[END_REF] A given time series can be equivalently expressed by the Fourier amplitude and phase spectra. The Fourier Amplitude Spectra (FAS) then directly shows the frequency content of the time series. Since ground displacement observed in farfield conditions is proportional to the moment rate function, the FAS of ground displacement is directly related to the characteristics of the source. By comparing seismograms obtained by the same seismograph at the same station from two earthquakes with the same epicenter, [START_REF] Aki | Scaling Law of Seismic Spectrum[END_REF]) and [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes[END_REF]) observed that the far-field displacement generated has a spectrum with a constant value at low frequencies, proportional to the seismic moment 𝑀 0 , and proportional to a negative power of two of frequency above a corner frequency 𝑓 𝑐 , defined as the intersection of the low-and high-frequency asymptotes in the spectrum. 𝑓 𝑐 is inversely proportional to the source duration or fault dimension, whence the 𝜔 -2 source model, where 𝜔 = 2𝜋𝑓 is the angular frequency (Figure 12345). [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes[END_REF] theoretically obtained such a simple 𝜔 -2 model for small earthquakes for S-waves, considering an instantaneous shear stress release for a circular crack in a homogeneous medium. Finally, [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes[END_REF]) also proposes a relation between the corner frequency and the source radius, or source dimension, given by equation (1-9):

𝒇 𝒄 = 𝟐. 𝟑𝟒 𝑽 𝒔 𝝅 𝒓 (1-9)
and between the seismic moment, the stress drop and the source radius, given by equation (1-10):

𝑴 𝟎 = 𝟏𝟔 𝟕 ∆𝝉 𝒓 𝟑 (1-10)
Therefore, for a given earthquake, the stress drop is proportional to the corner frequency to the power 3 (relation (1-11)):

∆𝝉 ∝ 𝒇 𝒄 𝟑 (1-11) Chapter 1: State of Art
As the stress pulse is applied instantaneously, no effects of rupture nucleation, propagation or stopping are considered. This point source representation is then a good approximation for stations located far from the fault. But as one get closer to the fault, or for large events for which finite-fault effects are significant, deviations from this simple representation may become important. Note the relation (1-11) stands for a constant speed of the rupture propagation. In order to consider variability in rupture speed, [START_REF] Causse | Are Stress Drop and Rupture Velocity of Earthquakes Independent? Insight from Observed Ground Motion Variability[END_REF] proposed another relation (1-12):

𝒇 𝒄 ∝ 𝑽 𝒓 ∆𝝉 𝟏 𝟑 ⁄
(1-12) More realistically, an earthquake occurs on a finite fault of a rupture length 𝐿: seismic waves are first generated by slip on the fault at the hypocenter, and propagates outward on the fault plane with a rupture speed V r that is typically 60-90% of the shear velocity of the rock [START_REF] Heaton | Evidence for and Implications of Self-Healing Pulses of Slip in Earthquake Rupture[END_REF]. Each point on the fault starts to slip when the moving rupture front arrives at that point, and it takes a finite amount of time (rise time) for that point to undergo slip, by that releasing the total energy described by the seismic moment M 0 . Dynamically, an earthquake results in a strain , assuming an elastic medium. Therefore, the stress drop is related to the energy released as a consequence of an earthquake rupture (expressed by its M 0 ) and to the dimensions of the rupture by: ∆τ = c M 0 L 3 . We describe below various approaches that can be used to model the rupture process.

Dynamic modeling

Earthquake ruptures are not random, but are dictated by physical properties of the faults, friction laws and stresses. So as to model the rupture propagation, one approach is to start from an initial state of stress on the fault plane (having some hypothesis about the material properties around the source, and the initial and boundary conditions on the fault plane), force the rupture to initiate at a particular point, and then let the rupture evolve freely. The distribution of slip on the fault plane is obtained by solving the elasto-dynamic equation of motion under a frictional failure model, considering essentially the energy balance at the crack tip during rupture growth. This is called dynamic rupture modeling [START_REF] Madariaga | High Frequency Radiation from Dynamic Earthquake[END_REF][START_REF] Pulido | Estimation of the High-Frequency Radiation of the 2000 Tottori (Japan) Earthquake Based on a Dynamic Model of Fault Rupture: Application to the Strong Ground Motion Simulation[END_REF][START_REF] Madariaga | Earthquake Dynamics on Circular Faults: A Review 1970-2015[END_REF]. Therefore, dynamic modeling requires knowing the initial stresses as well as the time varying stress and friction on the fault surface during the earthquake rupture process.

One of the first dynamic rupture model was proposed by [START_REF] Madariaga | High Frequency Radiation from Dynamic Earthquake[END_REF], who performed dynamic finite-difference simulations for a circular shear crack nucleating at its center and propagating at a constant rupture velocity until it suddenly stops at a given radius. The displacement amplitude spectra from the Madariaga's model also roughly shows 𝜔 -squared as observed by [START_REF] Aki | Scaling Law of Seismic Spectrum[END_REF]), but the corner frequencies are a factor of two smaller than expected from Brune's considerations. This illustrates that the relation between the corner frequency and the other source parameters is model dependent.

For many years, seismologists have modeled earthquake fault motion by a homogeneous slip and stress drop over the entire fault plane. After, they recognized the need for an irregular slip motion through the use of patches (asperities and barriers-Figure 123456). We expect that some of the patches on a fault behave as asperities and others as barriers. Asperities are patches where stress has accumulated before the earthquake and are characterized by a high stress drop. They explain the occurrence of the main shocks. Barriers are unbroken patches after an earthquake. They explain the occurrence of the aftershocks [START_REF] Aki | Asperities, Barriers, Characteristic Earthquakes and Strong Motion Prediction[END_REF]. The static and the dynamic stress drop may then vary greatly over the fault. As such, several studies have proposed to perform dynamic rupture simulations considering heterogeneous distributions of stress and friction parameters. Yet, the physical parameters governing such models (stress and strength distribution along space and time on the fault plane) are not well constrained.

Furthermore, these simulations are computationally expensive (in terms of CPU requirements) and cover a rather limited frequency range.

Kinematic and Pseudo-Dynamic modeling

Because friction and stress parameters are often hardly available and because dynamic modeling is computationally expensive, we refer to kinematic source model, which is another approach to describe the source. It consists of a priori prescribing the displacement discontinuity across the fault surface. The local sourcetime function needs to be specified [START_REF] Tinti | A Kinematic Source-Time Function Compatible with Earthquake Dynamics[END_REF][START_REF] Liu | Prediction of Broadband Ground-Motion Time Histories; Hybrid Low/High-Frequency Method with Correlated Random Source Parameters[END_REF] to describe the evolution of the slip with time by kinematic parameters: the attributed final slip for each fault point 𝐷, the time needed to reach the maximum slip (called the rise time 𝑇 𝑟𝑖𝑠𝑒 ) and the rupture velocity 𝑉 𝑟 to determine when each fault point ruptures. The simplest kinematic source model was developed by Haskell [START_REF] Haskell | Total Energy and Energy Spectral Density of Elastic Wave Radiation from Propagating Faults. Part II. A Statistical Source Model[END_REF]. It consists of representing the fault area as a rectangular fault plane of length 𝐿 and width 𝑊 , with 𝐿 ≫ 𝑊 . The rupture propagates unilaterally, assuming constant rupture velocity, rise-time and final slip. This simple model is consistent with the 𝜔 -2 model proposed by [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes[END_REF].

In this context, pseudo-dynamic means that the kinematic rupture process, quantified in terms of slip, slip rate, rise time, and rupture speed at each point of the fault are specified according to some of the principles of earthquake dynamics. The advantage of pseudo-dynamic (PD) models is to maintain the computational efficiency of kinematic models, preserving at the same time enough degree of freedom to represent models that are compatible with the physics of a dynamic rupture [START_REF] Mai | Toward a More Physical Basis for Strong-Motion Simulation[END_REF][START_REF] Guatteri | A Pseudo-Dynamic Approximation to Dynamic Rupture Models for Strong Ground Motion Prediction[END_REF][START_REF] Mena | Pseudodynamic Source Characterization for Strike-Slip Faulting Including Stress Heterogeneity and Super-Shear Ruptures[END_REF][START_REF] Schmedes | A Kinematic Rupture Model Generator Incorporating Spatial Interdependency of Earthquake Source Parameters[END_REF][START_REF] Schmedes | A Kinematic Rupture Model Generator Incorporating Spatial Interdependency of Earthquake Source Parameters[END_REF][START_REF] Song | Developing a Generalized Pseudo-Dynamic Source Model of Mw 6.5-7.0 to Simulate Strong Ground Motions[END_REF].

Despite the recent advancement in kinematic modeling, a major issue for simulating future earthquakes remains the choice of source parameters. The next section describes in details the different source parameters and the main observations that can help constraining them, and the major source of uncertainties.

Constraining source parameters from observations

The simple source models with uniform slip and slip velocity are a very strong simplification of the rupture process during earthquakes. Indeed, heterogeneities exist. The heterogeneity can originate from a heterogeneous medium, or nonuniform stresses from creeping fault regions and strength along the fault surface, or from non-planar fault geometry [START_REF] Ampuero | Heterogeneities and Complexity in Earthquake Dynamics[END_REF][START_REF] Irikura | Recipe for Predicting Strong Ground Motion from Crustal Earthquake Scenarios[END_REF]. In the following, we will then consider the source parameters at two different scales: Large-scale: Describes the average properties of the fault rupture. It comprises the macroscopic geometrical parameters (the rupture length L and width W) as well as the average values on the fault of stress drop ∆𝜏, slip D, rupture speed Vr and rise time Trise.

Local-scale: Describes the spatial heterogeneities of the source parameters along the fault area. It comprises the microscopic geometrical variations of the fault surface, as well as the spatial variation of ∆𝜏, D, Vr, Trise along the fault during the rupture.

Large-scale and local-scale processes involved in earthquake rupture are probably not independent. For instance, the heterogeneity of fault systems (roughness of fault surface topography, heterogeneity of source parameters, etc.) has a large impact on the rupture process, not only at local-scale (slip and rupture fluctuations…), but also at large scale (final rupture dimension, average stress drop, average rupture velocity). (Candela, Renard, Bouchon, et al. 2011) proposed that a controlling parameter of the average stress drop is related to the scale properties of the topography of the fault surface (i.e., fault roughness). Furthermore, [START_REF] Bouchon | Faulting Characteristics of Supershear Earthquakes[END_REF]) observed that simple geometry and homogeneity of faults lead to faster rupture speed. Also, [START_REF] Bydlon | Rupture Dynamics and Ground Motions from Earthquakes in 2-D Heterogeneous Media[END_REF] conducted 2D dynamic rupture simulations in a heterogeneous medium and observed that heterogeneities of the material structure decrease the average slip, rupture length (and therefore the magnitude) and rupture speed along the fault interface. In addition, it increases the fluctuations in slip and rupture speed. In case of strong heterogeneities, the rupture would stop. Likewise, [START_REF] Zielke | Fault Roughness and Strength Heterogeneity Control Earthquake Size and Stress Drop[END_REF] showed that earthquake's moment release and stress drop may vary widely depending on the geometric roughness of the rupture surface and the location of strength asperities (Figure 1234567(1) rough faults release less seismic moment per tectonic loading than smooth ones;

(2) ruptures with near-center strength asperities increase the seismic moment relative to the homogeneous strength case. Near-edge strength asperities have the opposite effect. A quantitative understanding of rupture complexity is essential for ground motion prediction studies that require an adequate characterization of the source parameters to generate realistic scenario earthquakes. Ideally, our understanding of earthquake source parameters and its complexity stems from surface rupture observation, finite source rupture models of past earthquakes or other seismological observations, and more recently from dynamic rupture simulations. Those techniques identify distribution characteristics of source parameters, scalingrelationships and cross-correlation relationships between source parameters, which are essential to design realistic rupture scenarios of potential future earthquakes.

Fault observation measurements

The ruptures of large earthquakes (𝑀 𝑤 > 6.5 -7) generally reach the free surface, and can then be directly observed on field. By compiling a large number of surface rupture observations, coupled with other techniques, several scaling relationships have been published to characterize the distribution of large-scale source parameters like the average slip 𝐷 or the rupture length 𝐿 [START_REF] Wells | New Empirical Relationships among Magnitude, Rupture Length, Rupture Width, Rupture Area, and Surface Displacement[END_REF]. More recent observations provide, for a few earthquakes, high-resolution profiles of the slip along the fault. For instance, [START_REF] Rockwell | Surface Rupture and Slip Distribution of the 1940 Imperial Valley Earthquake, Imperial Fault, Southern California; Implications for Rupture Segmentation and Dynamics[END_REF] observed the variability of the surface displacement over a fault (Figure 12345678). They analyzed aerial photography and provided a detailed distribution of the surface slip over 14 km for the 1940 Imperial Valley earthquake that ruptures with a magnitude of 7 producing nearly 60 km of surface rupture. This observation shows the variability of the slip over the fault. According to this distribution, the standard deviation is at the same order as the mean value. Large-scale geometrical features of fault surface ruptures may also help constraining the rupture properties of future events. [START_REF] Sagy | Evolution of Fault-Surface Roughness with Slip[END_REF] investigated the relationship between slip-surface roughness and fault accumulated fault displacement using laser-based methods to map exposed fault surfaces (e.g. Figure 1-9) and showed that fault-surface roughness evolves with increasing fault displacement. Slip surfaces of faults which have accumulated little displacement are relatively rough at all measured scales, whereas those of large-displacement faults are polished at local scales but contain bumps and depressions at scales of a few to several meters. As above-mentioned, [START_REF] Bouchon | Faulting Characteristics of Supershear Earthquakes[END_REF]) observed that fault segments with a simple geometry are associated with faster rupture speed, exceeding the shear wave speed (super-shear ruptures). Furthermore, the surfaces of some faults are exhumed, providing some direct observations of the fault topography. With the recent development of high-resolution distance meters, geometrical fault roughness is observable on the fault surface and can be accurately quantified using statistical approaches (Candela, Renard, Schmittbuhl, et al. 2011) -

Insights from seismological observations

While surface ruptures allow for direct observation of some rupture parameters, they do not provide information about the rupture in depth. This can be achieved from seismological observations.

Finite source rupture inversion models

Using strong ground motion recordings often enhanced with geodetic and/or tsunami measurements, finite source kinematic-inversion models provide the space-time evolution of co-seismic displacement that occurs on one or more fault segments. The resulting finite fault rupture models (also known as slip models) quantify the heterogeneous distribution of slip (and hence the static stress drop) and rise time (and thus, in combination with slip, the slip velocity), and how fast the rupture expanded over the fault surface (P.M. Mai et al. 2016). As inversion methods become more advanced, more real-time data is available, and computer capacities are enhanced, we see growing numbers of finite-fault rupture models. Models obtained for a given earthquake being by different scientific teams are however often remarkably dissimilar. The discrepancies may be attributed to 1) the differences in data selection and processing, 2) the methods used for computing the Green's functions for each dataset 3) imperfect knowledge of the Earth structure and fault geometry and 4) the method and parameterization for the inversion itself (linearized or fully nonlinear inversion; spatial and temporal discretization; applied smoothing and regularization). The source inversion problem is then non-unique (Figure 1 -11). Despite uncertainties in the inversion problem, finite-fault source models attempt to produce credible images of earthquake rupture processes by achieving consistency between observed data and geophysical model predictions up to ~1Hz [START_REF] Zhang | Pseudo-Dynamic Earthquake Source Inversion[END_REF][START_REF] Mai | Standards for Documenting Finite-Fault Earthquake Rupture Models[END_REF]). Thus, a database containing now a few hundreds of these finite source rupture models was compiled (http://equake-rc.info/SRCMOD/). Finite-fault rupture models, combined with field observations of surface rupture length and data from space distribution of aftershocks allow seismologists to study relations between the parameters of the subsurface faults. Thus, several source-scaling relationships between the observable large-scale source parameters (fault length, fault area, fault slip), measurable localscale source parameters (e.g. the correlation lengths along dip and strike directions -Figure 1-12) and rupture magnitude have then been published based on statistical observations [START_REF] Wells | New Empirical Relationships among Magnitude, Rupture Length, Rupture Width, Rupture Area, and Surface Displacement[END_REF][START_REF] Somerville | Characterizing Crustal Earthquake Slip Models for the Prediction of Strong Ground Motion[END_REF][START_REF] Papazachos | Global Relations between Seismic Fault Parameters and Moment Magnitude of Earthquakes[END_REF][START_REF] Goda | New Scaling Relationships of Earthquake Source Parameters for Stochastic Tsunami Simulation[END_REF][START_REF] Thingbaijam | New Empirical Earthquake Source-Scaling Laws[END_REF]. In addition, finite fault source inversions reveal the spatial complexity of earthquake slip over the fault plane. Images of the spatial and temporal evolution of earthquake slip on fault planes provide compelling evidence that fault displacement is spatially variable at all resolvable scales. [START_REF] Mai | A Spatial Random Field Model to Characterize Complexity in Earthquake Slip[END_REF] proposed a characterization of spatial complexity of earthquake slip as found in finite source slip inversions that follows the von karman autocorrelation function, with a decay of about 𝑘 -2 in the wavenumber domain. In this distribution, the slip correlation lengths increase with increasing earthquake magnitude. The resulting slip spectrum is in accord with previous theoretical slip models: [START_REF] Andrews | A Stochastic Fault Model: 1. Static Case[END_REF][START_REF] Andrews | A Stochastic Fault Model. 2. Time-Dependent Case[END_REF] showed that a slip spectrum that decays as 𝑘 -2 in the wavenumber domain leads to far-field displacements that follow the widely observed 𝜔 -2 spectral decay, assuming that stress drop is scale invariant. Based on this concept, [START_REF] Bernard | Slip Heterogeneity, Body-Wave Spectra, and Directivity of Earthquake Ruptures[END_REF] introduced the 𝑘 -2 model, in which the slip spectrum decays as 𝑘 -2 beyond a corner wavenumber, which is inversely proportional to the fault length. In this representation, slip is scale-invariant. Hence there are no characteristic length scales of small size of asperities. In their model, the rupture velocity is constant. Thus [START_REF] Hisada | A Theoretical Omega-Square Model Considering the Spatial Variation in Slip and Rupture Velocity[END_REF] proposed an adapted 𝑘 -2 source model by considering spatial variation not only in final slip, but also in rupture velocity. This model is also consistent with the commonly observed 𝜔 -2 model of displacement spectra.

Earthquake displacement spectra

Many authors looked at the amplitude spectra of the recorded ground motion in order to extract conclusions related to the source parameters. [START_REF] Allmann | Global Variations of Stress Drop for Moderate to Large Earthquakes (English)[END_REF] analyzed the corner frequencies values 𝑓 𝑐 obtained from Brune-type source model matching the 𝜔 -2 model. 𝑓 𝑐 scales with the seismic moment 𝑀 0 , and is generally directly related to the stress drop, and sometimes the rupture speed (equations (1-11) and (1-12)).

Apparent Source time functions

Another way to characterize source parameters is to analyze the moment rate function 𝑀 ̇(𝑡) (also called source time functions), which can provide information about integrated source parameters like source duration or stress drop. By studying apparent moment rate functions (that is moment rate functions observed at different stations around fault ruptures), one can also get information about the rupture propagation such as the average rupture velocity. This approach was conducted for example by [START_REF] Chounet | Global Catalog of Earthquake Rupture Velocities Shows Anticorrelation between Stress Drop and Rupture Velocity[END_REF]) who concluded that the distributions of stress drop and rupture velocity are not independent, but anticorrelated. Besides, [START_REF] Archuleta | Moment Rate Scaling for Earthquakes 3.3 ≤ M ≤ 5.3 with Implications for Stress Drop[END_REF]) also looked at the apparent moment rate functions and observed that their Fourier amplitude spectra must be characterized not by a single but two corner frequencies to match observations of PGA and PGV values.

Insight from dynamic simulations

Since in-situ measurements directly at the earthquake source are non-existent, seismologists often make use of dynamic simulations to get insights on the rupture source.

Link between frictional and kinematic rupture parameters

With dynamic simulations, seismologists have analyzed the relationship between the stresses and strengths along the fault, and the slip and the rupture velocity during the earthquake. [START_REF] Madariaga | High Frequency Radiation from Dynamic Earthquake[END_REF][START_REF] Pulido | Estimation of the High-Frequency Radiation of the 2000 Tottori (Japan) Earthquake Based on a Dynamic Model of Fault Rupture: Application to the Strong Ground Motion Simulation[END_REF][START_REF] Mena | Pseudodynamic Source Characterization for Strike-Slip Faulting Including Stress Heterogeneity and Super-Shear Ruptures[END_REF][START_REF] Madariaga | Earthquake Dynamics on Circular Faults: A Review 1970-2015[END_REF] dynamically simulated ground motion and showed that geometrical obstacles or barriers (that is area of high fault strength) cause abrupt changes in rupture velocity. Besides, [START_REF] Bouchon | The State of Stress on Some Faults of the San Andreas System as Inferred from Near-Field Strong Motion Data[END_REF] observed that the regions of the fault that break with a high stress drop are also the regions where slip is large. He also related the change of the rupture speed to the strength excess: a high strength excess value implies that either the fault strength was high or the shear stress was low at that specific area of the fault prior to the earthquake. Thus, the spatial distribution of the strength excess over the fault plane appears to be inversely correlated to the local rupture velocity. Rupture propagates slowly where the relative fault strength is high and accelerates over low strength regions. In addition, dynamic simulations can help representing more realistic ruptures in kinematic simulations: [START_REF] Mai | Accounting for Fault Roughness in Pseudo-Dynamic Ground-Motion Simulations[END_REF]) conducted a suite of dynamic earthquake simulations for various rough-fault realizations and suggested to take into account fault roughness in kinematic models behavior by retaining the moment tensor orientations but neglecting their off-fault positions.

1. 1.5.3.2 Correlation between kinematic rupture parameters Last but not least, simulating dynamic source models and collecting the resulting source parameters in a catalogue is a way to study statistical relationships between various rupture parameters (e.g. [START_REF] Song | Physics-Based Earthquake Source Characterization and Modeling with Geostatistics[END_REF]. [START_REF] Schmedes | A Kinematic Rupture Model Generator Incorporating Spatial Interdependency of Earthquake Source Parameters[END_REF] analyzed a database of dynamic strike-slip rupture models computed using different models of initial conditions of stresses and strength properties, and concluded that slip strongly correlates with rise time. Whereas they reported that the correlation pattern between the rupture velocity and the slip is unclear at least under certain conditions of dynamic rupture models. [START_REF] Bizzarri | Rupture Speed and Slip Velocity; What Can We Learn from Simulated Earthquakes?[END_REF]) considered a wide catalogue of dynamic 3D models, and observed that the peak slip rate and the rupture speed are positively correlated, which is supported by theory [START_REF] Ida | The Maximum Acceleration of Seismic Ground Motion[END_REF]. [START_REF] Oglesby | Stochastic Fault Stress; Implications for Fault Dynamics and Ground Motion[END_REF] used 3D dynamic models with variable assumptions on strength and stress heterogeneity and concluded that rupture velocity, rise time, and slip are associated with the fault strength and stress drop, as well as each other, however the connections between these quantities are not simple. On the other hand, [START_REF] Trugman | A 2D Pseudodynamic Rupture Model Generator for Earthquakes on Geometrically Complex Faults[END_REF]) presented a 2D pseudodynamic model that emulates earthquake source parameters on rough faults where final slip, local rupture speed and peak slip velocity are anti-correlated with the observed-fault roughness. Thus, the correlation patterns between the different kinematic rupture parameters remain unclear. This is mainly because dynamic simulations are sensitive to the assumptions on the input friction laws, which are poorly constrained. [START_REF] Song | Developing a Generalized Pseudo-Dynamic Source Model of Mw 6.5-7.0 to Simulate Strong Ground Motions[END_REF] showed that the source parameter correlation structures can be significantly affected by the input fracture energy distribution (Figure 1 -15).

A compilation of references to papers that studied the correlation between source parameters at different scales, between kinematic and dynamic characteristics, is presented in Table 12 1.1.6 From source rupture towards ground motion These local and large-scale rupture features largely affect the ground motion. The waves transport the complexity of the rupture process through the Earth structure to the ground surface in terms of displacement, velocity or acceleration ground motion. Each source parameter induces its own signature on the ground motion, in various frequency ranges. For instance, ground motion is strongly affected by the fault dimension for a given moment magnitude [START_REF] Aagaard | Characterization of Near-Source Ground Motions with Earthquake Simulations[END_REF], the hypocenter position [START_REF] Aagaard | Characterization of Near-Source Ground Motions with Earthquake Simulations[END_REF][START_REF] Somerville | Modification of Empirical Strong Ground Motion Attenuation Relations to Include the Amplitude and Duration Effects of Rupture Directivity[END_REF][START_REF] Ripperger | Variability of Near-Field Ground Motion from Dynamic Earthquake Rupture Simulations[END_REF], rupture speed and rise time especially in the fault vicinity [START_REF] Aagaard | Characterization of Near-Source Ground Motions with Earthquake Simulations[END_REF]. Recent near-fault ground motion simulation studies (e.g. [START_REF] Moschetti | 3D Ground-Motion Simulations of Mw 7 Earthquakes on the Salt Lake City Segment of the Wasatch Fault Zone: Variability of Long-Period (T≥1 s) Ground Motions and Sensitivity to Kinematic Rupture Parameters[END_REF] all show that earthquake ground motions and its variability are highly sensitive to the choice of slip distribution, rupture speed, slip velocities and hypocenter locations. It shows that there is a real need for further characterization of the kinematic source parameters for probabilistic seismic hazard analyses. 1.2 Impact of seismic rupture on surface ground motion

Understanding and mitigating earthquake risk depends critically on predicting the intensity of strong ground motion, including estimates of the aleatory variability, which remains a scientific challenge. Here, the term ground motion variability refers to the variability due to source effects only, that is the variability one would expect for repeating events of the same magnitude on a given fault, recorded at the same station. The fault rupture process that generates seismic waves is complex and incompletely understood. In this chapter, we aim to better quantify the link between the rupture properties and the Peak Ground Acceleration (PGA) and the Peak Ground Velocity (PGV), which are two commonly used measures of the ground motion intensity. Note that these two quantities are controlled by different frequency ranges. As explained in section 1.1.4.1, the frequency content of ground motion is, at the first order, controlled by the corner frequency fc. The ground displacement amplitude spectrum decays with a slope of -2 for frequencies larger than fc. Hence, the velocity spectrum increases with the frequencies up to fc and then decreases as f -1 , and the acceleration spectrum is flat above fc. This implies that PGA is essentially controlled by frequencies larger than fc, while the PGV is expected to be driven by lower frequencies (note that reported values of fc for various moment magnitudes are displayed on Figure 1-13). A more detailed analysis of the frequency range that mostly controls PGA and PGV values is proposed in chapter 2.

1.2.1 PGA controlled by the large-scale source parameters: role of average stress drop 1.2.1.1 Stress parameter controls the PGA The stress drop ∆τ has become a key parameter to measure the strength of the observed high-frequency ground motion. Many authors referred to stress drop Δτ and its variability driving the high-frequency ground motion: [START_REF] Lavallée | Coupling of the Random Properties of the Source and the Ground Motion for the 1999 Chi Chi Earthquake[END_REF] noticed that the variability of the PGA, which represents a measure of the high-frequency ground motion, is not that different from the variability of the Δτ. [START_REF] Cotton | What Is Sigma of the Stress Drop?[END_REF] identified that the PGA variability should translate directly into earthquake Δτ variability as described by equation (1-13), derived under the assumption of the classical omega squared model and the random vibration theory:

𝝈 𝐥𝐧(∆𝝉) = 𝟔 𝟓 𝝈 𝐥𝐧(𝑷𝑮𝑨)
(1-13) (Bindi, Spallarossa, and Pacor 2017a) developed local GMPE models from records in Central Italy and observed a clear dependence of the between-event residuals of PGA on the stress drop. [START_REF] Youngs | Magnitude-Dependent Variance of Peak Ground Acceleration[END_REF] suggested the reduction of the stress drop variability with increasing magnitude as possible explanation for the reduction of the PGA between event variability with increasing magnitude. [START_REF] Oth | On the Relation of Earthquake Stress Drop and Ground Motion Variability[END_REF]) studied a large dataset from Japan with 2 different families of stress drop values, computed from the corner frequency, and corresponding to 2 different regions. They found a clear correlation between the PGA between event residuals and the stress drop estimates (Figure 1-16). However, the determination of stress drop is not an easy task. Assuming an extended fault, ∆τ ∝ M 0 L 3 , this means that the estimation of the static stress drop is related to the estimation of source dimensions, which, in turn, depends on the methodology being adopted. The assumption of a Brune's source model [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes[END_REF] is the most commonly used assumption (section 1.1.5.2.2). The stress drop is then assumed to be proportional to the cube of the corner frequency ∆τ ∝ fc 3 . Thus the uncertainty in stress drop estimates is dominated by the uncertainty in the corner frequency fc. fc values are generally obtained from Brune's spectrum matching. This parameter is affected by measurement errors and possible errors due to the imperfect knowledge of path and site effects, which might be quite large too. In addition, the utilized model to relate fc and the source dimension must always be clearly stated. Average stress drop values can also be derived from static slip distributions, obtained from kinematic source inversions, from which the distribution of static stress drop can be inferred. It should be emphasized that stress drop is spatially variable (section 1.1.5), and hence any estimate of the average stress drop depends on the adopted definition of average. For instance, it makes of course a strong difference whether stress drop averaged over the asperity area or stress drop averaged over the entire fault plane is considered: the stress drop averaged over the entire fault plane underestimates the stress drop on the asperities. [START_REF] Zielke | Fault Roughness and Strength Heterogeneity Control Earthquake Size and Stress Drop[END_REF] shed the light on the question of how much of the variability from classical 𝑀 0 -𝑓 𝑐 is actually representative of true source variability.

Are stress drop and stress parameters equivalent?

Finally, it is important to note that in the framework of high-frequency ground motion simulation, the physical meaning of stress drop remains ambiguous. In the context of generation of acceleration time histories using stochastic methods, several authors refer to the stress parameter as the parameter controlling the highfrequency ground motion level [START_REF] Atkinson | Ground-Motion Relations for Eastern North America[END_REF]. In such simulations, the stress parameter is often used as if it were stress drop. Although originally derived from a relation between static stress drop, fault slip and fault size, stress parameter is best thought of as simply a parameter controlling the strength of the highfrequency radiation [START_REF] Boore | Stochastic Simulation of High-Frequency Ground Motions Based on Seismological Models of the Radiated Spectra[END_REF]. Static stress drops obtained from seismological spectra range from about 0.1 to 100 MPa and is function of the magnitude (Figure 1-13) [START_REF] Allmann | Global Variations of Stress Drop for Moderate to Large Earthquakes (English)[END_REF], however the stress parameter used in the stochastic models is about 10 MPa and is independent of the magnitude (Figure 1-17) [START_REF] Edwards | Measurements of Stress Parameter and Site Attenuation from Recordings of Moderate to Large Earthquakes in Europe and the Middle East[END_REF].

On the other hand, recent studies shows that using distributions of stress drop obtained from 𝑓 𝑐 measurements in stochastic simulations can lead to significant overestimation of the PGA variability [START_REF] Cotton | What Is Sigma of the Stress Drop?[END_REF][START_REF] Oth | On the Relation of Earthquake Stress Drop and Ground Motion Variability[END_REF]. [START_REF] Causse | Are Stress Drop and Rupture Velocity of Earthquakes Independent? Insight from Observed Ground Motion Variability[END_REF] explained this discrepancy by highlighting on the role of the average rupture velocity. Rupture velocity may be anti-correlated with the average stress drop, which tends to decrease the variability of PGA and makes it closer to the observations of between event variability of PGA. Recently, [START_REF] Archuleta | Moment Rate Scaling for Earthquakes 3.3 ≤ M ≤ 5.3 with Implications for Stress Drop[END_REF] proposed that the stress parameter and the stress drop are simply not the same quantity. They propose that the PGA is not controlled by the corner frequency (and hence by the stress drop) but by source processes at a larger frequency, represented by a second corner frequency. The latter may be due to the breakdown of asperities smaller than the overall source dimension, which excite high frequency ground motion.

PGA controlled by the local-scale heterogeneities

The corner frequency parameter connected to ∆τ and Vr might not be sufficient to relate the source process to the generation of high-frequency ground motion. [START_REF] Song | Developing a Generalized Pseudo-Dynamic Source Model of Mw 6.5-7.0 to Simulate Strong Ground Motions[END_REF][START_REF] Song | Pseudo-Dynamic Source Modelling with 1-Point and 2-Point Statistics of Earthquake Source Parameters[END_REF]) used a pseudo-dynamic source model, where source parameters are calibrated from a database of dynamic source models: they run a sensitivity analysis and showed that, for the same ∆τ and Vr, both larger standard deviation and correlation between source parameters produce stronger peak ground velocities near the source.

1.2.2.1 Insight from dynamic simulations [START_REF] Madariaga | High Frequency Radiation from Dynamic Earthquake[END_REF] dynamically simulated ground motion and showed that strong variations of the rupture velocity at the crack boundaries (generating the stopping phases) play a very important role in the radiation of high frequency from the source. High frequency energy is emitted only when the rupture front is accelerating or decelerating. Rupture velocity changes are largely induced by barriers (locally stronger fault sections) across the fault plane, and high frequency radiation mainly originates within asperities (large stress drop regions). Generally, high-frequency energy generates in areas where the product of dynamic stress drop and rupture velocity changes is maximum [START_REF] Pulido | Estimation of the High-Frequency Radiation of the 2000 Tottori (Japan) Earthquake Based on a Dynamic Model of Fault Rupture: Application to the Strong Ground Motion Simulation[END_REF]. Likewise, with dynamic simulations, seismologists can observe the importance of the geometrical complexities on the ground motion. For instance, [START_REF] Dunham | Earthquake Ruptures with Strongly Rate-Weakening Friction and Off-Fault Plasticity, Part 2: Nonplanar Faults[END_REF][START_REF] Bydlon | Rupture Dynamics and Ground Motions from Earthquakes in 2-D Heterogeneous Media[END_REF][START_REF] Zielke | Fault Roughness and Strength Heterogeneity Control Earthquake Size and Stress Drop[END_REF] showed the importance of geometric roughness on the rupture process and resulting high-frequency ground motions in the near-fault region. They performed 2D dynamic rupture simulations with respect to geologic observations of fault surface roughness. The observed levels of roughness introduce variations in slip and rupture velocity in a manner consistent with realistic high-frequency ground motions. (Archuleta and Crempien 2015) simulated pseud-dynamic model and performed a sensitivity analysis of the ground motion to the slip correlation lengths using kinematic rupture simulations: they simulated synthetic white noise scaled to the seismic moment over the fault plane, and filtered it in the wavenumber domain with a Von Karman power spectrum based on the values of the slip correlation lengths. They found that as the asperity size decreases, the ground motion amplitude and its variability decrease (Figure 1 -19). This is in contradiction with [START_REF] Causse | New Approach for Coupling K-2 and Empirical Green's Functions: Application to the Blind Prediction of Broad-Band Ground Motion in the Grenoble[END_REF], who obtained that large correlation length of the slip tend to generate smaller PGA values. It shows the need for further investigation on the link between the heterogeneity level of kinematic rupture parameters and the high-frequency ground motion PGA.

Insight from kinematic simulations

THE SOURCE PARAMETERS CONTROLLING THE HIGH-FREQUENCY GROUND MOTION

In this chapter we investigate the source parameters including average source parameters and heterogeneities that mostly contribute to the generation of the peak ground motions values. This Manuscript is submitted for publication to: Pure and Applied Geophysics (PAGEOPH).

- -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------2.1 Abstract During an earthquake, seismic waves carry the complexity of the rupture to the ground surface. Empirical ground motion prediction equations, calibrated by past earthquake seismic recordings, are often used to predict Peak Ground Acceleration (PGA) and its variability. However, the scarcity of near-fault recordings for large earthquakes prevents using such equations to predict near-fault PGA. Simulation of strong ground motion then offers an attractive alternative to assess seismic hazard in near-fault. In order to better understand the effects of rupture parameters on surface ground motion and to capture the key source ingredients that most impact ground motion variability, we simulate ground motions produced by various M7 rupture earthquake scenarios on a vertical strike-slip fault. We compute the ground motion up to 5 Hz at sites located at 5 km, 25 km and 70 km from the fault. The kinematic rupture parameters are modelled using a statistical source model generator as proposed in [START_REF] Song | Pseudo-Dynamic Source Modelling with 1-Point and 2-Point Statistics of Earthquake Source Parameters[END_REF]. We show that PGA is mainly generated by abrupt changes of the rupture propagation, that is, stopping phases at the fault boundaries or strong heterogeneities of rupture speed along the rupture. We observe that PGA is mostly controlled by the average rupture speed and the average stress drop (in the far-field) and to a lesser extent by the standard deviation of the rupture speed. Interestingly, correlation between source parameters and spatial correlation length do not affect average PGA and related variability significantly, for the set of stations in study.

Key words: source simulation, near-fault, far-fault, sensitivity, peak ground acceleration.

Introduction

Earthquakes have caused, and can cause in the future, enormous loss of life, injury, destruction of property, and economic and social disorder. Determining the earthquake intensity does help developers make decisions about where to build and what type of forces structures should withstand. Statistical analysis of strong ground motion databases helps providing quantitative estimates of expected ground-motion levels for a potential future earthquake. Due to the lack of recordings in the vicinity of faults, there is however the need to develop physicsbased simulation techniques incorporating the complexity of earthquake rupture to obtain reliable near-field ground motions (e.g. [START_REF] Moschetti | 3D Ground-Motion Simulations of Mw 7 Earthquakes on the Salt Lake City Segment of the Wasatch Fault Zone: Variability of Long-Period (T≥1 s) Ground Motions and Sensitivity to Kinematic Rupture Parameters[END_REF]. The present study focuses on the relationship between the rupture process and the highfrequency ground motion (average and variability) represented by the Peak Ground Acceleration (PGA) and the Peak Ground Velocity (PGV).

It is commonly claimed that the PGA is driven by the stress drop Δτ (e.g. [START_REF] Bindi | Between-Event and between-Station Variability Observed in the Fourier and Response Spectra Domains: Comparison with Seismological Models[END_REF][START_REF] Cotton | What Is Sigma of the Stress Drop?[END_REF][START_REF] Oth | On the Relation of Earthquake Stress Drop and Ground Motion Variability[END_REF], which is related to the available elastic energy during the rupture process on the rupture area. The stress drop is commonly supposed to be proportional to the cube of the corner frequency f 𝑐 , determined from the Fourier amplitude spectra of the displacement ground motion, under the assumption of a 𝜔 -2 source model [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes[END_REF]). In the framework of stochastic simulations of ground motion time histories, the 'stress parameter' is commonly used as a proxy for the high-frequency level of ground motion (e.g. [START_REF] Boore | Stochastic Simulation of High-Frequency Ground Motions Based on Seismological Models of the Radiated Spectra[END_REF][START_REF] Edwards | Measurements of Stress Parameter and Site Attenuation from Recordings of Moderate to Large Earthquakes in Europe and the Middle East[END_REF]. The 'stress parameter' is generally used as if it would be the stress drop. Using the random vibration theory, [START_REF] Boore | Stochastic Simulation of High-Frequency Ground Motions Based on Seismological Models of the Radiated Spectra[END_REF]) obtained that:

𝑃𝐺𝐴 ∝ Δτ 0.8 (2-1)
for a [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes[END_REF] source model. [START_REF] Oth | On the Relation of Earthquake Stress Drop and Ground Motion Variability[END_REF]) observed that PGA is effectively strongly correlated by the Brune's stress drop. Nevertheless, the PGA variability which is predicted from the Brune' stress drop distribution using equation (2-1) has been shown to be larger than the observed between-event variability of PGA [START_REF] Cotton | What Is Sigma of the Stress Drop?[END_REF][START_REF] Oth | On the Relation of Earthquake Stress Drop and Ground Motion Variability[END_REF]. This points out that the 'stress parameter' and the stress drop are two different physical quantities. Recently, [START_REF] Causse | Are Stress Drop and Rupture Velocity of Earthquakes Independent? Insight from Observed Ground Motion Variability[END_REF] pointed out the importance of considering the rupture velocity V 𝑟 in the corner frequency definition, and obtained:

𝑃𝐺𝐴 ∝ Δτ 0.8 V 𝑟 2.4
(2-2) [START_REF] Causse | Are Stress Drop and Rupture Velocity of Earthquakes Independent? Insight from Observed Ground Motion Variability[END_REF] suggested that anti-correlation between stress drop and rupture velocity may be a potential source mechanism that reduces the PGA variability and makes it closer to observations of between-event variability of PGA. Such anti-correlation has been recently observed from a database of apparent source time functions including ~100 shallow earthquakes (Chounet et al. 2018).

Note that Equation (2-2) is based on a simple point source representation, and Δτ and V 𝑟 are then considered as average source properties. As such, the PGA observed in far-field is mainly controlled by the corner frequency, which, in turn, depends on large-scale source parameters describing the macroscopic features of the rupture process.

On the other hand, several studies evidence that the high-frequency ground motion is controlled by smaller-scale processes at frequencies larger than the corner frequency. Recently, [START_REF] Archuleta | Moment Rate Scaling for Earthquakes 3.3 ≤ M ≤ 5.3 with Implications for Stress Drop[END_REF] proposed that earthquake momentrate functions are better fitted by a two-corner frequency spectrum model than by a classical Brune's source model. The largest corner frequency is associated with localscale source parameters and may drive the PGA. Which of those source processes mainly control the high-frequency ground motion remains however strongly debated. [START_REF] Madariaga | High Frequency Radiation from Dynamic Earthquake[END_REF] showed that strong variations of the rupture velocity are very efficient sources of high frequency radiation, especially rupture stopping phases at the fault boundaries. This has been observed for some large ruptures using 'back-projection' techniques of teleseismic high-frequency body waves (e.g. [START_REF] Vallée | Ten Year Recurrence Time between Two Major Earthquakes Affecting the Same Fault Segment[END_REF]. Using the ray theory and the concept of isochrones (that is, the locus of energy emissions arriving at a station simultaneously), [START_REF] Spudich | Use of Ray Theory to Calculate High-Frequency Radiation from Earthquake Sources Having Spatially Variable Rupture Velocity and Stress Drop[END_REF] demonstrated that the PGA is proportional to the temporal changes of isochrone velocity, which depends on spatial variations of rupture velocity and slip velocity function. Deploying kinematic rupture simulations, [START_REF] Schmedes | Near-Source Ground Motion along Strike-Slip Faults: Insights into Magnitude Saturation of PGV and PGA[END_REF] showed that for a strike slip homogeneous rupture, the strongest changes of isochrones velocity is at a 'critical point', which remains at an almost constant position on the top fault boundary for all near fault stations. The spatial variations of near-fault PGA values are then essentially controlled by the station positions with respect to the critical point. Besides, [START_REF] Causse | Constraining the Roughness Degree of Slip Heterogeneity[END_REF] claimed that PGA is directly connected to the characteristic length of static slip asperities, small slip asperities generating larger PGA values. This contradicts [START_REF] Crempien | Within-Event and Between-Events Ground Motion Variability from Earthquake Rupture Scenarios[END_REF], who showed that large asperities increase ground motion coherency and lead to higher PGA values. Furthermore, [START_REF] Beresnev | Factors Controlling High-Frequency Radiation from Extended Ruptures[END_REF] claims that high-frequency ground motion is much more sensitive to the peak sliprate (that is, the impulse character of the slip-rate function) than slip heterogeneities. These results shed the light on the need for further investigations on the link between the heterogeneity level of kinematic rupture parameters and high-frequency ground motion. Finally, small-scale source heterogeneities pertain seismic motion wavelengths that are difficult to model owing to uncertainties in the propagation medium. Small-scale heterogeneities remain then poorly resolved by source studies.

Based on heterogeneous spontaneous dynamic rupture simulations, several authors showed that local kinematic source parameters may not be independent but correlated (e.g. [START_REF] Schmedes | A Kinematic Rupture Model Generator Incorporating Spatial Interdependency of Earthquake Source Parameters[END_REF][START_REF] Song | Pseudo-Dynamic Source Modelling with 1-Point and 2-Point Statistics of Earthquake Source Parameters[END_REF][START_REF] Song | Developing a Generalized Pseudo-Dynamic Source Model of Mw 6.5-7.0 to Simulate Strong Ground Motions[END_REF]. Such correlations may also impact the ground motion and its variability. Thus, larger correlation between source parameters may produce stronger PGV near the source [START_REF] Song | Developing a Generalized Pseudo-Dynamic Source Model of Mw 6.5-7.0 to Simulate Strong Ground Motions[END_REF]. Nevertheless, the level of correlation between kinematic source parameters remains poorly constrained. (Schmedes, Archuleta, and Lavallee 2010) analyzed a database of dynamic strike-slip rupture models computed using various initial distributions of stresses and strength properties and concluded that slip and rise time are strongly correlated. Nevertheless, the correlation pattern between the rupture velocity and the slip is unclear, and strongly depends on the parameters and their correlation in linear slip weakening friction laws, often assumed in the dynamic simulations [START_REF] Song | Developing a Generalized Pseudo-Dynamic Source Model of Mw 6.5-7.0 to Simulate Strong Ground Motions[END_REF].

This study aims to identify the kinematic source parameters that mostly control the PGA, the PGV and their variability, deploying numerical ground motion simulations of M 7 vertical strike-slip ruptures at a set of 5 near-fault stations as well as in farfield conditions. The choice of the simulation maximum frequency is delicate because a large maximum frequency results in unreasonably long computation time, while a too small frequency may not catch the real characteristic length of physical processes involved in the generation of high-frequency ground motion. Table 2-1 reports the frequencies that mostly contribute to the PGA and PGV values for the strong motion data selected in the PEER database (https://ngawest2.berkeley.edu/, last access June 27th, 2018) for 6.6 < M < 7.3, Vs30 > 600 m/s ('soft rock') and recorded less than 50 km from the rupture (Rrup < 50 km), based on Stockwell-transform analysis (Stockwell, Mansinha, and Lowe 1996) (see Figure 2-1). Among the 15 available recordings, 9 have their PGA controlled by frequencies smaller than 5 Hz, while the PGV is systematically dominated by frequencies below 5 Hz. Although the reported PGA and PGV values are controlled by a mix of source and wave propagation effects, this simple real data analysis provides a gross estimation of the frequency range that mostly contributes to the PGA and PGV. In the present study, the ground motion simulations are performed in the frequency range [0-5 Hz]. First, we identify the mechanisms of PGA generation in homogeneous kinematic ruptures. Hence, we study the PGA generated in heterogeneous kinematic ruptures. Afterwards, we run a sensibility analysis to determine the kinematic source parameters that mostly contribute to the PGA and PGV and their variability. We consider not only "large-scale" source parameters (average stress drop and average rupture velocity) but also "local-scale" parameters, that is statistical parameters controlling the level of source heterogeneity as well as the level of correlation between source parameters. This latter point is important because the distributions of "large-scale" source parameters are easier to constrain than the level of source heterogeneity in the framework of ground motion simulation of potential future earthquakes. Earthquake ruptures generate seismic waves that travel from the source to the surface and cause ground motions over a wide range of frequencies. One approach to describe the source process is the so-called kinematic approach, which consists in a priori prescribing the displacement discontinuity across the fault surface. The local slip-rate function needs to be specified (e.g. [START_REF] Liu | Prediction of Broadband Ground-Motion Time Histories; Hybrid Low/High-Frequency Method with Correlated Random Source Parameters[END_REF][START_REF] Tinti | A Kinematic Source-Time Function Compatible with Earthquake Dynamics[END_REF]) to describe the space-time evolution of slip along the fault by means of kinematic parameters. We use the pseudo-dynamic source model developed by [START_REF] Song | Pseudo-Dynamic Source Modelling with 1-Point and 2-Point Statistics of Earthquake Source Parameters[END_REF] for a rectangular fault plane. In this model, kinematic source parameters are calibrated using a suite of spontaneous heterogeneous dynamic rupture simulations. The rupture starts from the hypocenter and expands over the fault plane with a rupture speed 𝑉 𝑟 . Each point on the fault slips as it is reached by the rupture front and is characterized by the final slip value (𝐷) and the peak slip velocity (𝑝𝑠𝑣) or the slip duration, also called the rise time (𝑇 𝑟𝑖𝑠𝑒 ). In order to characterize the spatial variability of the kinematic source parameters (𝑉 𝑟 , 𝐷 and 𝑝𝑠𝑣) over the fault area, two statistical properties are considered. First, the 1-point statistics is defined for a given fault point by the mean value (𝜇) and the standard deviation (𝜎) of the considered source parameter, considering a normal distribution. Second, the 2-point statistics is defined by the correlation lengths (𝑎 𝑥 and 𝑎 𝑧 , representing the characteristic length of heterogeneities along-strike and along-dip, respectively) and the spatial cross-correlation, defined by the correlation coefficient (𝜌) between any pair of kinematic parameters at a given point, and by a correlation function. We use a Von Karman autocorrelation function [START_REF] Mai | A Spatial Random Field Model to Characterize Complexity in Earthquake Slip[END_REF].

Note that our statistical model is stationary, which implies that the statistics of any parameter is constant over the fault plane.

Earthquake source parameterization

We generate rupture models equivalent to a moment magnitude M = 7 (Figure 2-2). The rupture length L = 70 km and width W = 14 km are derived from the 𝑀 -𝐿 scaling relationship by [START_REF] Papazachos | Global Relations between Seismic Fault Parameters and Moment Magnitude of Earthquakes[END_REF]). Note that the relationship provided by [START_REF] Thingbaijam | New Empirical Earthquake Source-Scaling Laws[END_REF] result in close values of the rupture dimension (L = 67 km and W = 19 km). The mean value of the slip 𝜇 𝐷 is then defined by:

𝜇 𝐷 = 𝑀 0 𝐺 𝐿 𝑊
, where 𝑀 0 is the seismic moment and 𝐺 is the shear modulus. We make sure that the maximum slip does not exceed the ceiling defined by [START_REF] Mcgarr | Maximum Slip in Earthquake Fault Zones, Apparent Stress, and Stick-Slip Friction[END_REF] as a function of magnitude (500 cm in our case). For a strike slip rupture with L>>W, the average stress drop ∆𝜏 is expressed as:

∆𝜏 = 𝐶 * 𝐺 * 𝜇 𝐷 𝑊 (2-3)
where 𝐶 is a shape factor with a value close to 1 [START_REF] Kanamori | Theoretical Basis of Some Empirical Relations in Seismology[END_REF]. The average value of the rupture speed 𝜇 𝑉𝑟 is chosen in the range of values commonly reported by source studies [START_REF] Heaton | Evidence for and Implications of Self-Healing Pulses of Slip in Earthquake Rupture[END_REF]). We then assume that 𝜇 𝑉𝑟 = 0.85 𝑉 𝑠 , where 𝑉 𝑠 is the shear wave speed. The mean value of the 𝑝𝑠𝑣 is chosen from the database of spontaneous dynamic rupture simulations developed by [START_REF] Song | Pseudo-Dynamic Source Modelling with 1-Point and 2-Point Statistics of Earthquake Source Parameters[END_REF], 𝜇 𝑝𝑠𝑣 = 115 ± 25 𝑐𝑚/𝑠. The slip duration, also called the rise time 𝑇 𝑟𝑖𝑠𝑒 , is calculated with respect to the peak slip velocity 𝑝𝑠𝑣 and the slip value 𝐷, for a regularized Yoffe slip-rate function [START_REF] Tinti | A Kinematic Source-Time Function Compatible with Earthquake Dynamics[END_REF]. While the rise time is one of the most difficult parameters to extract from source models [START_REF] Heaton | Evidence for and Implications of Self-Healing Pulses of Slip in Earthquake Rupture[END_REF]), 𝑇 𝑟𝑖𝑠𝑒 is allowed to vary between 0.

PGA computation in the far-field approximation

The ground displacement for a homogeneous elastic medium in the far-field approximation 𝑢 𝐹𝐹 (𝑡) is proportional to the source time function, also called the moment rate function 𝑀 ̇(𝑡) . We first use Equation (2-4) to compute the displacement in the frequency domain:

𝑢 𝐹𝐹 (𝑋, 𝑓) = 1 4𝜋𝜌 0 𝑉 𝑠 3 1 X RP 𝑀 ̇(𝑓) 𝑒 - 2𝜋𝑖𝑓𝑋 𝑉 𝑠 𝑒 - 𝜋𝑓𝑋 𝑉 𝑠 𝑄 𝑠 (2-4)
where 𝑋 is the distance to the rupture, assumed equal to 100 km, 𝜌 0 is the rock density (𝜌 = 2.7 g/cm 3 ), 𝑉 𝑠 is the shear wave speed (𝑉 𝑠 = 3.58 km/s), 𝑅𝑃 is the average radiation pattern of the shear waves (according to [START_REF] Boore | Average Body-Wave Radiation Coefficients[END_REF], we assume RP = 0.63), 𝑄 𝑠 is the attenuation factor (we choose 𝑄 𝑠 = 220 (e.g. Heacock, Research, and Mines 1977). The attenuation of ground motion is represented by the geometrical attenuation 1/𝑋 and by the anelastic attenuation 𝑒𝑥𝑝 (-𝜋𝑓𝑋 𝑉 𝑠 𝑄 𝑠

). We then use inverse Fourier transform to obtain ground motion in the time domain and compute PGA as the maximum absolute value of the displacement second derivative. Note that a quarter-period-cosine taper is applied to the first second of the acceleration to remove the strong phase due to the sharp increase of rupture velocity at the rupture nucleation. Theoretical studies show that the rupture velocity increases smoothly during nucleation (e.g. [START_REF] Latour | Effective Friction Law for Small-Scale Fault Heterogeneity in 3D Dynamic Rupture[END_REF] and such initiation phases are not observed on real seismograms.

Mechanism of PGA generation for homogeneous ruptures

We start by investigating homogeneous ruptures, in order to identify the mechanisms of the far-field PGA generation in a simple rupture case. The slip, the rupture speed and the rise time are then constant along the rupture (𝜎 = 0), except at the fault boundaries due to the applied tapering. The parameters used for simulations are summarized in Table 2-2 (simulation A). We use the concept of isochrones to extract the part of the rupture that produces the PGA [START_REF] Spudich | Use of Ray Theory to Calculate High-Frequency Radiation from Earthquake Sources Having Spatially Variable Rupture Velocity and Stress Drop[END_REF]. Isochrones are all the points on the fault that radiate elastic waves such that the waves arrive at a given station at the same time. In the case of the farfield approximation (Equation (2-4), the isochrone at the PGA time is simply the rupture front at the PGA time (Figure 2-3 a, b, c and d). [START_REF] Spudich | Use of Ray Theory to Calculate High-Frequency Radiation from Earthquake Sources Having Spatially Variable Rupture Velocity and Stress Drop[END_REF] demonstrated that ground acceleration is proportional to the variations of isochrones velocity. In the far-field approximation, ground motion is then proportional to the variations of rupture velocity. Thus, for homogeneous ruptures, ground acceleration is essentially dominated by four peaks (Figure 2-3 e-1) corresponding to the times where the rupture reaches the four fault boundaries. For the chosen rupture nucleation position and fault boundary tapering function, the stopping phase generated by the rupture arrest at the fault top is responsible for the PGA. Since the tapering function determines the sharpness of the rupture stopping, it highly controls the PGA value. Thus, increasing the tapering length from 20% to 30% (a factor of 1.5), while holding the mean values of the parameters (D, Vr, Trise, psv) unchanged tends to decrease the PGA by a factor of 1.4 (Figure 2-3 e-2). In the following we explore the impact of various kinematic parameters on the PGA.

By decreasing the length of the rupture L while preserving the magnitude, we increase the slip D and therefore the stress drop ∆𝜏. Considering a decrease of L from 71 km to 55 km (that is by a factor of 1.29), D and ∆𝜏 increase by a factor of 1.29 (note that the mean value of Trise is unchanged, hence psv also increases by a factor of 1.29). According to Equation (2-2), PGA should then increase by a factor of ~1.22. Though the PGA changes by a slightly higher factor of 1.30, similar to the slip and the stress drop increase (Figure 2-3 e-3). Equation (2-2) is derived assuming a simple Brune's source model and random phases for the source spectrum [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes[END_REF]. PGA is then estimated using the random vibration theory and depends only on the corner frequency (that is, the overall rupture duration). Our source model also matches a Brune's source model (Figure 2-3 f). However PGA is not driven by the corner frequency but is controlled locally by the phase generated by the rupture stopping at the fault top, which is proportional to the local slip and hence to the stress drop.

We next analyze the effect of the rupture velocity Vr and the rise time Trise.

According to Equation (2-2), increasing Vr by a factor of 1.2 should increase the PGA by a factor of ~1.55 (Figure 2-3 e-4). This is slightly higher than the results of our simulations (~1.48). Again, we note that the simulated PGA cannot be modeled using random phases but is generated by a local process. Furthermore, reducing the Trise value from 4 to ~1.5 s only increases the PGA by a factor of 1.11 (Figure 2-3 e-5). This is because the PGA is mainly sensitive to Tacc (duration of the local slip acceleration phase), which equals 10 % of the rise time values. For Trise = 4 s, Tacc = 0.4 s. This implies that a decrease of Tacc only affects frequencies larger than 1/0.4 = 2.5 Hz, which is close to the maximum frequency of our simulations (5 Hz). Accelerogram obtained by increasing the psv to 277 cm/s resulting in a shorter rise time of 1.5 s. f) and g). Ground displacement spectra and the corresponding ground displacement time series for all the scenarios presented in e).

Sensitivity of Peak ground motions to source parameters

The above-mentioned tests illustrate some basic source mechanisms involved in the PGA generation. We now aim to quantify the PGA sensitivity to the source parameters, by perturbing one parameter at a time. We still compute the PGA in the far-field approximation (Equation (2-4)), but also for a network of near-field stations (Figure 2-2). We investigate the sensitivity to the 1-point and 2-point statistical parameters describing the rupture heterogeneity, which are very poorly constrained (e.g. [START_REF] Song | Developing a Generalized Pseudo-Dynamic Source Model of Mw 6.5-7.0 to Simulate Strong Ground Motions[END_REF], as well as the sensitivity to large-scale source parameters (average rupture velocity and average stress drop). We also include the PGV sensitivity in our analysis.

Computation of near-fault PGA and PGV

We synthesize near-fault ground motions in a 1D layered medium (Appendix 2-1) for stations located at rupture distances 𝑅 𝑟𝑢𝑝 of 5 km (station S1 and S2), 25 km (station S3 and S4) and 70 km (station S5) (Figure 2-2), using the representation theorem:

𝑈 𝑠𝑡𝑎𝑡𝑖𝑜𝑛 𝑖 (𝑡) ≈ ∑ 𝐺 𝑖 (𝑥, 𝑦, 𝑡) * 𝐹 𝑀 (𝑥, y, 𝑡) 𝑆 (2-5)
where * is the convolution operator. The summation over space integrates the contributions from the finite rectangular fault plane, discretized into a 2-D grid of subfaults. 𝐹 𝑀 (𝑥, 𝑦, 𝑡) is the slip rate function at position (𝑥, 𝑦) computed using the source model defined above, while 𝐺 𝑖 (𝑥, 𝑦, 𝑡) represents the Green's functions calculated using the discrete wavenumber technique in the frequency range 0-5 𝐻𝑧 [START_REF] Bouchon | The State of Stress on Some Faults of the San Andreas System as Inferred from Near-Field Strong Motion Data[END_REF][START_REF] Cotton | Dynamic Stress Variations Due to Shear Faults in a Plane-Layered Medium[END_REF].

Finally, the PGA and the PGV are computed using an orientation-independent measure proposed by [START_REF] Boore | Orientation-Independent Measures of Ground Motion[END_REF]) (GMRoTD50). This measure comprises a rotation of the two orthogonal components from 1 to 90, and evaluates the peak ground motion from the geometric mean of the rotated time series.

Appendix 2-2 illustrates a realization of kinematic rupture from case 5 (Table 2-2) with the calculated accelerations for both the EW and the NS components (black and red solid lines, respectively). The far-field acceleration computed according to Equation (2-4) is also shown.

Computation of the PGA and PGV sensitivity

We consider 8 rupture scenarios, reported in Table 2-2. In order to consider potential variations of the ground motion sensitivity due to the hypocenter position, we implement 3 different locations of the hypocenter (along-strike coordinates of 20 %, 50 % and 80 % of the rupture length). We generate 50 rupture realizations for each nucleation position, leading to 150 rupture realizations for each of the 8 scenarios.

The first scenario (referred to as case 1) assumes 𝜇 𝑉𝑟 = 3 𝑘𝑚/𝑠, 𝜎 𝐷 = 0.5 𝜇 𝐷 and 𝜎 𝑉𝑟 = 0.25 𝜇 𝑉𝑟 . The correlation lengths in the along-strike and down-slip directions are derived from the 𝑀 𝑤 -[𝑎 𝑥 ; 𝑎 𝑧 ] scaling relationship proposed by [START_REF] Mai | A Spatial Random Field Model to Characterize Complexity in Earthquake Slip[END_REF]. The final slip and the rupture velocity are supposed to be uncorrelated (D-Vr = 0) (Schmedes, Archuleta, and Lavallee 2010), while the final slip and the peak slip rate are assumed to be positively correlated (D-psv = 0.8 [START_REF] Song | Pseudo-Dynamic Source Modelling with 1-Point and 2-Point Statistics of Earthquake Source Parameters[END_REF]). Since such a positive correlation is found in most of the published studies (e.g. [START_REF] Schmedes | A Kinematic Rupture Model Generator Incorporating Spatial Interdependency of Earthquake Source Parameters[END_REF][START_REF] Song | Exploring Spatial Coherence between Earthquake Source Parameters[END_REF], the value of D-psv = 0.8 is kept for all the considered cases. The parameters of the reference case 1 are then perturbed independently to generate 7 additional cases of rupture scenarios. The sensitivity of the peak ground motion (PGM, which can refer to PGA or PGV) 𝑆 PGM 𝑘 to a given source parameter is computed as:

𝑆 PGM 𝑘 = (PGM k ̅̅̅̅̅̅̅̅ -PGM r ̅̅̅̅̅̅̅̅ ) PGM r ̅̅̅̅̅̅̅ ⁄ (𝑃 k -𝑃 r ) 𝑃 r ⁄ (2-6)
where PGM k ̅̅̅̅̅̅̅̅ and PGM r ̅̅̅̅̅̅̅ are the average PGM values computed for the case k and the reference case r, 𝑃 k and 𝑃 r are the values of the perturbed parameter for case k and r, respectively. 𝑆 PGM 𝑘 refers to the sensitivity of the peak ground motion to 𝜇 𝑉𝑟 (k=2, r=1), 𝜎 𝐷 (k=3, r=1), 𝜎 𝑉𝑟 (k=4, r=1), 𝑎 𝑥 (k=5, r=1), D-Vr + (k=6, r=4), D-Vr -(k=7, r=4) and ∆𝜏 (k=8, r=1).

In Figure 2

Results: sensitivity of PGA and PGV

We first consider the PGA sensitivity computed in the far-field approximation (Figure 23456). The first striking observation is that the PGA sensitivity is essentially controlled by the large-scale source parameter Vr and ∆τ and to a lesser extent to the amplitude of the rupture velocity fluctuation 𝜎 𝑉𝑟 . The sensitivity to Vr and ∆τ are similar as for an homogenous rupture. Note that the sensitivity to 𝜎 𝑉𝑟 can be negative. This is for instance the case if the PGA value is generated by the rupture arrest at the top fault (above the rupture nucleation) when the local rupture velocity is smaller than the median value. The PGA sensitivity to 𝜇 𝑉𝑟 has the largest average value (1.5) and is the least variable (16 th and 84 th percentiles equal to 0.9 and 1.9, respectively). This illustrates that a Vr increase always leads to a PGA increase whatever the mechanism of PGA generation is (by enhancing a rupture stopping phase or making the breaking of a slip asperity shorter). We also note that the average sensitivity is smaller than the value of 2.4 predicted by Equation (2-2). Furthermore, the average sensitivities to the statistical source parameters (𝜎 𝐷 , 𝜎 𝑉𝑟 , 𝑎 𝑥 , 𝜌 + , 𝜌 -) range between -0.17 and 0.14. In addition to these smaller values, the ratios between the average absolute and the 68 % confidence interval are also smaller, indicating that there is no clear tendency. The sensitivity patterns are similar for the PGA and PGV.

The results obtained at the near-field stations also reveal a strong dominance of the average rupture velocity effect, for both PGA and PGV values. The sensitivity values obtained at station S5, located at a Rrup distance equal to the rupture length, are very close to the value obtained in the far-field approximation, using Equation (2-4). As expected, the sensitivity to the rupture velocity is stronger at stations S2 and S4 due to forward directivity effects. Note that we do not compute the PGA sensitivity to the stress drop because it controls the fault dimension and hence modifies the position of the source heterogeneities with respect to the stations, making the sensitivity calculation meaningless.

Furthermore, it is interesting to note that the average PGA sensitivity to 𝜎 𝑉𝑟 (amplitude of the rupture velocity fluctuation) is significant, ranging from ~0.3 to ~0.5 at stations S3, S4 and S5, while it remains close to 0 at station S1, with the 16 th and 84 th percentiles of the order of -0.3 and 0.5, respectively. One explanation is that the PGA values at station S1 are controlled by a small portion of the rupture area located at the top fault center (close to the critical point defined by [START_REF] Schmedes | Near-Source Ground Motion along Strike-Slip Faults: Insights into Magnitude Saturation of PGV and PGA[END_REF] for homogeneous vertical strike slip ruptures), whatever the position of the nucleation. This area is statistically equally occupied by negative or positive fluctuations of rupture velocity around the median value and hence, the sensitivity of the peak values is centered around 0. The other stations, located at larger distances, "see" the whole rupture and the peak values are then controlled by the patches of maximum rupture velocity or rupture velocity drop. The same tendency is observed for PGV. We also note that the PGA and PGV sensitivities to ax and  have very small average values but are strongly variable, ranging between -1.1 and +1.7, and depend on the location of the station with respect to the rupture.

Finally, it is important to mention that the sensitivity alone is not sufficient to quantify the importance of a given source parameter because each parameter has its own range of variability. Defining the physically possible range of source parameters is important. For example, analyses of past earthquakes show that the average rupture speed can vary between 0.6 and 0.85 Vs [START_REF] Heaton | Evidence for and Implications of Self-Healing Pulses of Slip in Earthquake Rupture[END_REF], representing a potential ratio change of about 1.4. On the other hand, several studies report that the stress drop variability is well characterized by a log normal distribution with 𝜎 ln(∆𝛕) ≈ 0.8 (e.g. [START_REF] Mai | A Spatial Random Field Model to Characterize Complexity in Earthquake Slip[END_REF][START_REF] Courboulex | Stress-Drop Variability of Shallow Earthquakes Extracted from a Global Database of Source Time Functions[END_REF][START_REF] Cotton | What Is Sigma of the Stress Drop?[END_REF], which means a potential factor change of 4.3 considering the 68 % confidence interval.

Effect of the nucleation position

Our sensitivity study mixes rupture scenarios with various hypocenter positions.

Here we analyze the impact of the rupture nucleation position on the ground motion peak values. Figure 2-8(a) shows the coefficient of correlation between the PGA values for cases 'L' and 'C', 'L' and 'R', and 'R' and 'C' for all the ruptures cases and for all the stations. Figure 2345678shows the PGA values for the 50 realizations computed for the rupture scenario case2 at stations S1 and S2, considering hypocenter positions at distances of 0.2L ('L'), 0.5L ('C') and 0.8L ('R') from the fault left hand side corner. Station S2 shows a significant level of correlation for all the cases. The correlation is slightly smaller for station S1. This implies that the area which contributes to the PGA does not strongly depend on the hypocenter alongstrike position and that this area has a size smaller than the correlation length of rupture parameters (that is ax = 16 km and az = 5 km along strike and along dip respectively). This statement is in agreement with the study by [START_REF] Schmedes | Near-Source Ground Motion along Strike-Slip Faults: Insights into Magnitude Saturation of PGV and PGA[END_REF]) performed for a homogeneous rupture with a similar rupture size and mechanism. Based on their analysis, the locus of highest isochrones velocity variation (that is, the area generating the PGA) for station S1 would be located at the fault top edge 10-20 km away from the hypocenter surface projection for 'L' and 'R', and at fault top edge at mid-distance from the fault length for 'C', both corresponding to the same location on the fault top edge. For station S2, located beyond the rupture termination, the PGA is rather generated by the rupture arrest at the fault right hand corner. The correlation becomes small for the farthest stations S3, S4 and S5. The correlation at stations S3, S4 and S5 is especially for small values of the correlation length ax, meaning that the area controling the PGA generation for the 3 nucleation position is much larger than the characteristic size of source heterogeneities.

Peak Ground motions variability

Using our synthetic ground motion database, we next compute the between-event variability of PGA and PGV, referred to as 𝜎 ln(𝑃𝐺𝐴) and 𝜎 ln(𝑃𝐺𝑉) respectively. They are defined as the standard deviation of the natural logarithm PGA and PGV values, computed over 150 rupture realizations (50 realizations for each of the 3 nucleation positions). Appendix 2-3 (resp. Appendix 2-4) shows the values of the PGA (resp. PGV) variability obtained for each of the 8 simulation cases at stations S1, S3, S5 and in the far-field approximations with hypocenters on the right hand side and left hand side of the fault. Given the symmetry in these source-station configurations, the variability should converge to the same value for the two-hypocenter positions. We obtain discrepancies ranging from 1 % to 29 %, showing that 50 rupture realizations per hypocenter and per scenario provide a fairly good approximation of the ground motion variability, while preserving a reasonable number of rupture simulations. The values of 𝜎 ln(𝑃𝐺𝐴) and 𝜎 ln(𝑃𝐺𝑉) for each station and for the 8 analyzed cases are displayed in Figure 2 First, we notice that 𝜎 ln(𝑃𝐺𝐴) and 𝜎 ln(𝑃𝐺𝑉) are generally lower for small values of the correlation lengths ax and az (case 5: ax decreases from 16 to 4 km and az decreases from 5 to 2 km). This may arise because the characteristic size of heterogeneities gets smaller than the fault area controlling the peak values. Thus, the values of the rupture parameters averaged over this area are expected to be less variable over the different rupture realizations. Second, 𝜎 ln(𝑃𝐺𝐴) and 𝜎 ln(𝑃𝐺𝑉) are in overall higher for larger standard deviation of the rupture velocity (cases 4, 6 and 7). This is especially striking for stations S2 and S4. Indeed, the PGA values at stations S2 and S4 are controlled by the rupture stopping phase generated at the right hand side fault corner, which is very sensitive to the rupture velocity drop. Thus, a larger value of 𝜎 𝑉𝑟 gives rise to a larger variability of this stopping phase. We also observe larger values of 𝜎 ln(𝑃𝐺𝑉) for larger values of the standard deviation of slip (case 3). Finally, considering the aligned stations S1, S3 and S5, we observe a decrease of the PGA and PGV variability as the distance from the rupture increases. Again, this is because ground motion at distant stations is controlled by fault area larger than the characteristic length of rupture heterogeneities. The value of 𝜎 ln(𝑃𝐺𝐴) reaches 0.6 for the closest station. Note that the variability at station S5 is almost the same as the far-field approximation case.

As a matter of calibration of our chosen scenarios, we compare the variability issued from the source to the between-event variability τ reported in GMPEs. τ is a measure of the average deviation of the observed ground motion of any individual earthquake from the model median prediction, and therefore represents the variability in source parameters contributing to the ground motion but not modeled by GMPE such as ∆τ, Vr, source heterogeneities (L. Al [START_REF] Al Atik | The Variability of Ground-Motion Prediction Models and Its Components[END_REF]. As such, it represents an upper bound of the variability observed at one station for repeating ruptures on the same fault and for a given magnitude. The between-event variability reported in Ground Motion Predictions Equations (GMPEs) ranges between 0.23 and 0.42 for 𝜎 ln(𝑃𝐺𝐴) (see [START_REF] Causse | Are Stress Drop and Rupture Velocity of Earthquakes Independent? Insight from Observed Ground Motion Variability[END_REF], and between 0.25 and 0.4 for 𝜎 ln(𝑃𝐺𝑉) (L. Al [START_REF] Al Atik | The Variability of Ground-Motion Prediction Models and Its Components[END_REF]. These values constitute an upper bound for the ground motion variability due to source effects, and hence for our numerically computed 𝜎 ln(𝑃𝐺𝐴) and 𝜎 ln (𝑃𝐺𝑉) . Note that the between-event variability provided in GMPE studies is not valid in the near-fault area, because GMPE are essentially derived from far-field strong motion data. The simulated values are slightly above the values reported by GMPEs at the farthest station S5 and in the far-field approximation. The 𝜎 ln(𝑃𝐺𝐴) values reach ~0.6 at the near-fault station S1 and ~0.7 at the near-fault directive station S2, but these values may then be overestimated. As proposed by [START_REF] Causse | Are Stress Drop and Rupture Velocity of Earthquakes Independent? Insight from Observed Ground Motion Variability[END_REF] and (Chounet et al. 2018), anticorrelation between stress drop and rupture velocity may be a potential mechanism of PGA reduction, which would make the simulated PGA variability closer to observations.

Conclusion

Deploying ground motion simulations based on kinematic rupture models, we investigate the origin of the PGA and its variability. From the analyzed stations, we obtained the following results:

(1) the PGA and PGV values generated by heterogeneous ruptures are mainly controlled by the average rupture velocity and the average stress drop, and to a lesser extent by the standard deviation of rupture velocity (controlling the amplitude of rupture velocity fluctuations). The other statistical source parameters (correlation length and correlation between slip and rupture velocity) have a very weak impact on the average PGA values;

(2) in the far-field approximation the PGA values are mainly carried by the rupture stopping phase generated at the top fault edge above the rupture nucleation, like for simple homogeneous ruptures. Thus, the interaction between the rupture heterogeneities and the rupture top edge plays a key role in the PGA generation. This also implies that the choice of the function used to taper the distributions of final slip and rupture velocity strongly impact the simulated PGA values.

(3) the ground motion variability tends to increase with increasing amplitude of rupture velocity fluctuations, amplitude of slip fluctuations and characteristic size of source heterogeneities;

(4) 𝜎 ln(𝑃𝐺𝐴) and 𝜎 ln(𝑃𝐺𝑉) tend to increase as the distance to the rupture Rrup decreases. The obtained variability is slightly larger than the reported betweenevent variabilities at distances of the order of the rupture dimension and reaches 0.6 for Rrup = 5 km for PGA. The PGV variability is lower than the PGA variability.

Appendix

EARTHQUAKE

We have seen in chapter 2 that the rupture speed and the location of the site with respect to the hypocenter highly affect the peak ground motion in the near fault. When the rupture propagates from the hypocenter towards a site near the fault, the site is said to be located in the forward directivity region; when the rupture front propagates toward the site and at a velocity almost equal to the shear-wave velocity of the ground, all the seismic energy radiated from the fault rupture arrives at the site in a single, short-duration pulse. Structures located in the forward directivity region and having their natural periods close to the period of the directivity pulse are facing an amplification of their seismic response. Predicting the pulse period is therefore important, and it is commonly related to the magnitude of the earthquake. Even though pulse-like recordings are rare, however, there is a large variability of the pulse period values recorded for a given event. Can we do better to estimate the period of the pulse?

This paper is published in Bulletin of the Seismological Society of America (BSSA), 107, 1, pp. 308-318.

-

------------------------------------------------------------------------------------------------------

Introduction

The directivity of the earthquake rupture propagation gives rise to a large variability of the ground motions recorded at a given distance from the source over various source-receiver azimuths (e.g. [START_REF] Somerville | Modification of Empirical Strong Ground Motion Attenuation Relations to Include the Amplitude and Duration Effects of Rupture Directivity[END_REF][START_REF] Spudich | Directivity in NGA Earthquake Ground Motions: Analysis Using Isochrone Theory[END_REF]. In particular, the energy of the seismic waves successively released from the fault constructively interferes in the forward direction of the rupture, which makes the amplitude of the ground shaking large, especially when the rupture speed approaches the shear wave speed. This results in a large amplitude S-wave, called the pulse. Such pulses are essentially observed in case of forward directivity and in the near-fault region, the probability that ground-velocity is pulse-like being maximum at the vicinity of the rupture termination [START_REF] Shahi | An Empirically Calibrated Framework for Including the Effects of Near-Fault Directivity in Probabilistic Seismic Hazard Analysis[END_REF]; see also Figure 3-1). The pulse is of particular interest from a structural earthquakeengineering point of view, because the demand on the structure is amplified when the natural period of the structure equals the pulse period (e.g. [START_REF] Biggs | Introduction to Structural Dynamics[END_REF](Veletsos, Newmark, and Chelapati 1965) [START_REF] Anderson | Uncertainties in Establishing Design Earthquakes[END_REF][START_REF] Hall | Near-Source Ground Motion and Its Effects on Flexible Buildings[END_REF]. In particular, the pulse period has been shown to be a critical parameter for design spectra, strength reduction factors, damping modification factors, residual displacements and ductility demands [START_REF] Alavi | Effects of Near-Fault Ground Motions on Frame Structures[END_REF][START_REF] Hubbard | Damping Coefficients for Near-Fault Ground Motion Response Spectra[END_REF][START_REF] Liossatou | Near-Fault Effects on Residual Displacements of RC Structures[END_REF][START_REF] Mavroeidis | Near-Fault Ground Motions, and the Response of Elastic and Inelastic Single-Degree-of-Freedom (SDOF) Systems[END_REF][START_REF] Ruiz-Garcia | Inelastic Displacement Ratios for Seismic Assessment of Structures Subjected to Forward-Directivity Near-Fault Ground Motions[END_REF]. Accurate predictions of directivity pulse periods are then crucial for near fault seismic risk assessment. Several studies showed that the pulse period scales with moment magnitude, and proposed empirical relationships to relate the pulse period to earthquake magnitude (e.g. [START_REF] Somerville | Development of an Improved Representation of near Fault Ground Motions[END_REF]2003;[START_REF] Shahi | An Empirically Calibrated Framework for Including the Effects of Near-Fault Directivity in Probabilistic Seismic Hazard Analysis[END_REF][START_REF] Bray | Characterization of Forward-Directivity Ground Motions in the near-Fault Region[END_REF][START_REF] Mavroeidis | A Mathematical Representation of Near-Fault Ground Motions[END_REF]. Recently, [START_REF] Cork | Effects of Tectonic Regime and Soil Conditions on the Pulse Period of Near-Fault Ground Motions[END_REF] claimed that the pulse period may be related to other source features, like the tectonic regime or the stress drop.

In this article, we investigate the spatial dependency of the pulse period and its dependency on some basic parameters describing the rupture process. We propose that the pulse period can be approximated from a simple equation including the relative location of the observation points with respect to the rupture (𝐷, 𝑐𝑙𝑠𝐷), the ratio between the rupture speed and the shear wave speed (𝑉 𝑟 / 𝑉 𝑠 ) and the duration of the local slip (𝑇 𝑟𝑖𝑠𝑒 ). Our simple model is first validated on a synthetic dataset based on simulations of extended strike-slip ruptures. The model is then compared to real data selected from the NGA-West2-database, compiling 110 observations of pulses periods from 10 strike-slip events and 6 non-strike slip events.

3.3 Relationship between pulse period, rupture parameters, and station position based on analysis of synthetic velocity time series

Simulation of velocity time series

First, we simulate a suite of velocity time series for vertical strike slip ruptures in a homogeneous medium at a set of 12 stations (Figure 3-1). The stations are located in the near-fault region (distance from the surface fault projection smaller than half the rupture length). The rupture initiates at one edge of the fault and propagates at a constant rupture velocity, so that each station sees part or totality of the rupture arriving. For comparison, Figure 3-1 displays the contour lines of the probability that the observed velocity is pulse-like, as computed by [START_REF] Shahi | An Empirically Calibrated Framework for Including the Effects of Near-Fault Directivity in Probabilistic Seismic Hazard Analysis[END_REF] empirical equation, derived from a real dataset of strike-slip velocity time series. The Green's functions are computed using the discrete wavenumber technique [START_REF] Bouchon | The State of Stress on Some Faults of the San Andreas System as Inferred from Near-Field Strong Motion Data[END_REF], up to a frequency of 3 𝐻𝑧, considering a shear wave velocity 𝑉 𝑠 = 3 𝑘𝑚/𝑠. The rupture area is 40 * 14 𝑘𝑚 2 and the average slip is 1 𝑚, so that the simulated ruptures correspond to 𝑀~6. The rise time (i.e. the local slip duration) and the final slip are assumed to be constant over the fault plane. In order to analyze the effect of slip heterogeneity, we also consider heterogeneous slip distributions (Figure 3-2a). Those distributions are generated assuming a 𝑘 -2 slope in the wavenumber domain beyond a corner wavenumber 𝑘 𝑐 (e.g. [START_REF] Causse | New Approach for Coupling K-2 and Empirical Green's Functions: Application to the Blind Prediction of Broad-Band Ground Motion in the Grenoble[END_REF]. The inverse of 𝑘 𝑐 is proportional to the slip correlation length (characteristic size of slip heterogeneity), while the level of the high frequency slip spectrum is related to the standard deviation of slip (amplitude of slip heterogeneity). The value of the pulse period 𝑇 𝑝 is next computed using the algorithm developed by [START_REF] Baker | Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis[END_REF], based on continuous wavelet transform. The pulse period is equal to the period of the wavelet associated with the largest coefficient, considering Daubechies wavelet of order 4 as a mother wavelet. Note that the [START_REF] Baker | Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis[END_REF] algorithm is used here for the sake of consistency, since it has been used to complete the NGA-West2 database analyzed in the subsequent section.

Figure 3-2b displays the fault-normal components of the simulated velocity at station 4 (left) and the corresponding extracted pulses (right) for constant or heterogeneous slip distributions. The result indicates that the pulse period is almost not sensitive to the level of heterogeneity. Figure 3-2b also shows that considering the whole rupture area or only the area defined by the shaded rectangle (Figure 3-2a) does not modify the velocity pulse. The only noticeable difference between the two synthetics is a more pronounced stopping phase when the rupture is shortened and stops aside the station. This clearly illustrates that the pulse essentially arises from a coherent summation of waves emitted by the fault area that ruptures towards the site. This is also pointed out by [START_REF] Mavroeidis | Effect of Fault Rupture Characteristics on Near-Fault Strong Ground Motions[END_REF], who used the concept of isochrones to relate the characteristics of near-fault directivity pulses to the rupture properties. They showed that velocity pulses observed at stations located close to the fault surface projection are associated with a fault area of large isochrones velocities, which extends from the hypocenter to the top of the fault.

Simple relation between pulse period, rupture parameters and station position

Based on the previous analysis, we propose that the pulse period can be approximated by the following simple equation:

𝑇 𝑝𝑢𝑙𝑠𝑒 = ( 𝐷 𝑉 𝑟 + 𝑐𝑙𝑠𝐷 𝑉 𝑠 - ℎ𝑦𝑝𝐷 𝑉 𝑠 ) + 𝑇 𝑟𝑖𝑠𝑒 , (3-1)
in which 𝐷 is the length of the fault area that ruptures toward the site, measured between the hypocenter and the closest point from the fault to the site, 𝑐𝑙𝑠𝐷 is the closest distance from the recorded site to the ruptured fault area, and ℎ𝑦𝑝𝐷 is the hypocentral distance (as shown in Figure 3-2a). Those parameters are illustrated in a more general case (for any fault mechanism) in Figure 3-3. Note that Equation (2-1) is valid for subshear ruptures only, that is when the ratio between the rupture speed and the shear wave velocity 𝑉 𝑟 /𝑉 𝑠 is smaller than 1. To test the robustness of this simple approximation, we compare the pulse period values extracted from our simulations using the wavelet algorithm [START_REF] Baker | Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis[END_REF]) and the values obtained from Equation (2-1), for different values of 𝑉 𝑟 /𝑉 𝑠 and 𝑇 𝑟𝑖𝑠𝑒 at the whole set of stations (Figure 34). The standard deviation of the natural logarithm residuals equals 0.27 , and the coefficient of correlation between observations and predictions equals 0.91, indicating that Equation (2-1) provides a fairly good approximation of the pulse period. However, because the shape of the considered wavelet does not systematically match the shape of the synthetic velocity pulses, Equation (2-1) sometimes underestimates the pulse periods obtained from the wavelet algorithm. This is the case for stations located next or beyond the rupture termination (stations 6, 7, 9, 10 and 12) and when the 𝑉 𝑟 /𝑉 𝑠 ratio is lower than ~0.8. This is illustrated in Figure 5, showing that the duration of the synthetic velocity pulse at station 12, well delimited by the first S-wave arrival and a stopping phase, is well predicted by Equation (2-1). Nevertheless, the wavelet used to approximate the pulse has a significantly larger period. This raises the question of the real meaning of the pulse period, its relation to the pulse duration and how to properly measure it. Recently, [START_REF] Cork | Effects of Tectonic Regime and Soil Conditions on the Pulse Period of Near-Fault Ground Motions[END_REF] pointed out that the use of various techniques can results in significantly different values of the pulse period. This issue remains, however, beyond the scope of this study.

Comparison between predicted pulse period (Equation

(2-1)) and real observations (NGA-West2 Database)

Data selection

To test Equation (2-1) with real data, we refer to the NGA-West2 database. Earthquake data are selected based on the availability of the velocity pulse period and the rupture parameter values. Two additional earthquakes, not fully described in the database, are however included in our dataset: the 2003 Bam, Iran earthquake and the 2004 Parkfield, California earthquake. For the Bam earthquake, rupture parameters are determined from the study by [START_REF] Bouchon | Some Insight on Why Bam (Iran) Was Destroyed by an Earthquake of Relatively Moderate Size[END_REF]. For the Parkfield event, we refer to [START_REF] Twardzik | Robust Features of the Source Process for the 2004 Parkfield, California, Earthquake from Strong-Motion Seismograms[END_REF] for the 𝑉𝑟 and 𝑉𝑠 values. Furthermore, according to (Custodio et al., 2005), we choose 𝑇 𝑟𝑖𝑠𝑒 = 0.88 𝑠. The source parameters of the considered events are reported in Table 3-1, and information about each station (velocity pulse period and station position) can be found in Appendix 3-1.

Results

Figure 3-6 displays the values of the extracted (Tp(NGA)) and the calculated (Tp(eq1)) velocity pulse periods at all the stations. It illustrates the large variability of the recorded pulse period for a given earthquake. In order to quantify the misfit between the observations and the predictions, we compute the standard deviation of the natural logarithm residuals. We obtain 𝜎 ln 𝑇 = 0.58 using the whole dataset. Furthermore, we note that the distribution of the natural logarithm residuals is not centered around 0 but around 0.2, which means that Equation (2-1) slightly overestimates the observed Tp(NGA) values (~+ 20%). One possible explanation is that the pulse period may be controlled by asperity dimensions, which are typically less than the distance D used in Equation (2-1). This potential behavior is not captured by our k -2 heterogeneous rupture simulations, which show that the pulse period is almost insensitive to the level of slip heterogeneity, but may be revealed by more complex (non-stationary) slip distributions. Finally, the correlation coefficient between the predicted and the observed pulse periods is 𝑐𝑜𝑟𝑟 = 0.82, indicating that even if significant discrepancies can be observed for certain events, Equation (2-1) explains most of the spatial variability of the pulse period.

The largest discrepancies are observed at some of the stations that recorded the 1994 Northridge earthquake, for which Equation (2-1) overestimates the pulse period by a factor larger than 4. These large discrepancies may be explained by the simplicity of Equation (2-1), in which the source process is modeled by means of a single fault plane. Using teleseismic waveform inversion and analysis of aftershock distribution, [START_REF] Thio | Source Complexity of the 1994 Northridge Earthquake and Its Relation to Aftershock Mechanisms[END_REF] showed that the 1994 Northridge rupture was however complex, consisting of 3 sub-events that ruptured with a time-shift about 2𝑠. We also note that Equation (2-1) overestimates the pulse period for the 1995 Kobe earthquake (by a value of about 2𝑠). This arises because the average rupture speed of 2.7 𝑘𝑚/𝑠 provided in the NGA-West2 database may be underestimated. Table 3-2 lists the rupture speed values as reported by several published finite-source inversion models (available at SRCMOD database, [START_REF] Mai | SRCMOD: An Online Database of Finite-fault Rupture Models[END_REF], indicating the rupture speed is probably closer to 3 𝑘𝑚/𝑠. Using 𝑉 𝑟 = 0.9 𝑉 𝑠 = 3.1 𝑘𝑚/𝑠 (e.g. Yujia Guo et al., 2013) results in a better fit with observations, as indicated by triangles in Figure 3456. Finally, the pulse period observed at station BOL during the 1999 Duzce earthquake (Figure 3-6, filled circle) is strongly overestimated. This may arise because the 𝑉 𝑟 /𝑉 𝑠 value reported in the NGA-West2 database is 0.8, while the rupture towards BOL station propagated at a supershear speed (that is 𝑉 𝑟 /𝑉 𝑠 >1) [START_REF] Bouchon | Faulting Characteristics of Supershear Earthquakes[END_REF]). As reported in several studies (e.g. [START_REF] Bernard | Shear Mach Wave Characterization for Kinematic Fault Rupture Models with Constant Supershear Rupture Velocity[END_REF], supershear ruptures are associated with a shock wave propagating in the near source region, called Mach front. This shock-wave is characterized by a pulse of large amplitude and a short duration. Note that supershear rupture was also observed for the 1999 Koaceli earthquake, but on a fault segment that ruptured beyond the 4 stations considered in our study (Arcelik, Gebcze, Izmit and Yarimca) [START_REF] Bouchon | Faulting Characteristics of Supershear Earthquakes[END_REF]). Thus, these 4 stations were not affected by the Mach front.

After screening out the data from the Northridge and the Duzce earthquakes, and considering 𝑉 𝑟 / 𝑉 𝑠 = 0.9 for the Kobe earthquake, the value of 𝜎 ln 𝑇 drops from 0.58 to 0.47 considering the whole data set, from 0.63 to 0.55 considering strike-slip earthquakes only, and from 0.52 to 0.38 for non-strike slip earthquakes. The correlation coefficient remains however almost unchanged (𝑐𝑜𝑟𝑟 = 0.84). Separating the events according to their mechanism, we obtain 𝑐𝑜𝑟𝑟 𝑆𝑆 = 0.86 for strike-slip events and 𝑐𝑜𝑟𝑟 𝑁-𝑆𝑆 = 0.68 for non-strike-slip events. This seems to indicate that while the 𝜎 ln 𝑇 is smaller for non-strike-slip events, Equation (2-1) still better predicts the spatial variability of the pulse period for strikeslip events.

Discussion

Several studies proposed empirical models relating the pulse period to earthquake magnitude, assuming a linear relationship between log(𝑇 𝑝 ) and 𝑀 (e.g. [START_REF] Cork | Effects of Tectonic Regime and Soil Conditions on the Pulse Period of Near-Fault Ground Motions[END_REF][START_REF] Shahi | An Empirically Calibrated Framework for Including the Effects of Near-Fault Directivity in Probabilistic Seismic Hazard Analysis[END_REF][START_REF] Somerville | Magnitude Scaling of the near Fault Rupture Directivity Pulse[END_REF]. All these models show evidence of an increase of 𝑇 𝑝 with 𝑀. Note that our model (Equation (2-1)) does not include magnitude or seismic moment explicitly (or any parameter that scales with seismic moment like average slip or rupture length). The increase of 𝑇 𝑝 with 𝑀 is however implicitly expressed through the parameter 𝐷 (length of the fault area rupturing toward a given station). This is simply because as magnitude increases, the length of the rupture increases, and therefore the range of potential values of 𝐷, and hence 𝑇 𝑝 ), also increases (as shown in Figure 3-7a). Figure 3-7b clearly indicates that the 𝑇 𝑝 values follow a similar trend, and that large events can also have small 𝑇 𝑝 values. Thus, our model provides a simple physical basis to explain how 𝑇 𝑝 is linked to magnitude. Although 𝑇 𝑝 -𝑀 empirical models give a practical and direct way to predict 𝑇 𝑝 for a potential scenario earthquake, our approach requires the knowledge of some rupture parameters (rupture velocity and rise time). Though these parameters only describe the basic features of the rupture propagation, they remain difficult to predict a priori. After [START_REF] Heaton | Evidence for and Implications of Self-Healing Pulses of Slip in Earthquake Rupture[END_REF] and [START_REF] Somerville | Characterizing Crustal Earthquake Slip Models for the Prediction of Strong Ground Motion[END_REF], the values of the 𝑉𝑟/𝑉𝑠 ratio observed for most ruptures are in the range [0.6 -0.9]. The rupture velocity can also vary locally within a single rupture (e.g. [START_REF] Archuleta | A Faulting Model for the 1979 Imperial Valley Earthquake[END_REF]. In order to quantify the variability of the velocity pulse period that would be predicted from Equation (2-1) for a future earthquake, we compute 𝑇 𝑝 for 10 000 rupture scenarios for a 𝑀6 strike-slip rupture (Figure 3-2a). We assume that the 𝑉 𝑟 /𝑉 𝑠 ratio is uniformly distributed in the range [0.6 -0.9]. We also assume that the rise time values are uniformly distributed in the range [0.1 -1]. Figure 3-8a and Figure 3-8b represent the spatial distribution of the mean pulse period and logarithm standard deviation, respectively. Standard deviation varies from 0.09 to 0.56, depending on the location of the station with respect to the fault. In overall, stations located close to the fault surface projection or beyond the fault termination are more sensitive to uncertainties in the source parameters. The highest uncertainty is observed at some specific locations at the vicinity of the hypocenter, for which the pulse period equals the smallest possible value, i.e. the rise time value. The pulse period uncertainty is then directly controlled by the rise-time uncertainty. To reduce uncertainty in the pulse period predictions for a potential future event, it is then essential to constrain the range of physically realistic rupture parameter values. Some studies established links between the rupture speed and some fault properties or other physical parameters describing the rupture, which may help refining the a-priori estimation of the rupture speed. [START_REF] Bouchon | Faulting Characteristics of Supershear Earthquakes[END_REF] observed the fault rupture surface of several earthquakes with supershear ruptures and concluded that the rupture may propagate at a supershear speed only when the geometry of the fault is simple. Besides, [START_REF] Manighetti | Earthquake Scaling, Fault Segmentation, and Structural Maturity[END_REF] and [START_REF] Radiguet | Dependency of Near-Field Ground Motions on the Structural Maturity of the Ruptured Faults[END_REF] analyzed stress drop, another important source parameter, with respect to the so-called 'maturity' of faults. Maturity includes fault features like age, length and cumulative displacement on the fault. The authors conclude that mature faults are associated with low stress drop. Finally, [START_REF] Causse | Are Stress Drop and Rupture Velocity of Earthquakes Independent? Insight from Observed Ground Motion Variability[END_REF] propose that average stress drop and average rupture velocity may not be independent but anticorrelated, and propose joint distributions of these parameters. A-priori knowledge of the stress drop may then also help constraining the rupture velocity for a-priori estimations of the pulse periods.

Conclusions

This article presents a simple equation to predict the period of the velocity pulses that can be observed in the near-fault region and in the forward rupture direction. This equation is based on a few basic parameters: the location of the station with respect to the rupture, the velocity of the rupture propagation, the rise time and the shear wave velocity of the medium around the fault. Our approach is first validated by analyzing a suite of synthetic velocity time series of strike slip extended ruptures. The velocity pulse periods are computed from the [START_REF] Baker | Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis[END_REF] algorithm, based on wavelet transform. This analyses shows that (1) the pulse period is sensitive to the rupture length toward the station rather than the whole length of the fault; (2) the pulse period is not sensitive to the heterogeneity of the slip distribution on the fault plane; (3) the value of the pulse period as computed from the wavelet analysis can differ from the real duration of the directivity pulse at stations located next or beyond the rupture termination.

Our equation is then tested on a dataset build from NGA-West2 database, consisting of 110 observations of pulses periods from 10 strike-slip events and 6 nonstrike slip events. The standard deviation of the natural logarithm residuals between observations and prediction is ~0.5, and the correlation coefficient between observations and predictions is 0.84 (0.86 considering strike-slip events only, and 0.68 for non-strike slip events). This indicates that despite significant discrepancies are observed at some stations, our simple model fairly well reproduces the spatial variability of the pulse periods recorded during an earthquake, especially for strike slip events. As mentioned above, some of these discrepancies can be explained by unsuitable values of the rupture velocity, complexity of the fault geometry or inadequacy between the pulse duration and the value of the pulse period, which depends on the technique used to extract the pulse. Note that some of these discrepancies may also be explained by the variability of the site conditions, which is not considered in Equation (2-1). In overall, soil site are generally characterized by larger values of the pulse period than rock sites (e.g. [START_REF] Bray | Characterization of Forward-Directivity Ground Motions in the near-Fault Region[END_REF][START_REF] Cork | Effects of Tectonic Regime and Soil Conditions on the Pulse Period of Near-Fault Ground Motions[END_REF]P. G. Somerville, 2003).

Finally, it is important to mention that proper a priori estimations of the pulse period for a potential future earthquake rely on a proper a priori knowledge of the location of the hypocenter, the rupture velocity and the rise time. The fault system accommodates the left-lateral motion between the Sinai-Nubai (African) plate and the Arabian plate (Figure 4-1) with a total of 105 km of left lateral transform motion between the plates since early Miocene (20 Myr).

Appendix

The DSTF, illustrated in Figure 4-2, is divided into three main segments: the ~520km long N-S trending southern section, linking the Red Sea rift to the Hula basin; the ~170km long NE-SW trending section in Lebanon; and the ~340km long N-S trending northern section from the Ghab basin in Syria into southern Turkey. Indeed, both Sinai and Arabian plates drift north at a rate of 18 mm/year for Arabia and at about 6 mm/year for Africa-Nubia. GPS campaigns from 1991 to 2011 provide an average relative plate motion of ~5 mm/year, ~5 mm/year, and ~2 mm/year for the south, central, and north segments, respectively [START_REF] Reilinger | GPS Constraints on Continental Deformation in the Africa-Arabia-Eurasia Continental Collision Zone and Implications for the Dynamics of Plate Interactions[END_REF].

The DSTF poses a seismic threat to the population centered in its vicinity, affecting present day Lebanon, Cyprus, Syria, Jordan, Israel and Palestine. The historical and recent earthquake catalogue that is used for evaluating the seismic hazard of the region was compiled from different sources: the region has long been inhabited and it provides a rare historical archive including numerous earthquake witnesses. Historical earthquake information was compiled from different Arabic, Islamic, Jewish and Christian historians who assembled descriptions of earthquakes mentioned in ancient literature [START_REF] Ambraseys | Value of Historical Records of Earthquakes[END_REF]Ben-Menahem 1991a;[START_REF] Shapira | A New Estimate for the Epicenter of the Jericho Earthquake of 11 July 1927[END_REF]. Large earthquakes are known to have occurred along this fault zone with recurrent magnitudes 6-7.5 earthquakes over the historical period.

At the end of the 19 th century, instrumental seismic monitoring in the Middle East started with the installation of station Helwan HLW in Cairo, Egypt. From 1898 until 1912, HLW was the only seismological station that operated in the region, the closest stations being remote and including stations at Athens and Istanbul. Another important station that started operating in 1912 was Ksara, which is located in central-east part of Lebanon [START_REF] Allen | Earthquake Early Warning for Israel: Recommended Implementation Strategy[END_REF]. Nowadays, earthquakes are monitored and analyzed through a network of seismological stations all over the region. Seismological stations located all over the world could also detect large magnitude earthquakes.

All together, the seismicity of the DSTF is a unique example of an area for which information and documentation on historical earthquakes cover a time span of more than four millennia. Nevertheless, among all continental transform systems, the DSTF is an exception because of its apparent last-century seismic quiescence and, therefore, variability of earthquake activity and faulting behavior. The instrumental seismicity (younger than 1900 A.D.) of the continental DSTF is of a relatively low level compared to other faults, such as the nearby East Anatolian fault or North Anatolian fault. Moreover, the apparent quiescence and the lack of major seismic events with M>6.0 on most fault segments in the last centuries are in contradiction with the historical catalogue and related report of faulting events over the last 3,000 years or so along the continental DSTF [START_REF] Guidoboni | Catalogue of Ancient Earthquakes in the Mediterranean Area up to the 10th Century[END_REF][START_REF] Ambraseys | Faulting Associated with Historical and Recent Earthquakes in the Eastern Mediterranean Region[END_REF][START_REF] Sbeinati | The Historical Earthquakes of Syria: An Analysis of Large and Moderate Earthquakes from 1365 BC to 1900 AD[END_REF]. [START_REF] Elnashai | Earthquake Hazard in Lebanon[END_REF] found that the frequency-magnitude distribution indicates that the recorded earthquake activity on the DSTF during the instrumental period (i.e. the past century) has been lower than expected to accommodate the Africa-Arabia relative plate motion. (Ambraseys and Barazangi 1989) analyzed the macroseismic data on ~350km long segment of the northern part of the DSFT for the period 1100-1988. The 10 major historical earthquakes (M>6.5) struck in three relatively short periods [1157-1202; 1404-1407; 1759-1796] with repeat times of 200-350 years. The lack of such large events during the past 100 years should not be interpreted to minimize potential earthquake hazard in this region.

A recent study published in Nature analyzed the past earthquakes on the DSFT south of Lebanon over the last 1600 years. [START_REF] Lefevre | Slip Deficit and Temporal Clustering along the Dead Sea Fault from Paleoseismological Investigations[END_REF] suggested that temporal clustering of earthquakes is a common behavior over the entire region. During each earthquake cluster, the entire fault eventually ruptures, although the spatial and temporal distributions of earthquakes seem to be random. [START_REF] Lefevre | Slip Deficit and Temporal Clustering along the Dead Sea Fault from Paleoseismological Investigations[END_REF]) also revealed that more than 2 m of slip are accumulated during the last 1600 years all along the DSFT, which would correspond to a magnitude Mw 7.3 event on each section to fully release the current accumulated slip. Nevertheless, in 1995, a large earthquake of Mw 7.3 occurred in the Gulf of Aqaba (Figure 4-2a), identified as the large earthquake along the entire DSFT for more than 200 years, raising the doubtful possibility of a new earthquake series clustering in the upcoming decades.

Overview of the seismicity within Lebanon

Five earthquakes stand out in the history of seismic activities in Lebanon: the earthquakes of 551, 1202, 1759 (two events) and 1837. Within the Lebanese bend, the DSTF fault splits into four main branches, with clear evidence for slip partitioning: the Yammouneh YF, the Roum RF and Rachaya RcF-Serghaya SF are left-lateral strike-slip faults, and the previously unrecognized offshore fault, the Mount-Lebanon Thrust MLT. Of which, YF constitutes the main fault branch of the DSTF within the Lebanese restraining bend (Figure 4-2b).

The Yammouneh fault (YF), which bisects the length of Lebanon and bends eastward to compress and uplift the Mount-Lebanon Range, has been linked to the historical 1202 earthquake of magnitude Ms 7.6, and is slipping at an estimated rate of 4-5 mm/year [START_REF] Plassard | Seismicité Du Liban[END_REF][START_REF] Ambraseys | Faulting Associated with Historical and Recent Earthquakes in the Eastern Mediterranean Region[END_REF]Ben-Menahem 1991b;[START_REF] Ellenblum | Crusader Castle Torn Apart by Earthquake at Dawn[END_REF][START_REF] Daëron | Constraints on the Post∼ 25-Ka Slip Rate of the Yammoûneh Fault (Lebanon) Using in Situ Cosmogenic 36Cl Dating of Offset Limestone-Clast Fans[END_REF][START_REF] Daeron | 12,000-Year-Long Record of 10 to 13 Paleoearthquakes on the Yammouneh Fault, Levant Fault System, Lebanon[END_REF][START_REF] Gomez | Global Positioning System Measurements of Strain Accumulation and Slip Transfer through the Restraining Bend along the Dead Sea Fault System in Lebanon[END_REF]. According to recent GPS data acquired in Lebanon (1999Lebanon ( , 2002Lebanon ( , 2010)), small lateral fault slip rates (2-4 mm/year) are detected with a slight slip rate decrease from south to north [START_REF] Vergnolle | GPS Measurements of Crustal Deformation in Lebanon: Implication for Current Kinematics of the Sinaï Plate[END_REF].

The Roum fault (RF) branches from the DSTF in south Lebanon, runs along the south-western boundary of the Mount-Lebanon range and is around 35 km in length. The RF produced the 1837 (MS 7.1) earthquake [START_REF] Nemer | Coseismic Growth of Sedimentary Basins along the Yammouneh Strike-Slip Fault (Lebanon)[END_REF][START_REF] Sigbjörnsson | Uncertainty Analysis of Strong-Motion and Seismic Hazard[END_REF].

The coupled Rachaiya (RcF) -Serghaya (SF) fault system is a left-lateral strainpartitioning complex through the Anti-Lebanon mountain Range: The RcF is approximately 45 km long and traces along the western flank of the Anti-Lebanon range. The SF is 100-150 km in length and traces from the south-eastern to the north western flank of the Anti-Lebanon range. The RcF-SF system shows 1.4 ± 0.1 mm/year of movement [START_REF] Gomez | Holocene Faulting and Earthquake Recurrence along the Serghaya Branch of the Dead Sea Fault System in Syria and Lebanon[END_REF] and produced estimated M6.7 and M7.2 earthquakes in 1759 that destroyed much of the region [START_REF] Elnashai | Earthquake Hazard in Lebanon[END_REF].

The recently identified offshore Mount Lebanon Thrust (MLT) is undergoing compressional displacements of 1.0-2.0 mm/year. The MLT trace at the surface lies mainly offshore, between Tripoli in the north and Saida in the south, cutting the seabed at not more than 8km from the coast of central Lebanon between Beirut and Anfeh [START_REF] Elias | Active Thrusting Offshore Mount Lebanon; Source of the Tsunamigenic A.D. 551 Beirut-Tripoli Earthquake[END_REF]). The most significant historical earthquake along the Mount-Lebanon Thrust was the tsunamigenic ~M7.5 earthquake in 551 [START_REF] Plassard | Seismicité Du Liban[END_REF][START_REF] Ambraseys | Faulting Associated with Historical and Recent Earthquakes in the Eastern Mediterranean Region[END_REF][START_REF] Guidoboni | Catalogue of Ancient Earthquakes in the Mediterranean Area up to the 10th Century[END_REF][START_REF] Darawcheh | The 9 July 551 AD Beirut Earthquake, Eastern Mediterranean Region[END_REF][START_REF] Elias | Active Thrusting Offshore Mount Lebanon; Source of the Tsunamigenic A.D. 551 Beirut-Tripoli Earthquake[END_REF]. Scientists suggest that the Mount-Lebanon Thrust might connect to the Yammouneh fault in the south and north of Lebanon.

The Yammouneh fault (YF) is through-going across the bend and it connects the southern and northern sections of the plate boundary, whereas the other fault branches are limited in extent [START_REF] Nemer | Coseismic Growth of Sedimentary Basins along the Yammouneh Strike-Slip Fault (Lebanon)[END_REF].

Overall, the long historic record of the region reveals a pattern of large earthquakes clustered in time. In more recent times, several earthquakes were recorded in the country with M<6. A M5.8 double shock struck Lebanon in 1956, followed by a M5.3 earthquake in 1997. The epicenter was estimated close to the northern inland segment of the Roum Fault [START_REF] Sadek | Updated Seismic Hazard for Lebanon and Implications on Micro-Zonation of the Greater Beirut Area[END_REF]. In early 2008, parts of South Lebanon were shaken by series of earthquakes of which the largest had a reported magnitude of 5.1 [START_REF] Huijer | Re-Evaluation and Updating of the Seismic Hazard of Lebanon[END_REF]. Present day, in Lebanon, instrumentally recorded seismicity is generally sparse within the Lebanese restraining bend. Despite the apparent lack of present-day small-seismicity on the north part of the YF (Figure 4-3), and large events on overall the fault systems, geodetic measurements and faulted landforms indicate that the YF accommodates most of the Africa-Arabia relative plate motion within the restraining bend [START_REF] Daëron | Constraints on the Post∼ 25-Ka Slip Rate of the Yammoûneh Fault (Lebanon) Using in Situ Cosmogenic 36Cl Dating of Offset Limestone-Clast Fans[END_REF]). This implies that the YF is at present accumulating strain, which underscores the concern for the earthquake potential of the YF. The mean recurrence period of 990-1260 yr for large earthquakes along of the YF [START_REF] Daëron | 12,000-Year-Long Record of 10 to 13 Paleoearthquakes on the Yammoûneh Fault, Levant Fault System, Lebanon[END_REF], and the last large earthquake occurred in 1202. It remains the most dangerous one because of its length as it crosses the Lebanese territory from south to north.

Seismic risk in Lebanon

Lebanon is located in an active tectonic environment where the seismic hazard is considered moderate to high though strong motion has never been recorded in Lebanon till now due to the presently infrequent large-magnitude seismicity. Hence the seismicity aspects need to be assessed with care. However, people are not ready for earthquakes and their consequences [START_REF] Beck | Crisis Mobility of Pedestrians: From Survey to Modelling, Lessons from Lebanon and Argentina[END_REF][START_REF] Cartier | Les Représentations Du Risque Sismique à Beyrouth et Bourj Hammoud[END_REF]; Beirut resilience plan 2017; [START_REF] Beck | Vulnerability to Earthquake of Beirut Residents (Lebanon): Perception, Knowledge, and Protection Strategies[END_REF]. As the recent instrumental seismicity is only moderate even in regions where the historical seismicity witnessed strong earthquakes, people ignore earthquake risks. Convincing the government and decision makers to establish a seismic policy and to invest in the assessment of the seismic risk is not an easy task, especially as Lebanon had unfortunately to cope over recent decades with other kinds of natural and non-natural hazards [START_REF] Verdeil | Atlas Du Liban[END_REF].

Moreover, during the Lebanese war that lasted many years, many buildings were constructed without any consideration of the building codes (Brax, Causse, and Bard 2016a). Despite the high seismic hazard that threatens Lebanon, until very recently there was no official building code taking into consideration this risk.

In order to assess the seismic hazard, the scientists often use the ground motion attenuation relationships GMPEs calibrated by a large set of recorded earthquakes worldwide. The GMPEs are function of the magnitude of the earthquake, the distance between the site and the rupture, and the characteristics of the site represented by the Vs30. [START_REF] Huijer | Upgrading the Seismic Hazard of Lebanon in Light of the Recent Discovery of the Offshore Thrust Fault System[END_REF] and [START_REF] Huijer | Re-Evaluation and Updating of the Seismic Hazard of Lebanon[END_REF] studied the seismic hazard of Lebanon, using the developed GMPEs applicable to the Mediterranean region [START_REF] Akkar | Empirical Equations for the Prediction of PGA, PGV, and Spectral Accelerations in Europe, the Mediterranean Region, and the Middle East[END_REF][START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods between 0.01 s and 10.0[END_REF][START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10[END_REF]. In their study, the attenuation characteristics of the ground motion are isotropic, i.e. they are independent of the location of the site relative to the source of energy release. They concluded that the peak ground acceleration on rock corresponding to 10 percent probability of exceedance in 50 years exposure time varies between 0.2 g and 0.3 g (Figure 4-4). They thus proposed to divide Lebanon into two seismic zones for regulatory applications:

-Zone I is the coastal zone between Saida and Tripoli that includes Beirut and the east part of Lebanon with effective ground acceleration on rock equal to 0.3 g.

-Zone II corresponds to the remaining part of the country, with an effective ground acceleration on rock of 0.25 g.

Consequently, The Lebanese Council of Ministers issued a Public Safety Decree in 2005 (decree No. 14293) that was revised in 2012 (decree No. 7964) to regulate safety procedures in buildings, installations, and elevators for the protection against fires and earthquakes; the seismic hazard is characterized by minimum horizontal peak ground acceleration on rock equal to 0.25 g. Accordingly, the Lebanese Standards Institution (LIBNOR) published the Lebanese norm relative to the protection from earthquakes NL135 (2012) that presents the seismic rules for the design of buildings and civil work structures. As in Lebanon there are two major structural engineering schools (French and American), the norm NL135 gives the possibility either to adopt the former French code PS92 (1995) or the American codes UBC97 (1997) or IBC ( 2009), for a level of seismic coefficient of (0.25 g) (Brax, Causse, and Bard 2016a). After 2012, Structures are designed to withstand a peak ground acceleration of PGA = 0.25 g. For large epicentral distances in comparison with the dimensions of the surface rupture (far-fault region), the magnitude and distance are dominating, and most probably enough to describe the source effect, and therefore GMPEs considering only the Mw as the source parameter are a good approach. However, Lebanon has small territory (10 452 km 2 ), within its ~200 km length and 30 to 80 km width. Hence, when dimensions of the surface rupture are of the same order as the distance to the location of interest (near-fault), the magnitude alone is not enough to characterize the earthquake source and the generated strong ground motion at a given site. The effect of other parameters such as the fault rupture speed, the average displacement, the rupture directivity and stress drop, and the heterogeneities along the fault area can be significant and have to be considered in strong ground motion predictions. Near-fault rupture simulations thus constitute a complimentary approach for the GMPEs in the far-fault region in order to assess seismic hazards for a large region.

Overview of the geology in Lebanon

The structure, topography, and history of the entire region are magnificent expressions of continental plates moving along a transform system. Lebanon is a mountainous country in the Levant region. It extends along the eastern coastline of the Mediterranean Sea. The landscape rises steeply from a narrow (5-15 km) shelf along the coast at the west where major cities in Lebanon are found (Beirut, Tripoli, and Tyre) to elevations of 2.5-3 km in the Mount-Lebanon Range (Figure 45, symbol 1). The Mount-Lebanon chain trends NNE and encloses the highest point in the Levant region (at 3,088 m asl). East of the Mount-Lebanon Range, the landscape drops steeply over sparsely vegetated slopes and reaches the Bekaa Valley (Figure 4-5, symbol 2), a major agricultural region. The altitudes of the Bekaa valley floor range from 500 to 1,000 m asl. Continuing east, the Anti-Lebanon Range trends NE, and rises to ~2,600-2,800 m asl over rolling hill slopes that form Lebanon's border with Syria. Hence, the Bekaa Valley separates the two mountain chains, the Mount-Lebanon chain to the west and the Anti-Lebanon chain to the east. A major structural feature on the eastern Mount-Lebanon mountain slope is the Yammouneh Faullt, which seperates a steep narrow strip on high mountain area from mountain slopes with less pronounced relief east of the fault. In the southern Mount-Lebanon range, two faults branch out from the Yammouneh Fault system in northwest and northeast direction respectively: the Roum Fault crossing the Lebanon mountains toward the Mediterrranean Sea coast and the Rachaya and Serghaya Faults extending into the Anti-Lebanon massif. Along these fault systems, several intermountain basins subsided between the mountain chains.

The landscape of Lebanon is principally dominated by limestone and sandstone, dating back to the Early Jurassic (Figure 45). The early Jurassic Kesrouane Limestone is a structurally complex and often thick limestone unit that underlies much of the mid-elevation Mount-Lebanon region. The Chouf Sandstone unit, which overlies Jurassic limestone, is widespread and highly fractured. Limestones (Sannine and Mdairej) form the modern Mount-Lebanon and Anti-Lebanon ranges [START_REF] Grant | Multimodal Method for Coseismic Landslide Hazard Assessment[END_REF]. The relatively high precipitation on the mountains and the high infiltration rates on widely exposed karst surfaces and karstified Jurassic limestones make Mount-Lebanon range certainly with the most productive aquifer system of the Arabian Plate. Surface drainage from streams on the western slopes of the Mount-Lebanon mountains builds a dense network of coastal rivers that is directed to the Mediterranean Sea. The 170 km long Litani river (Figure 4-6) runs from its headwaters near Baalbek in the Bekaa plain to the Mediterranean Sea coast in southern Lebanon. The 487km Orontes (Nahr el Aasi) river begins in the northern Bekaa plain and flows through Syria and Turkey before entering the Mediterranean Sea. The rivers drain from the eastern slopes of the Anti-Lebanon mountains into the closed basin of the Damscus plains in Syria.

Apart from the surface geology and active faults mapping [START_REF] Grant | Multimodal Method for Coseismic Landslide Hazard Assessment[END_REF][START_REF] Daeron | 12,000-Year-Long Record of 10 to 13 Paleoearthquakes on the Yammouneh Fault, Levant Fault System, Lebanon[END_REF], we only know little about the geological and seismic structure at depth down to the Moho, which limit our understanding of the interaction of major faults at depth, the thickness of the seismogenic zone, the seismic cycle, and therefore limit our capacity to predict the strong ground motion and consequently the seismic hazards.

Two onshore seismic cross-sections up to 5 km in depth, the Batroun-Ainata profile trending almost east -west across northern Mount-Lebanon, and the Aley-Barr Elias profile trending almost east-west across central Mount-Lebanon have been recently performed by [START_REF] Nader | The Geology of Lebanon[END_REF] and are shown in Figure 4-7. These first 2D seismic reflection profiles have demonstrated that the North Mount Lebanon and the Beqaa plain are well distinguished, and that there is a continuity of the exposed structures into the subsurface.

The Earth structure and extent of the underlying sedimentary succession onshore Lebanon can only be estimated from regional correlations and thickness extrapolations from the surrounding countries. In Israel, there is a strong interface at 3-5km; in Jordan, the interface is shallower (at 2km). Both countries have another interface at 8-10 km depth. Deeper in the Earth, the structure is homogeneous, until reaching the Moho (Aldersons et al. 2003;[START_REF] Pinsky | Imaging Heterogeneity of the Crust Adjacent to the Dead Sea Fault Using Ambient Seismic Noise Tomography[END_REF]. Those interfaces of the stratigraphy in depth are characterized by seismic wave-speed changes. The structure of the offshore lithosphere in the Levant basin located offshore Egypt, Cyprus, Israel, Lebanon and Syria showed that the Moho at the Mediterranean coast is about 20-23 km depth [START_REF] Inati | Lithospheric Architecture of the Levant Basin (Eastern Mediterranean Region): A 2D Modeling Approach[END_REF]. The Moho depth increases strongly from west (~25 km) toward the east (~39 km), indicating a small asymmetric topography in the Moho discontinuity [START_REF] Koulakov | Moho Depth and Three-Dimensional P and S Structure of the Crust and Uppermost Mantle in the Eastern Mediterranean and Middle East Derived from Tomographic Inversion of Local ISC Data[END_REF][START_REF] Khair | Crustal Structure of the Northern Levant Region: Multiple Source Werner Deconvolution Estimates for Bouguer Gravity Anomalies[END_REF].

TOMOGRAPHY OF LEBANON USING SEISMIC AMBIENT NOISE

A seismic network is primarily installed to detect an earthquake when it occurs. Even outside the periods of seismicity, seismic stations register ambient noise as seismic waves generated permanently at the surface of the globe, due to the interaction between the atmosphere, the oceans and the solid earth. Seismic noise records are not unusable: they allow imaging the Earth. Those ambient seismic waves propagating from one point in the Earth crust towards the surface hold the fingerprints of the medium of propagation. Technically, the time cross-correlation function of the ambient noise computed at two distant receivers converges to the complete Green's Function between these two receivers. Green's function contains the complete information that characterizes the medium of propagation, and thus is the essential ingredient for imaging. This passive imaging technique is applied to see inside the Earth crust in Lebanon. This chapter is in preparation for publication.

-

------------------------------------------------------------------------------------------------------ ------------------------------------------------------------------------------------------------------- R. FAYJALOUN 1 , L. YANG 1 , L. STEHLY 1 , C. CORNOU 1 , C. VOISIN 1 1 Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, 38000 Grenoble, France -------------------------------------------------------------------------------------------------------

Abstract

The Dead Sea Fault is a 1000 km transform plate boundary in the Middle East, presenting an impressive tectonic feature in the region, however posing a considerable seismic risk. Within Lebanon, the fault splays into four fault branches, and one of them crosses through the whole country. The crustal model is a vital prerequisite to predict the seismic risks in Lebanon. Because of its location on the coast and its mountainous terrain, Lebanon's physical and geological landscape is complex and varies markedly with short distances. What about deep inside the earth? In this work, we present the first 3D imaging of the earth crust (in terms of shear wave velocity VS) in Lebanon, using the seismic ambient noise propagating in the medium.

We use cross-correlations of the ambient seismic noise recorded by 21 stations. The cross-correlation functions converge to the Green's Function between any pair of stations. Rayleigh waves group velocity is therefore extracted at different periods in the 1-25 s period band using the frequency-time analysis procedure and then inverted to compute Rayleigh wave velocity maps using the adaptive grid inversion approach. Finally, these maps are inverted using iterative methods to compute the VS maps in the depth range 1-40 km. The VS model is discussed in the light of known geological units in Lebanon and compared to the VS models of the neighboring countries. A first order estimation of the Moho interface is also presented.

Key words: Seismic ambient noise; Passive imaging; 3D VS profile, Lebanon, Moho.

Introduction

The Dead Sea fault transform (DSFT) is the most impressive tectonic feature in the Middle East. It is a left lateral transform fault, separating the Arabian plate and the Sinai plate, which transfers sea floor spreading in the Red Sea in the south to the Taurus-Zagros collision zone in Turkey and Iran in the north (Figure 5-1). The DSFT has been active since the Miocene [START_REF] Garfunkel | Internal Structure of the Dead Sea Leaky Transform (Rift) in Relation to Plate Kinematics[END_REF][START_REF] Garfunkel | The Structure of the Dead Sea Basin[END_REF] with movement continuing today over a length of more than 1000 km with a cumulative displacement of 105 km. Large earthquakes are known to have occurred along this fault zone with recurrent magnitudes in the range of 6-7.5 over the historical period [START_REF] Ambraseys | Value of Historical Records of Earthquakes[END_REF]Ben-Menahem 1991a;[START_REF] Shapira | A New Estimate for the Epicenter of the Jericho Earthquake of 11 July 1927[END_REF][START_REF] Daeron | 12,000-Year-Long Record of 10 to 13 Paleoearthquakes on the Yammouneh Fault, Levant Fault System, Lebanon[END_REF][START_REF] Beon Maryline | Slip Rate and Locking Depth from GPS Profiles across the Southern Dead Sea Transform[END_REF][START_REF] Lefevre | Slip Deficit and Temporal Clustering along the Dead Sea Fault from Paleoseismological Investigations[END_REF]. The DSFT poses a considerable seismic hazard to Syria, Lebanon, Israel, Jordan, and the Palestine Territories. A detailed knowledge concerning both the regional motions of the tectonic plates and the crustal structure allows to understand the source rupturing and the propagation properties in the crustal model and thus, to better estimating the seismic hazard. Previous seismic imaging has established solid baseof knowledge of crustal structure for different regions along the fault, mainly in Jordan, Israel, the Dead Sea and offshore, with some seismic profiles in Lebanon (e.g. [START_REF] Khair | Crustal Structure of the Northern Levant Region: Multiple Source Werner Deconvolution Estimates for Bouguer Gravity Anomalies[END_REF]Aldersons et al. 2003;[START_REF] Koulakov | Moho Depth and Three-Dimensional P and S Structure of the Crust and Uppermost Mantle in the Eastern Mediterranean and Middle East Derived from Tomographic Inversion of Local ISC Data[END_REF][START_REF] Pinsky | Imaging Heterogeneity of the Crust Adjacent to the Dead Sea Fault Using Ambient Seismic Noise Tomography[END_REF][START_REF] Nader | The Geology of Lebanon[END_REF][START_REF] Inati | Lithospheric Architecture of the Levant Basin (Eastern Mediterranean Region): A 2D Modeling Approach[END_REF]. Lebanon is one of the countries crossed by the DSFT and faces large seismic risk. Within the Lebanese bend, the DSTF fault splits into four main branches, with clear evidence for slip partitioning: the Yammouneh, the Roum and Rachaya -Serghaya are left-lateral strike-slip faults, and the offshore Mount-Lebanon Thrust fault.

Yet, it lacks the crustal model to predict seismic waves propagation, a vital prerequisite to seismic hazard and risk assessment. How do the rock elastic properties change along the fault and with the depth in Lebanon? At what depth does the crust-mantle boundary manifest?

Over the last two decades, the ambient noise cross-correlation technique was developed, allowing to retrieving the elastic properties at the crustal and shallow structures between pairs of receivers at a low cost tool, using only randomly scattered ambient noise data recorded for many months. Indeed, cross-correlation of the ambient noise recorded at receivers returns an approximation to the Green's function characterizing the medium between the receivers [START_REF] Weaver | Ultrasonics without a Source: Thermal Fluctuation Correlations at MHz Frequencies[END_REF][START_REF] Shapiro | Emergence of Broadband Rayleigh Waves from Correlations of the Ambient Seismic Noise[END_REF][START_REF] Draganov | Seismic Interferometry: Reconstructing the Earth's Reflection Response[END_REF][START_REF] Wapenaar | Green's Function Representations for Seismic Interferometry[END_REF]. At frequencies less than 1 Hz, most of ambient seismic noise is generated by natural atmospheric and oceanic forces at the surface [START_REF] Longuet-Higgins | A Theory of the Origin of Microseisms[END_REF][START_REF] Friedrich | Ocean-Generated Microseismic Noise Located with the Gräfenberg Array[END_REF][START_REF] Bonnefoy-Claudet | The Nature of Noise Wavefield and Its Applications for Site Effects Studies: A Literature Review[END_REF][START_REF] Stehly | A Study of the Seismic Noise from Its Long-range Correlation Properties[END_REF][START_REF] Landès | Origin of Deep Ocean Microseisms by Using Teleseismic Body Waves[END_REF]. Therefore, the surface waves dominate the Green's function extracted from the noise cross-correlations and gives us the possibility of measuring the dispersive characteristics of surface waves between any pair of stations, with the lower frequency component characterizing deeper structures. The passive seismic imaging of crustal structure by using ambient noise cross-correlation was first applied for Southern California [START_REF] Shapiro | High-Resolution Surface-Wave Tomography from Ambient Seismic Noise[END_REF][START_REF] Karim | Extracting Time-domain Green's Function Estimates from Ambient Seismic Noise[END_REF], and was applied since then in various parts of the world: e.g. the southeastern Tibetan plateau in China [START_REF] Yao | Surface-Wave Array Tomography in SE Tibet from Ambient Seismic Noise and Two-Station Analysis -I. Phase Velocity Maps[END_REF], in western Europe [START_REF] Stehly | Tomography of the Alpine Region from Observations of Seismic Ambient Noise[END_REF]), northern Finland [START_REF] Poli | Noise Directivity and Group Velocity Tomography in a Region with Small Velocity Contrasts: The Northern Baltic Shield[END_REF], in the dead sea fault area [START_REF] Pinsky | Imaging Heterogeneity of the Crust Adjacent to the Dead Sea Fault Using Ambient Seismic Noise Tomography[END_REF]) and even deep 3D structure of the Earth [START_REF] Nishida | Global Surface Wave Tomography Using Seismic[END_REF].

In this paper, we use the cross-correlation technique to study the lithosphere of Lebanon. The pseudo 3D Shear wave velocity structure is obtained in a two-step approach, the first step comprises calculating fundamental group velocity dispersion curves in the period range of 1-25 s; the second step involves an inversion of the dispersion curves in order to obtain local 1D S-wave velocity depth profiles.

The final pseudo 3D S-wave velocity model is obtained by interpolation between the individual 1D S-wave velocity profiles and compared to the existing profiles of the surrounding regions.

Geological background

Lebanon extends along the eastern coastline of the Mediterranean Sea and forms the eastern margin of the northern part of the offshore Levant Basin. The area of Lebanon is roughly rectangular in shape, with approximately 200 km long, becoming narrower toward the south (Figure 4-5). Its widest point is 88 km, and its narrowest is 32 km. A major feature of Lebanese topography is the alternation of lowland and highland that runs generally parallel with a north-to-south direction.

Because of its mountainous terrain, Lebanon's physical geography is complex and varied. Landforms, soil and vegetation change markedly within short distances. The country is divided into four physiographic units, from west to east:

1) the coastal plain runs along the Mediterranean shore, forms a narrow and discontinuous strip of river-deposited alluvium and marine sediments which alternate with rocky beaches and sandy bays; 2) the mountain range called the Mount-Lebanon range rises to alpine heights southeast of Tripoli (up to 3088 m). Many springs emerge from the permeable limestone of the mountains; 3) the Beqaa valley is a central highland separating the Mount-Lebanon and the Anti-Lebanon mountains to the east. It is narrow in the south and wider northward, and is composed mainly of alluvial deposits from mountains on either side; 4) and the Anti-Lebanon mountains form the eastern border with Syria.

The geology of Lebanon is principally dominated by limestone and sandstone, dating back to the Early Jurassic (Figure 45). The early Jurassic Kesrouane Limestone is a structurally complex and often thick limestone unit at mid-elevation of the Mount-Lebanon region. The Chouf Sandstone unit, which overlies Jurassic limestone, is widespread and highly fractured. Limestones (Sannine and Mdairej) form the modern Mount-Lebanon and Anti-Lebanon ranges [START_REF] Grant | Multimodal Method for Coseismic Landslide Hazard Assessment[END_REF]. The relatively high precipitation on the mountains and the high infiltration rates on widely exposed karst surfaces and karstified Jurassic limestones make Mount-Lebanon range certainly with the most productive aquifer system of the Arabian Plate. Surface drainage from streams on the western slopes of the Mount-Lebanon mountains builds a dense network of coastal rivers that is directed to the Mediterranean Sea. The 170 km long Litani river runs from its headwaters near Baalbek in the Beqaa plain to the Mediterranean Sea coast in southern Lebanon. The 487 km Orontes (Nahr el Aasi) river begins in the northern Beqaa plain and flows through Syria and Turkey before entering the Mediterranean Sea. The rivers drain from the eastern slopes of the Anti-Lebanon mountains into the closed basin of the Damscus plains in Syria. Another major feature of Lebanon is that it is crossed by the continental plate boundary-the Dea Sea Transform Fault DSFT. The DSTF N-S striking fault system separates the Arabian plate from the African (or Sinai) plate (Figure 5-1-a). Within the Lebanese region, the DSTF splays into four fault branches (Figure 5-1-b), and one of them, the strike slip Yammouneh Fault, is through-going across the whole country (e.g. [START_REF] Walley | A Braided Strike-Slip Model for the Northern Continuation of the Dead Sea Fault and Its Implications for Levantine Tectonics[END_REF][START_REF] Gomez | Holocene Faulting and Earthquake Recurrence along the Serghaya Branch of the Dead Sea Fault System in Syria and Lebanon[END_REF][START_REF] Nemer | Evidence of Coseismic Ruptures along the Roum Fault (Lebanon): A Possible Source for the AD 1837 Earthquake[END_REF]. The Yammouneh fault (~170 km) trends NE-SW and constitutes the main continuity of the Dead Sea Fault DSTF separating the Arabian and the African plates. Rachaya-Serghaya Fault traces (~45 km and ~100-150 km respectively) run almost parallel to each other along the Anti-Lebanon mountain range. The Roum fault (~35 km) branches from the DSTF in south Lebanon and runs along the south-western boundary of the Mount Lebanon range. And last but not least, the newly identified Mount Lebanon Thrust (~150 km), plunges under the western side of the Mount Lebanon range. Its surface trace lies mainly offshore, cutting the seabed at not more than 8 km from the coast of central Lebanon between Tripoli in the North and Saida in the South [START_REF] Elias | Active Thrusting Offshore Mount Lebanon; Source of the Tsunamigenic A.D. 551 Beirut-Tripoli Earthquake[END_REF][START_REF] Huijer | Upgrading the Seismic Hazard of Lebanon in Light of the Recent Discovery of the Offshore Thrust Fault System[END_REF].

Ambient noise cross-correlation and 3D tomography of Lebanon

In the following, the location of the seismic stations used for this study as well as the period at which the ambient noise was recorded is presented. Then, the different steps detailing the pre-processing of the recorded noise, the cross-correlation, the extraction of the surface Rayleigh wave group velocity and finally the inversion of the shear-wave velocity model is detailed.

Stations distribution and period of recordings

Within the framework of the LIBRIS project funded by the French National Research Agency, twenty-one broadband stations were distributed in Lebanon from January 2011 to mid September 2012: 18 seismological stations were provided by the French mobile seismological pool INSU/RESIF (LU and LUK stations in Figure 4-5) and 3 by the national permanent Lebanese network operated by CNRS-L (BHL0, HWQL, BEYL stations in Figure 45). These seismological stations were equipped with broadband three-component velocimeters (CMG40 sensors having a cut-off period of 30s) connected to Nanometrics Tauris, Kephren or Geosig digitizers (Appendix 5-1). These stations continuously recorded ground motion at a sampling rate of 200 Hz. During the experimental period, some stations failed or were stolen.

Appendix 5-2 shows the period of the data available for each station.

Ambient noise data processing and cross-correlation

The data were next prepared to reduce the effects of transient events and enhance the ambient noise. The continuous Z-component records are chopped into 24 hours time windows. Then the Fourier amplitude spectra of the records are whitened using the comb filter in the period bands of [0.25-1], [1-2], [2-4], [4-10], [10][11][12][13][14][15][16][17][18][19][20], s, in order to prevent the records to be dominated by energetic events such as storms or earthquakes. One example of raw ambient noise as well as its preprocessed resultant for day Julian day 200 in 2011 at stations LU12 are shown in Appendix 5-3. We then correlate velocity time series recorded on the vertical components day-by-day in 4-hour windows. All correlation functions for each pair are then stacked. This is equivalent to directly cross-correlating the whole period of records yielding to the Rayleigh wave Green's function. The correlation 𝐶 𝐴𝐵 between the stations 𝐴 and 𝐵, normalized by the energy of the 2 segments 𝑆 𝐴 and 𝑆 𝐵 yields the same Green's Function which would be gleaned at A if B were a receiver, 𝑇 𝑖 being the integration time, 𝑡 is the lag time [START_REF] Snieder | Imaging with Ambient Noise[END_REF]:

𝐶 𝐴𝐵 (𝑡) = ∫ 𝑆 𝐴 (𝑡) 𝑆 𝐵 (𝑡 + 𝜏) 𝑑𝜏 𝑇 𝑖 0 √ ∫ 𝑆 𝐴 2 (𝜏) 𝑑𝜏 ∫ 𝑆 𝐵 2 (𝜏) 𝑑𝜏 𝑇 0 𝑇 𝑖 0
(5-1)

Figure 5-3 shows the ZZ cross-correlation of records between the stations LU12 and LU08, the stations being separated by an interstation distance of 143 km. The noise cross-correlation is shown at different range of periods. For each range of period, the positive time of the correlation corresponds to the causal Green's function of the medium between LU12 and LU08, and the negative time corresponds to its anticausal counterpart (i.e. the Green's function between LU08 and LU12).

Rayleigh waves group velocity measurements

Rayleigh waves dispersion curves are evaluated from the emerging Green's functions using the frequency-time analysis FTAN [START_REF] Herrmann | Some Aspects of Band-Pass Filtering of Surface Waves[END_REF][START_REF] Bhattacharya | Higher Order Accuracy in Multiple Filter Technique[END_REF][START_REF] Ritzwoller | Eurasian Surface Wave Tomography: Group Velocities[END_REF], by measuring the surface wave group velocity at different periods between all cross-correlated station pairs. The FTAN consists of Gaussian narrowband (bandpass) filtering the signal around each targeted period T, T is the central period of the filter. We identify at each period the maximum of the signal envelope amplitude (highlighted in blue in Figure 5-3) and find the corresponding group velocity u = D/t max , where D is the distance between the stations and t max is equal to the travel time of the maximum of signal amplitude envelope, within the range of 1-5 km/s (highlighted in red). Beyond the surface waves (precisely starting at t max + 3 * T), exists the noise (highlighted in yellow) that is used next to compute signal-to-noise ratio and therefore defines the quality of the constructed Green's functions.

Criteria for selecting the measured Rayleigh wave group velocity measured

We carefully select the group velocity based on the following quality tests. The first criterion is the signal-to-noise ratio (SNR), defined as the ratio of the amplitude of the Rayleigh wave and the standard deviation of the noise that follows the reconstructed Rayleigh wave. The SNR is measured separately on the positive and negative correlation time. Waveforms with SNR ≤3 are discarded from our analysis.

The second criterion is the discrepancies of the values of the velocities found in the causal part and acausal part of the cross-correlation. Both, they sample the same media and are expected to exhibit the same velocities and dispersion characteristics assuming an isotropic distribution of sources. For each path, we have one correlation corresponding to Rayleigh waves that leads to 2 measurements of the Rayleigh wave group velocity since we consider both the positive and negative correlation times. We keep only the correlations if the Rayleigh waves group velocities measured on the positive and negative correlation time differ by less than 0.5 km/s. The retained measurements are then averaged for every path to get the final group velocity.

The third criterion is based on the distance between the stations. We reject group velocity estimates that correspond to interstation distances smaller than one wavelength (λ=VR x T) otherwise the Rayleigh waves are not fully constructed.

The selected Rayleigh wave group velocities between the pairs of stations are illustrated in Figure 5-4 at different periods. The Rayleigh dispersion curve, that is the surface wave group velocities as function of the period between 1 and 25 s, is consequently constructed (Appendix 5-4).

Robustness of Rayleigh wave group velocity measurements

Northeast and south of Lebanon do not benefit from a good coverage of rays representing the computed group velocity. In the northwest, there are no stations recording in the area other than LU02. In the south, stations LU18 and LUK5 recorded for only 3 months; their periods of recordings do not always correspond to the same period of recordings of other stations. However, a fair convergence of the SNR for most of the Green's functions constructed is reached within 3 months (Appendix 5-5).

The azimuthal distribution of VR would depend on the geometry of the network and on the distribution of the noise sources. The geometry of the network depends on the distribution of the stations which, itself, depends on the geometry of the country and its orientation that is already NE-SW oriented with ~30 degrees from the North. In our case, the geometry of the network governs the azimuthal distribution of VR; Even though the sea constitutes a major source of the noise generation at the westside coast, its impact on the azimuthal distribution of VR is not clear (Appendix 5-6).

2D Rayleigh waves group velocity maps

The Rayleigh wave dispersion curves are inverted to obtain discretized group velocity maps at different periods. We implement an adaptive multi-grid approach for the discretization of our area, in order to account for non-uniformity in data coverage. The discretization is fully controlled by the ray distributions between the stations. The area is discretized with a rough grid size R; if the cell contains at least one ray, it is subdivided into 4 cells; each cell that contains a minimum of 5 rays/cell is divided into 4 cells, for 2 iterations, and therefore the finest grid size would be m = R/2 3 .

The standard forward problem can be written in tensor notation:

𝐃 = 𝐏 * 𝐌 (5-2)
where D is the data vector whose elements are the difference between the measured group travel times (t meas ) and the group travel times computed from the initial model for each path (t 0 ). P is the matrix that represents the travel path of the surface wave within each cell of the initial model, and vector M = (1/u 0 -1/u) represents the slowness perturbations, u being the velocity obtained after inversion, and u 0 the initial group velocity; the starting group velocity model is 2.7 km/s. The process of inversion is repeated for all periods of interest. The seismic inverse problem is smoothed (regularized by a first-order roughness damping), in order to suppress instabilities in the solution caused by noisy and incomplete data [START_REF] Schaefer | Adaptively Parametrized Surface Wave Tomography: Methodology and a New Model of the European Upper Mantle[END_REF]. This means that strong small-scale heterogeneities would not be properly accounted for [START_REF] Boschi | High-and Low-resolution Images of the Earth's Mantle: Implications of Different Approaches to Tomographic Modeling[END_REF]. Different discretization grid sizes and roughness parameters are tested (Appendix 5-7 and Appendix 5-8) to choose the parameters that give robust results. In the following, R = 0.4 degrees and the corresponding m = 0.05 degrees discretization parameters are adopted for a roughness damping coefficient of d = 0.4.

Figure 5-5 portrays the corresponding results of the group velocity tomography at different periods. The group velocity ranges between 1.5 and 3.5 km/s. The tomography reveals relatively two different zones persisting for the period range 1-7 s: lower velocity zone in the southwestern part of Lebanon with respect to the northern part. For periods larger than 7 s, NE-SW strips are well defined up to 20 s. For larger periods, a large group Rayleigh velocity appears in the south of Lebanon and extends to the north.

Resolution test of the tomography of Rayleigh wave velocities

The main factors that influence the results of the tomography of VR at different periods are: 1) the density of paths, 2) their azimuthal distribution and 3) the size of the grid and the smoothing parameter.

In order to evaluate the resolution of the 2-D group velocity maps, the checkerboard test is a useful alternative that depicts the ability of tomographic inversion to resolve structural details in the earth at different periods [START_REF] Boschi | Multiple Resolution Surface Wave Tomography: The Mediterranean Basin[END_REF][START_REF] Spakman | A Tomographic View on Western Mediterranean Geodynamics[END_REF]. The checkerboard test procedure is to superimpose a small perturbation signal onto the initial model, compute synthetic arrival-time data, and then invert the synthetic arrival-time in the same manner as the actual data. The ability of the tomographic method to quantitatively recover the perturbed model is then an estimate of the sensitivity of the original inversion of real data to recover similar details in real earth. For different periods, we use a synthetic velocity model deviating by ±10% from an input constant model. Synthetic Rayleigh wave travel times over the interstation paths are now computed according to this synthetic model and used for the tomography procedure. Appendix 5-9 shows the results of the checkerboard tests for group velocity maps at periods 2, 8 and 20 s using anomalies of size 0.2• (∼22 km) and 0.4• (∼44.4 km). The velocity anomalies of 0.4• are better recovered. At each period, the best resolution is achieved around Beirut and in the North. Towards peripheral areas, the resolution degrades and the smearing increases due to lower ray-paths density. 2.9

V R (g) (km/s)

3D Shear waves velocity inversion and related uncertainty

The area in Lebanon is discretized into cells of dimensions 0.05 x 0.05 degrees for the tomography mapping at different periods. For each cell, we construct the dispersion curves, i.e. Rayleigh group velocity as function of periods, as shown in Appendix 5-10. The dispersion curves are then inverted to obtain corresponding shear velocity profile, i.e. the shear wave speed as function of depth, at each grid point, using the conditional neighborhood algorithm [START_REF] Wathelet | An Improved Neighborhood Algorithm: Parameter Conditions and Dynamic Scaling[END_REF].

Ground model parameterization and the inverted 3D shear waves model

The ground model parameterization used in the inversion is defined according to known crustal characteristics (Primary wave speed Vp and shear wave speed VS) of the surrounding regions. The measured VS profiles (inferred from Vp profiles assuming a poisson's ratio of 0.33 at the surface and 0.25 otherwise) from Israel, Jordan and offshore tomography (Figure 5-6) show that the wave speed is low in the first 2 to 5 km, then becomes relatively constant until reaching the Moho at 20-35 km depth (Aldersons et al. 2003;[START_REF] Pinsky | Imaging Heterogeneity of the Crust Adjacent to the Dead Sea Fault Using Ambient Seismic Noise Tomography[END_REF][START_REF] Khair | Crustal Structure of the Northern Levant Region: Multiple Source Werner Deconvolution Estimates for Bouguer Gravity Anomalies[END_REF][START_REF] Koulakov | Moho Depth and Three-Dimensional P and S Structure of the Crust and Uppermost Mantle in the Eastern Mediterranean and Middle East Derived from Tomographic Inversion of Local ISC Data[END_REF][START_REF] Inati | Lithospheric Architecture of the Levant Basin (Eastern Mediterranean Region): A 2D Modeling Approach[END_REF].

Consequently, the parameterization of the ranges of Vp and VS values is indicated in Appendix 5-11. The thickness of the first layer is fixed to be larger than 1 km to account for the lack of dispersion estimates at frequencies larger than 1 Hz. As for the VS calculation, we proceed in two steps. The first one consists in searching for the "best misfit" profile, letting the inversion algorithm looking for the "absolute" minimum misfit between actual and forward modeled dispersion curves. Best misfit shear-wave profiles are obtained at each point of the grid after 50050 iterations. 3D representation of the VS model is shown in Surface wave inversion provides very non-unique ground models, especially in the absence of already a priori robust knowledge of the VS ground structure, as it is the case for Lebanon. In order to evaluate the uncertainty of the inverted VS profiles, we perform a new inversion by arbitrary considering standard deviation of 5% on the group velocity estimates and using the acceptable misfit concept [START_REF] Lomax | Identifying Sets of Acceptable Solutions to Non-Linear, Geophysical Inverse Problems Which Have Complicated Misfit Functions[END_REF][START_REF] Hollender | Characterization of Site Conditions (Soil Class, V S30, Velocity Profiles) for 33 Stations from the French Permanent Accelerometric Network (RAP) Using Surface-Wave Methods[END_REF]) in order to get a suite of VS profiles that explain the measured group velocities within their uncertainty bound. From these sets of equivalent VS profiles, we then compute at each depth the coefficient of variation (COV). Small COV witness well constrained VS estimates while large COV denotes lack of resolution or uncertainty in the seismic VS interface depth.

The COV of the proposed 3D VS model is shown in

Comparison of inverted V s model to previous research outcomes

The spatial variation of VS at various depths, the Moho interface depth, as well as the spatial variation VS at two vertical cross-sections are compared to the geology of Lebanon and to the known structure models in the surrounding region. In the south, a coastal plain of about 1-2 km width extends over a long stretch of the coast between Saida and Tyr (Figure 5-1). These coastal plains are composed of unconsolidated alluvium deposits that explain low Vs at ~2.1 km/s. Deeper in the crust between 2 and 3 km, a high shear-wave speed zone appears at the very north limits of the country. VS is around 3.2 km/s and could be associated to the Basalt patch of the volcanic massif that extends over the Lebanese-Syrian border area [START_REF] Grant | Multimodal Method for Coseismic Landslide Hazard Assessment[END_REF]. At z = 2-3 km, there is a transition of VS value along the Yammouneh fault trace. At z = 3 km, the westside of the Yammouneh Fault has VS ~3 km/s versus VS ~2.8 km/s on the eastside (the corresponding COV values are low). At largest depths (between 5 and 15 km), velocities are rather homogenous throughout the country and well resolved (COV less than 5%, Figure 5-9) with some changes however at the westside marking the intrusion of a new zone with larger VS (3.5 km/s), however with large COV up to 10%.

First order estimation of the Moho depth?

Although our inversion is poorly constrained at large depth (>15-20 km) due to less density of ray-paths at larger periods (Figure 5-4), and to the inverted VS profiles with larger COV values at large depths (Figure 5-9), our pseudo 3D model (Figure 5-7) clearly indicates a west-east dipping of a seismic interface from 20 km to 40 km, which could be associated to the Moho. This asymmetry is also reflected in the variation of the VS profile in depth getting larger . [START_REF] Segev | The Structure, Isostasy and Gravity Field of the Levant Continental Margin and the Southeast Mediterranean Area[END_REF]) studied the deep lithosphere of the Levant area, and showed that the Moho depth is increasing from 20 km offshore, to 40 km to the east. [START_REF] Khair | Crustal Structure of the Northern Levant Region: Multiple Source Werner Deconvolution Estimates for Bouguer Gravity Anomalies[END_REF] and [START_REF] Koulakov | Moho Depth and Three-Dimensional P and S Structure of the Crust and Uppermost Mantle in the Eastern Mediterranean and Middle East Derived from Tomographic Inversion of Local ISC Data[END_REF] found that the Moho increases strongly from west toward the east from 26 to 39 km, indicating a small asymmetric topography in the Moho discontinuity, and that Lebanon region has the Moho depth at ~30 km. Likewise, [START_REF] Inati | Lithospheric Architecture of the Levant Basin (Eastern Mediterranean Region): A 2D Modeling Approach[END_REF]) studied the structure of the lithosphere underlying the easternmost Mediterranean region (Egypt, Cyprus, Israel, Lebanon and Syria), and showed that the Moho offshore under the Arabian plate is 35-40 km deep, mainly in the south of Lebanon, and becomes shallower towards the Mediterranean coast.

Vertical cross-sections of VS

Finally, we study the lateral variation of VS along two cross-section profiles, labeled 'Profile 1' and 'Profile 2' in Figure 5-10 (a). The two vertical profiles are compared to the geological cross-sections 'A' and 'B' of [START_REF] Nader | Geological Interpretation of 2d Seismic Reflection Profiles Onshore Lebanon: Implications for Petroleum Exploration[END_REF] delimited in white in Figure 5-10 (a) and presented in (f) and (g). The lateral change under the Yammouneh Fault trace observed in the tomography is clearly perceived using the seismic reflection in profile 'B'. The interface between the Jurassic and the Triassic geological units is portrayed in black dots in Figure 5-10 (f) and (g) and is superposed on the VS section in Figure 5-10 (d) and (e) as the dashed black curves. Another dashed curves are also illustrated on each profile that refers the interface between the Jurassic and the Triassic units corrected by the surface topography. The dashed curves are in good agreement with VS interface from 2.65 km/s to 3 km/s inferred from tomography.

Even though the two methods, i.e. the seismic reflection and the ambient noise tomography, image the lithosphere at different spatial resolution scales, the overall variation of underground structure in both studies match for the important interfaces between different geological units. VS, that depends on the type of rocks, but also on the compaction, and the geological units are usually correlated, however there is no well-defined connection between them. Even though the comparison is not forthright between the two methods of tomography, the matching of the interfaces adds a validation value on the results of our tomography.

Conclusion

We present the 3D VS model of Lebanon's crust up to 40 km using passive seismic imaging technique, consisting of using seismic ambient noise recorded between 2011 and 2012 at 21 broadband stations in Lebanon. Rayleigh wave Green's functions were retrieved by correlating the noise records between all station pairs. We used these measurements to perform Rayleigh wave tomography at periods ranging from 1 to 25 s, and finally inverted a pseudo 3D VS model. Important interfaces are depicted at 2-3 km, 8-15 km, and larger than 20 km. The traces of the Yammouneh fault and the coastal basins are depicted at subsurface up to ~4 km. Homogeneous structures throughout the country are outlined deeper up to ~ 20 km. The deepest interface is propably the Moho, dipping from the west to the east between 20 and 40 km. The VS profile is consistent in average with VS profiles from the surrounding countries. The 3D VS pseudo model would help to, first, better comprehend the geodynamics in the region, and second, to better model wave propagations for seismic hazard assessment.

The surface topography is not taken into account during the inversion processes.

The resolution of the obtained VS model depends mostly on the density of the stations, the period of their recordings, and the inherent properties of the surface waves. Although our tomography allows retrieving the main geological and tectonic features of Lebanon, it would be interesting to install denser coverage stations to get a finer spatial resolution over the shallow structural part, as well as initiate body waves tomography to depict more precisely the Moho depth.

Appendix

Appendix 5-1: The characteristics of the seismological stations used for this study. * CMG40 are velocity-meters having their greater sensitivity in the frequency band from 0.03 Hz to 40 Hz.

YAMMOUNEH FAULT (LEBANON)

We have come full circle! In chapter 2, we inferred the source parameters that most affect the ground motion in the near fault, and in chapter 5, we derived the elastic properties of the propagation medium in Lebanon. We have now the main ingredients to simulate ground motion up to 1 Hz in Lebanon in the near fault. To simulate ground motion in a broad frequency range (up to 10 Hz), a hybrid stochastic model for near fault ground motion taking into account the directivity effects is presented in this chapter.

-------------------------------------------------------------------------------------------------------

Abstract

This work evaluates the ground motion in near-fault regions in Lebanon, due to possible rupture scenarios on the northern part of the Yammouneh Fault for a Mw7 using a hybrid ground motion simulation technique presented herein. First, the source rupture is simulated up to 1 Hz, considering subshear and supershear rupture speed (see chapter 2 for the rupture simulation method), and different nucleation locations of the rupture. The ground motion is then computed up to 1 Hz, using a 1D medium representative for Lebanon (see chapter 5). The broadband ground motion up to 10 Hz is then generated using a stochastic empirical model that is calibrated to worldwide recordings of large earthquakes in the near-fault area in the contrary to classical hybrid approaches, in which the low frequency (i.e. deterministic) and high frequency (i.e. stochastic) ground motions are computed independently and next combined using low and high pass filters; this hybrid model takes into account the characteristic of the forward directivity pulses simulated in the first stage to generate a suite of broadband acceleration time histories. This hybrid model allows computing the acceleration time series for critical source scenarios reflecting the directivity effects for near-fault locality. Acceleration response spectra are next calculated and compared to a set of Ground Motion Prediction Equations (GMPEs) and to the design spectrum for Lebanon.

Introduction

Accumulated ground motion data have been providing us very important knowledge about rupture processes of earthquakes, propagation characteristics, ground motion amplification due to site effects, relation between ground motion and damage, and so on. Luckily, the occurrence of destructive large seismic events is less frequent than the occurrence of small and moderate earthquakes. On the other hand, the library of existing recordings only samples a small subset of possible earthquake scenarios of large earthquakes, making the understanding and the prediction of future strong ground motion uncertain. As an alternative, advances in the understanding of fault rupture processes, wave propagation and site response characterization, coupled with the tremendous growth in computational power and efficiency has made the prospect of large-scale ground motion time series generation for future earthquakes much more feasible. Therefore, when strong ground motion recordings are not available, which is the case mainly for large earthquakes at near fault localities or in low seismicity region, strong motion simulations may be used instead for low frequencies (up to ~ 1-3 Hz) (e.g. [START_REF] Mena | Pseudodynamic Source Characterization for Strike-Slip Faulting Including Stress Heterogeneity and Super-Shear Ruptures[END_REF][START_REF] Moschetti | 3D Ground-Motion Simulations of Mw 7 Earthquakes on the Salt Lake City Segment of the Wasatch Fault Zone: Variability of Long-Period (T≥1 s) Ground Motions and Sensitivity to Kinematic Rupture Parameters[END_REF][START_REF] Ramirez-Guzman | Ground-motion Simulations of 1811-1812 New Madrid Earthquakes, Central United States[END_REF]. At higher frequencies, seismological observation showed that source radiation and wave propagation effects tend to become stochastic [START_REF] Boore | Stochastic Simulation of High-Frequency Ground Motions Based on Seismological Models of the Radiated Spectra[END_REF], primarily reflecting our relative lack of knowledge about the details of these phenomena at higher frequencies.

Hybrid broadband simulation techniques have been developed in which the lowfrequency (LF) and the high-frequency (HF) motions are generated separately, and then superposed using a highcut filter for the former and a lowcut filter for the latter, to produce broadband synthetics for the entire frequency band of interest (~ 0.1-20 Hz). The low-frequency motions are deterministically generated by modeling the earthquake source process and the propagation of energy in the medium. The source and propagation modeling require a detailed velocity structure of the region, as well as the description of the fault rupture parameters. On the other hand, the HF motions generation uses the physics of wave scattering methods (e.g. [START_REF] Zeng | A Composite Source Model for Computing Realistic Synthetic Strong Ground Motions[END_REF][START_REF] Hartzell | Calculation of Broadband Time Histories of Ground Motion, Part II: Kinematic and Dynamic Modeling Using Theoretical Green's Functions and Comparison with the 1994 Northridge Earthquake[END_REF] or the stochastic approaches (e.g. [START_REF] Graves | Broadband Time History Simulation Using a Hybrid Approach[END_REF][START_REF] Liu | Prediction of Broadband Ground-Motion Time Histories; Hybrid Low/High-Frequency Method with Correlated Random Source Parameters[END_REF][START_REF] Frankel | A Constant Stress-Drop Model for Producing Broadband Synthetic Seismograms; Comparison with the next Generation Attenuation Relations[END_REF]. The stochastic simulations have the considerable advantage of being simple and versatile and requiring little advance information on the slip distribution or details of the Earth structure. In a stochastic simulation technique, the motions are treated as a random Gaussian signal (white noise) superimposed by the theoretical spectrum defined by simple seismological model of source and propagation filters, typically as a function of magnitude and distance [Acc (M 0 ; R; f) = Source (M 0 ; f) * Path (R; f) * Site (f)] (e.g. [START_REF] Hanks | The Character of High-Frequency Strong Ground Motion[END_REF][START_REF] Boore | Stochastic Simulation of High-Frequency Ground Motions Based on Seismological Models of the Radiated Spectra[END_REF][START_REF] Atkinson | Ground-Motion Relations for Eastern North America[END_REF][START_REF] Boore | Stochastic Simulation of High-Frequency Ground Motions Based on Seismological Models of the Radiated Spectra[END_REF].

However, unlike real seismic ground motions, the filtered white-noise process lacks nonstationarity in both the time and frequency domains [START_REF] Sabetta | Estimation of Response Spectra and Simulation of Nonstationary Earthquake Ground Motions[END_REF][START_REF] Pousse | Nonstationary Stochastic Simulation of Strong Ground Motion Time Histories Including Natural Variability: Application to the K-Net Japanese Database[END_REF][START_REF] Rezaeian | Stochastic Modeling and Simulation of Ground Motions for Performance-Based Earthquake Engineering[END_REF][START_REF] Yamamoto | Stochastic Model for Earthquake Ground Motion Using Wavelet Packets[END_REF]. The nonstationarity in the time domain refers to the variation of the intensity of the ground motion in time, which gradually increases from zero to achieve a nearly constant intensity, representing the strong-shaking phase of an earthquake, and then gradually decays back to zero. The nonstationarity in the frequency domain refers to the variation of the frequency content of the motion in time. Typically, high-frequency P-waves dominate the initial few seconds of the motion. These are followed by moderate-frequency S-waves, which dominate the strong-motion phase of the ground motion. Towards the end of the shaking, the ground motion is dominated by low-frequency surface waves. ( constructed and correlations between parameters of the two components are empirically determined, using a random-effects regression analysis method. This method reflects the weighing observations and the statistical dependence between the data within an earthquake cluster. For a given earthquake and site characteristics, the stochastic model reproduces in the synthetics the variability present in real ground motions, conserving the intensity, duration and frequency content.

Following the same procedure, [START_REF] Dabaghi | Stochastic Model for Simulation of Near-fault Ground Motions[END_REF]2018) presented a stochastic model, however this time, for near-fault ground motion, taking into account the near fault effects, mainly the forward directivity. [START_REF] Dabaghi | Stochastic Model for Simulation of Near-fault Ground Motions[END_REF]2018b) studied 441 near-fault ground motions from moderate to large earthquakes (5.5 < M < 8) recorded at sites with the closest distance to the fault rupture Rrup < 31 km from the NGA-West2 database, to obtain empirical observations of the model parameters (the Arias intensity, the frequency content, the duration of the recordings, the cumulative number of zero-level crossings and the cumulative number of positive minima and maxima). In this database, ~30% of the recorded ground-motion are "pulse-like" (see chapter 3 for pulse-like definition), ~15% are pulse-like and associated to strike-slip events [START_REF] Dabaghi | Stochastic Modeling and Simulation of Near-Fault Ground Motions for Performance-Based Earthquake Engineering[END_REF], and 6% are issued from events where supershear rupture occurred.

The objective of this study is to generate broadband synthetic seismograms in the near fault, that are consistent with the overall characteristics of strong ground motions expected to be observed in the near fault region in Lebanon for a Mw7 rupture on the Yammouneh fault. Therefore, a new technique for the hybrid approach is presented and used to simulate the strong ground motion in the near-fault. At frequencies below 1 Hz, the rupture is kinematically simulated using a statistical-based kinematic source model developed in chapter 2. The generated seismic waves are propagated in a 1D medium representative of the crustal structure in Lebanon (chapter 5). The ground motion is computed on the surface up to 1 Hz and the directivity pulse is extracted when it exists. The low frequency ground motion is complemented by high-frequency ground motion up to 10 Hz using a stochastic model fitted to near fault observations developed by (Dabaghi and Der Kiureghian 2018a), taking into account the directive pulse-like ground motions. The peak ground acceleration and response spectra mean value and variability are computed for a number of realistic rupture realizations and compared to the Lebanese seismic design value and to empirical Ground Motion Prediction Equations (GMPEs) for a Mw 7 earthquake.

Tectonic setting and seismic hazard

Lebanon is located in an active tectonic environment where the seismic hazard is considered moderate to high; it is crossed by the continental plate boundary-the Dead Sea Transform Fault. Within the Lebanese bend, the fault splits into four main branches that generated Mw larger than 7 earthquakes in the past, of which the Yammouneh fault constitutes the main fault branch since it bisects the length of Lebanon. The Yammouneh Fault bends eastward, is slipping at an estimated rate of 2-6 mm/year, and has been linked to many large earthquakes in the past, for e.g. the historical 1202 earthquake of magnitude Ms 7.6 [START_REF] Plassard | Seismicité Du Liban[END_REF][START_REF] Ambraseys | Faulting Associated with Historical and Recent Earthquakes in the Eastern Mediterranean Region[END_REF]Ben-Menahem 1991b;[START_REF] Ellenblum | Crusader Castle Torn Apart by Earthquake at Dawn[END_REF][START_REF] Daëron | Constraints on the Post∼ 25-Ka Slip Rate of the Yammoûneh Fault (Lebanon) Using in Situ Cosmogenic 36Cl Dating of Offset Limestone-Clast Fans[END_REF][START_REF] Daeron | 12,000-Year-Long Record of 10 to 13 Paleoearthquakes on the Yammouneh Fault, Levant Fault System, Lebanon[END_REF][START_REF] Vergnolle | GPS Measurements of Crustal Deformation in Lebanon: Implication for Current Kinematics of the Sinaï Plate[END_REF]). The mean recurrence period for large earthquakes along of the Yammouneh Fault is 990-1260 yr [START_REF] Daëron | 12,000-Year-Long Record of 10 to 13 Paleoearthquakes on the Yammoûneh Fault, Levant Fault System, Lebanon[END_REF]. Present day, in Lebanon, instrumentally recorded seismicity of M < 5 is generally sparse within the Lebanese restraining bend. Strong motion has never been recorded in Lebanon till now due to the presently infrequent large-magnitude seismicity.

In order to assess the seismic hazard in Lebanon, scientists often use the ground motion attenuation relationships GMPEs calibrated by a large set of recorded earthquakes worldwide (e.g. [START_REF] Huijer | Re-Evaluation and Updating of the Seismic Hazard of Lebanon[END_REF]. One of the first attempt to simulate the ground motion for Lebanon due to a rupture on the Yammouneh Fault was the work of (Brax, Causse, and Bard 2016b). They simulated a Mw6.5 (L ~ 23 km) rupture using empirical Green's function's technique developed by [START_REF] Causse | New Approach for Coupling K-2 and Empirical Green's Functions: Application to the Blind Prediction of Broad-Band Ground Motion in the Grenoble[END_REF], and predicted a mean pga of 0.08 g seismic ground at Bhannes, located around 23 km from the rupture area. However, their simulations were limited up to Mw6.5, considering the limits of the method in the near field of extended sources. Nevertheless, because of its potential to generate large earthquakes (larger than Mw7), and because of its length as it crosses the Lebanese territory from South to North; consequently, the Yammouneh Fault poses a seismic threat to the population centered in its vicinity. The seismic hazard needs to be addressed using simulation for the near fault ground motion.

Fault rupture segment and target stations

According to [START_REF] Manighetti | Evidence for Self-Similar, Triangular Slip Distributions on Earthquakes; Implications for Earthquake and Fault Mechanics[END_REF], most earthquakes are observed to nucleate where two strongly oblique (commonly roughly perpendicular) faults intersect. Then they propagate and grow unilaterally until they touch another structural discontinuity. In this study, the considered 70 km rupture segment is located in the northern part of the Yammouneh Fault, in between two intersections with two other faults (the Mount Lebanon thrust in the north of the segment, and the Mid-Bekaa Fault in the south of the rupture segment (Figure 5-2). This segment has been seismically calm and records a very few small earthquakes (GRAL Network-CNRS Lebanon). Next, the location of the site to compute the near-fault ground motion is chosen.

Most of the population in Lebanon resides on the coast (Figure 6-1); more than half of the population resides in Beirut and its suburbs, Tripoli and its suburbs, and Zahle.

The imbalance in population distribution on Lebanese territory, mainly in Beirut, is up to the fact that in addition to being the capital, it has the political, administrative and economic centralizations. Besides from being highly populated, those concentrations face the high risk of near-fault earthquake ruptures. Therefore, we synthesize near-fault ground motions for stations S1 and S2 located at rupture distances Rrup of 5 km and for stations S3 and S4 at 25 km, located at mid-length of the rupture (S1 and S3) and at 10 km beyond the rupture length (S2 and S4), shown in Figure 6-1. Given the potential rupture position, it is then essential to analyze the effect of the rupture directivity toward the area of Beirut.

Ground-motion simulation methodology

Hybrid broadband ground-motion is generated using a combination of two different approaches at different ranges of frequencies, less than 1 Hz and larger than 1 Hz, to produce broadband synthetics for the entire frequency of interest (e.g. [START_REF] Graves | Broadband Time History Simulation Using a Hybrid Approach[END_REF][START_REF] Liu | Prediction of Broadband Ground-Motion Time Histories; Hybrid Low/High-Frequency Method with Correlated Random Source Parameters[END_REF][START_REF] Frankel | A Constant Stress-Drop Model for Producing Broadband Synthetic Seismograms; Comparison with the next Generation Attenuation Relations[END_REF]. The low-frequency ground motion is first generated via simulation of the source rupture and propagated to the surface via the crustal medium. The directivity pulse -when it exists-is extracted from the low frequency ground motion. The broadband ground motion is then generated using a calibrated stochastic model for given earthquake, at a given site, knowing the characteristics of the pulse or the non-pulse low frequency ground motion. The procedure is detailed in the following and is summarized in the flowchart presented in Figure 6-2. 6.5.1 Low-frequency ground motion (f ≤ 1 Hz)

The low-frequency ground motion (GM1Hz (t) in Figure 6-2) at the surface is generated via 2 steps.

Step 1 consists of simulating the source rupture on a vertical strike slip fault for critical rupture scenarios, and step 2 consists of propagating the seismic waves into the medium representative of Lebanon.

Pseudo-dynamic source model

One approach to generate the source rupture is the pseudo-dynamic source modeling developed in this work (see chapter 2). It consists of a priori prescribing the displacement discontinuity across the fault surface. The rupture starts from the hypocenter and expands over the fault plane with a rupture speed (Vr). Each point on the fault slips as it is reached by the rupture front and is characterized by the source function, also called the slip velocity function (SVF). SVF describes the evolution of the slip with time and is defined by the following parameters: the final slip value (D), the peak slip velocity (psv) and the rise time (Trise) that represents the time needed to reach the final slip. The statistical properties of the source parameters have been introduced in details in chapter 2 (section 0). For generating suites of rupture realizations for our ground motion predictions in Lebanon, these statistical parameters presented in Table 6-1, are constrained by the analysis of a database of dynamic models [START_REF] Song | Importance of 1-Point Statistics in Earthquake Source Modelling for Ground Motion Simulation[END_REF]. Using the pseudo-dynamic source model, we generate rupture models for a Mw7. The rupture length L = 70 km and width W = 14 km are derived from the Mw-L scaling relationship by [START_REF] Papazachos | Global Relations between Seismic Fault Parameters and Moment Magnitude of Earthquakes[END_REF]). The mean value of the slip 𝜇 𝐷 is then defined by: 𝜇 𝐷 = 𝑀 0

𝐺 𝐿 𝑊

, where 𝑀 0 is the seismic moment and 𝐺 is the shear modulus. We make sure that the maximum slip does not exceed the ceiling defined by [START_REF] Mcgarr | Maximum Slip in Earthquake Fault Zones, Apparent Stress, and Stick-Slip Friction[END_REF] as a function of magnitude. The mean value of the psv is chosen from the database of spontaneous dynamic rupture simulations developed by [START_REF] Song | Pseudo-Dynamic Source Modelling with 1-Point and 2-Point Statistics of Earthquake Source Parameters[END_REF]. The slip duration Trise is calculated as a function of the psv and the D, for a regularized Yoffe slip-rate function [START_REF] Tinti | A Kinematic Source-Time Function Compatible with Earthquake Dynamics[END_REF]. Trise is allowed to vary between 0.1 and 5 s. The fault area is embedded in at 0.5 km below the surface. The spatial distribution of final slip and rupture speed are tapered so as to avoid stress singularities at the fault boundaries. Therefore, the values of the source parameters decrease as they reach 20% of the fault dimension to reach zero at each side, with a quarter circular taper. We fix the hypocenter depth to 80% along the dip [START_REF] Mai | Hypocenter Locations in Finite-Source Rupture Models[END_REF]. Among the many source parameters, the sensitivity analysis performed in chapter 2 shows that the rupture speed and the hypocenter location are highly affecting the surface ground motion in the near-fault [START_REF] Somerville | Modification of Empirical Strong Ground Motion Attenuation Relations to Include the Amplitude and Duration Effects of Rupture Directivity[END_REF][START_REF] Aagaard | Characterization of Near-Source Ground Motions with Earthquake Simulations[END_REF][START_REF] Schmedes | Near-Source Ground Motion along Strike-Slip Faults: Insights into Magnitude Saturation of PGV and PGA[END_REF][START_REF] Ripperger | Variability of Near-Field Ground Motion from Dynamic Earthquake Rupture Simulations[END_REF][START_REF] Fayjaloun | CAPTURING THE SOURCE PARAMETERS CONTROLING THE HIGH-FREQUENCY GROUND MOTION[END_REF]. Both parameters affect the rupture directivity and it is crucial to take into account the range of their variability in our study for Lebanon. In addition, we have seen that for a given set of statistical source parameter, various rupture realizations can lead to highly different ground motions values. We then generate 50 rupture realizations for each set of source parameters.

Table 6-1: The source parameters: D, Vr, ax, az, and psv, stand for the slip, rupture speed, spatial correlation lengths along the strike and along the dip directions, and peak slip velocity, respectively. 𝝁 represents the mean value, 𝝈 is the standard deviation and 𝝆 is the coefficient of correlation. μD (cm) σD/μD μVr (km/s) σVr/μVr ρD-Vr ax (km) az (km) μpsv (cm/s) σpsv ρD-psv 120 0.5 3.12*0.85 3.12*1.50 0.2 0 16 5 160 80 0.8

Rupture speed

Seismological studies report that the rupture front typically propagates at ~80% of the shear-wave speed for crustal earthquakes [START_REF] Heaton | Evidence for and Implications of Self-Healing Pulses of Slip in Earthquake Rupture[END_REF][START_REF] Mai | SRCMOD: An Online Database of Finite-fault Rupture Models[END_REF]. However, the rupture speed may exceed the shear wave speed, as shown by theoretical and observational studies. The first earthquake for which supershear wave rupture speed was inferred was the 1979 Imperial Valley, California Earthquake, studied by [START_REF] Archuleta | A Faulting Model for the 1979 Imperial Valley Earthquake[END_REF][START_REF] Spudich | Direct Observation of Rupture Propagation during the 1979 Imperial Valley Earthquake Using a Short Baseline Accelerometer Array[END_REF]. Then, in the late 1990 to early 2000s, a few additional earthquakes with supershear wave rupture speeds were reported. Supershear ruptures where observed during the 1999 Izmit and the 1999 Duzce earthquakes on the strike-slip North Anatolian Fault with Mw = 7.6 and 7.2 respectively [START_REF] Bouchon | How Fast Is Rupture during an Earthquake? New Insights from the 1999 Turkey Earthquakes[END_REF], the 2001 Kunlun earthquake on the strike-slip Kunlun fault of Mw = 7.8 [START_REF] Walker | Illuminating the Near-sonic Rupture Velocities of the Intracontinental Kokoxili Mw 7.8 and Denali Fault Mw 7.9 Strike-slip Earthquakes with Global P Wave Back Projection Imaging[END_REF][START_REF] Vallée | Observation of Far-field Mach Waves Generated by the 2001 Kokoxili Supershear Earthquake[END_REF], the 2002 Denali earthquake on the strike-slip Denali fault of Mw = 7.9 [START_REF] Ellsworth | Near-Field Ground Motion of the 2002 Denali Fault, Alaska, Earthquake Recorded at Pump Station 10[END_REF][START_REF] Aagaard | Near-Source Ground Motions from Simulations of Sustained Intersonic and Supersonic Fault Ruptures[END_REF][START_REF] Dunham | Evidence for a Supershear Transient during the 2002 Denali Fault Earthquake[END_REF]), the 2010 Qinghai-China earthquake on the strike-slip Yushu fault of Mw = 6.9 (Wang, Mori, and Uchide 2012) and the 2013 Mw = 7.5 Craig-Alaska earthquake [START_REF] Yue | Supershear Rupture of the 5 January 2013 Craig, Alaska (Mw 7.5) Earthquake[END_REF]. Supershear rupture speed was also inferred for the 1906 San Francisco earthquake on the strike-slip San Andreas Fault of Mw = 7.8 [START_REF] Song | Evidence for Supershear Rupture during the 1906 San Francisco Earthquake[END_REF]. The supershear was observed to occur on segments where the fault was straight, long, narrow, and with a simple geometry; a lack of aftershocks and a low background seismicity of these segments were also noticed [START_REF] Bouchon | Faulting Characteristics of Supershear Earthquakes[END_REF]. [START_REF] Ansal | Perspectives on European Earthquake Engineering and Seismology[END_REF] identifies the Dead Sea Fault as a major active-slip fault system with long straight portions capable of sustained supershear rupture speeds. For our simulations, we consider two different cases: subshear (Vr = 0.85 * Vs) and supershear (Vr = 1.5 * Vs) rupture speeds.

Hypocenter location

The ground motion for a given moment magnitude is strongly affected by the hypocenter position [START_REF] Aagaard | Characterization of Near-Source Ground Motions with Earthquake Simulations[END_REF][START_REF] Somerville | Modification of Empirical Strong Ground Motion Attenuation Relations to Include the Amplitude and Duration Effects of Rupture Directivity[END_REF][START_REF] Ripperger | Variability of Near-Field Ground Motion from Dynamic Earthquake Rupture Simulations[END_REF]. 80% of shallow strong ruptures are predominantly unilateral [START_REF] Mcguire | Predominance of Unilateral Rupture for a Global Catalog of Large Earthquakes[END_REF], and the hypocenter is located at 20-30% of the total length of the rupture from the fault lateral edge [START_REF] Manighetti | Evidence for Self-Similar, Triangular Slip Distributions on Earthquakes; Implications for Earthquake and Fault Mechanics[END_REF]. We assume 3 possible positions for the hypocenter: at 20% ('L'), 50% ('C') and 80% ('R') respectively from the length of the rupture, represented by the 3 stars in Figure 6-1, therefore 75% of the ruptures are unilateral. The 3 possible positions for the hypocenter, and the 2 values for the rupture speed, produce 6 different source scenarios that are annotated as the following: L, C, R, LS, CS and RS, (S) standing for supershear. For each scenario, 50 realizations are run to take into account the spatial variability of the source parameters. Then, we propagate the seismic waves from the source to the surface. According to chapter 5, the propagation medium (resolved up to 1 km depth) can be sufficiently assumed to be 1D since the crust under Lebanon for at least the first 15 km is plane stratified. For the first km, we perform surface waves inversion combining Rayleigh wave group velocities (chapter 5) for frequencies lower than 1 Hz and Rayleigh and Love waves phase velocities for higher frequencies inferred from small aperture seismic ambient noise arrays [START_REF] Cornou | Shear-Wave Velocity Structure and Correlation with n-Spt Values in Different Geological Formations in Beirut[END_REF]Brax, Causse, and Bard 2016b) In a perfect elastic medium, the total energy of the wavefield is conserved. However, in reality, the medium is not perfectly elastic; the medium has a seismic intrinsic attenuation property characterized by a quality factor Qi, which is defined as the ratio of wave energy to the energy dissipated per cycle of oscillation. We use studies in the surrounding region to estimate Qi-s for the shear waves in Lebanon. Qis = 100 [START_REF] Meirova | Seismic Wave Attenuation in Israel Region Estimated from the Multiple Lapse Time Window Analysis and S-Wave Coda Decay Rate[END_REF] for the shear waves in Galilee Lebanon region. The low value of Q reflects the strong attenuation that is mainly due to the region being tectonically active to the south of Lebanon [START_REF] Meirova | Seismic Wave Attenuation in Israel Region Estimated from the Multiple Lapse Time Window Analysis and S-Wave Coda Decay Rate[END_REF]. Following [START_REF] Coulomb | Mantle and Core in Planetary Physics[END_REF], Qi-P = 2.25 Qi-s. We used this relationship to estimate the quality factor of the attenuation factor for the compressional waves. We synthesize near-fault ground motions in a 1D layered medium for stations located at rupture distances Rrup of 5 km (stations S1 and S2) and 25 km (stations S3 and S4) using the Representation theorem (see chapter 1). The broadband near-fault ground motion (GM10Hz (t) in Figure 6-2) is computed using a stochastic model that reproduces the statistical parameters of the near-fault recordings for a given earthquake and a specific site, knowing the pulse characteristic issued from the directivity effects. The procedures are detailed hereafter.

The simulated low-frequency horizontal ground motions are classified as pulse-like or non-pulse-like following [START_REF] Shahi | An Empirically Calibrated Framework for Including the Effects of Near-Fault Directivity in Probabilistic Seismic Hazard Analysis[END_REF], based on the significance of the extracted pulse relative to the original ground motion. The pulse is obtained using wavelet analysis; it corresponds to the wavelet associated with the highest coefficient (see [START_REF] Baker | Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis[END_REF]) for details). The pulse-like behavior is then determined by a pulse indicator (PI) defined by [START_REF] Baker | Quantitative Classification of Near-Fault Ground Motions Using Wavelet Analysis[END_REF]) and computed for each velocity time series. The PI depends on the PGV ratio (ratio of the PGV of the residual to the PGV of the original ground motion) and the energy ratio (ratio of the energy of the residual to the energy of the original ground motion, where energy is computed as the cumulative squared velocity of the record). As well, in order to be classified as pulse-like due to the forward directivity, the pulse should arrive early in the ground motion record. The pulse period is therefore the period associated with the maximum Fourier amplitude spectra of the wavelet having the largest wavelet coefficient. An example of pulse detection is shown in Figure 6-4 for 2 different source scenarios shown in Figure 6-5 and Figure 6-6. The 5 pulse parameters are then extracted: the pulse amplitude Vp, the pulse period Tp, parameter ϒ characterizing the number of oscillations within the pulse, the phase angle ν, and the time of the peak of the modulating envelope tmax,p. Consequently, if the low-frequency ground motion is classified as pulse-like, 50 stochastic two-components broadband ground motion are generated conditioned by the characteristics of the pulse extracted; if the low-frequency ground motion is classified as non pulse-like, 50 stochastic two-components broadband ground motion are generated conditioned by the fact that they don't have a pulse. The stochastic broadband time series are generated using a stochastic model with evolving temporal and spectral characteristics, calibrated with near-fault recordings for the given earthquake source and site characteristics (type of faulting F, earthquake magnitude Mw, depth to top of rupture plane ZTor, source-to-site distance Rrup, site characteristics Vs30, and directivity parameters sord and ΘorΦ), developed by (Dabaghi and Der Kiureghian 2018a), but also knowing the extracted properties of the pulse (Vp, Tp, ϒ, ν, and tmax,p) when it exists. Therefore, the stochastic broadband ground motion is generated, based on the properties of the conditional simulation and the partitioning of the mean vector and the covariance matrix [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF]. The mathematical method is detailed in [START_REF] Rezaeian | Stochastic Modeling and Simulation of Ground Motions for Performance-Based Earthquake Engineering[END_REF]). An example of source subshear scenario nucleating at hypocenter 'L' generating a pulse for station S3 is illustrated in Figure 6-5, with the corresponding pulse-like ground motion (time series represented in red) for frequencies up to 1 Hz. Two aleatory broadband ground motion are also shown in grey. Figure 6-6 illustrates an example for supershear rupture. Figure 6-7 illustrates a source subshear scenario nucleating at hypocenter 'L', that is at a fault edge, for a non-pulse low-frequency ground motion at station S3 (shown in grey), and the corresponding broadband ground motion.

In this work, we assume it is possible to use the stochastic model for subshear ruptures, as well as for supershear models, although the stochastic model is Subshear generated ground motion Subshear detected pulse Supershear generated ground motion Supershear detected pulse essentially calibrated with events characterized by subshear ruptures. Indeed, [START_REF] Bizzarri | Coherence of Mach Fronts during Heterogeneous Supershear Earthquake Rupture Propagation: Simulations and Comparison with Observations[END_REF] studied ground motions issued from supershear ruptures, and found that there is no average elevation of ground motion acceleration relative to what ground motion prediction equations (GMPEs) could predict. [START_REF] Vyas | Mach Wave Properties in the Presence of Source and Medium Heterogeneity[END_REF]) justify the ground shaking amplitudes for supershear ruptures not being elevated due to wavefield scattering and rupture complexity.

The amplitude, frequency content and duration of the motion give engineers valuable information related to the extent of damage of structures. A convenient measure of the response of the structures can be obtained by evaluating the response of linear elastic 1-D systems using Duhamel's integration method, as function of the damping ratio ζ=5%, and the natural spectral period. Finally, the PGA and the response spectra at different spectral periods are computed using an orientation-independent measure proposed by [START_REF] Boore | Orientation-Independent Measures of Ground Motion[END_REF]) (GMRoTD50). This measure comprises a rotation of the two orthogonal components from 1 to 90, and evaluates the peak ground motion from the geometric mean of the rotated time series.

Simulation results

The plots in [Figure 6-8, Figure 6-9, Figure 6-10, Figure 6-11] corresponding to the different stations S1 through S4, show the values of the pulse periods for all the scenarios, as well as the values of the corresponding pga of GM1Hz generated using the low-frequency deterministic simulation, and the pga of GM10Hz generated using the hybrid stochastic broadband model up to 10 Hz. There is a positive correlation between the pga of GM1Hz and GM10Hz, with a factor larger than 0.5 for stations S1, S2 and S3, and about 0.3 for S4, where less pulses are recorded. This positive correlation reflects the fact that a pulse-like motion with a large pulse amplitude (associated to pga1Hz) tends to have high arias intensity (associated to pga10Hz), according to (Dabaghi and Der Kiureghian 2018b) who found a correlation coefficient of 0.4 between Vp and Ia parameters.

Peak ground acceleration at different stations

The median value of the PGA of low frequency ground motions at S1, S3 and S4 increases, and the variability decreases, when comparing the subshear and the supershear scenarios (Figure 6-12); however, the median value decreases at S2 when the rupture speed regime changes from subshear to supershear to reach a minimum value of 0.03 g for all three hypocenter locations. By adding the HF ground motion using the stochastic model calibrated by the near-fault recordings, the impact of the rupture speed regime on the median value is also observed for the broadband ground motion. This is because the high-frequency stochastic part is empirically calibrated and conditioned by the deterministic pulse features. Considering all the scenarios, and looking at stations located at 5 km from the rupture, the broadband ground motion at S1 situated at mid-length of the rupture has PGA median value of 0.43 g, with values around 0.28 g and 0.67 g considering one standard deviation (Figure 6-8). S2 positioned 10 km beyond the rupture length has PGA median value of 0.25 g, between 0.17 g and 0.35 g, considering the standard deviation (Figure 6789). Looking at stations located at 25 km from the rupture, the ground motion at S3 situated at mid-length from the rupture has PGA median value of 0.26 g, with values around 0.16 g and 0.35 g, considering one standard deviation (Figure 6-10). S4 located 10 km beyond the rupture length has PGA median value of 0.13 g, between 0.1 g and 0.2 g, considering one standard deviation . At 5 km, ASK14 and BSSA14 predicted PGA of ~0.35 g, however CB14 predicts a higher value of 0.67 g. At 25 km, ASK14, BSSA14 and CB14 predicted PGA of ~0.3, 0.4 and 0.47 g, respectively, all showing larger values than our simulations do for S3 and S4. GMPEs predict PGA varying between 0.2 and 1.2 g for stations at 5 km from the fault, and between 0.2 and 0.8 g for stations at 25 km from the fault. The results are summarized in Table 6-3. The pga values of the broadband ground motion issued from our hybrid model are lower than predicted by GMPEs at stations S2, S3 and S4. However, the difference of the pga values issued from the two methods is not shocking. Even though both the stochastic model developed by (Dabaghi and Der Kiureghian 2018b) and the GMPEs are calibrated using the NGA-West2 Database, however, the two different methods apply different criteria for the selection of ground motions to be used in their calibration. (Dabaghi and Der Kiureghian 2018b) used the ground motions issued from moderate to large earthquakes and recorded in the near-fault, and therefore used a subset (2% only) from the NGA-West2 database. GMPEs however used more than 70% of the database ground motion recordings (e.g. [START_REF] Gregor | Comparison of NGA-West2 GMPEs[END_REF]Campbell and Bozorgnia 2014), containing near-fault but mainly far-fault data.

The cities of Bcharre (at 8 km from the Yammouneh Fault) and Harissa (at 10 km) are located between S1 and S3, and therefore would be subjected to an expected PGA value of 0.3-0.4 g. Zahle (at 5 km) and Hammana (at 6 km) are comparable to S2, and would have a PGA value of 0.25 g. Bhamdoun (at 10 km), Hermel (at 10 km), Aley (at 15 km), Halba (at 20 km), Baabda (at 22 km), Mansourieh (at 22 km), Tripoli and Beirut… all are near-fault station between S2 and S4, and would be subjected to PGA values of 0.13-0.25 g as median values.

Among the scenarios tested, the worst rupture scenarios can also lead to more disastrous results: S1 and S2 can experience pga up to 1.3 g. S3 and S4 can experience respectively 0.6 g and 0.43 g ground shaking pga. Note that, in Lebanon, the design acceleration, is 0.25 g (decree 7964), except for the engineering projects located in Beirut Central District (BCD) falling under the responsibility of Solidere (Lebanese Company for the Development and Reconstruction of BCD) are designed for 0.3 g [START_REF] Brax | Aléa et Microzonage Sismiques à Beyrouth[END_REF]).

Response spectra

The response spectra are computed for all the 15,000 simulations for spectral periods between 0.1 and 10 s for the 4 different stations. Their median plus/minus standard deviation are plotted in Figure 6-13 and compared to GMPEs response spectra computed at distance 5 km and 25 km respectively. In general, our simulated response spectra fit within the GMPEs. However, our simulations show an amplification of the response spectra at low frequencies caused by directivity pulses in the near fault [START_REF] Spudich | Final Report of the NGA-West2 Directivity Working Group[END_REF]. At S1, the amplification is centered on ~3 s, and around 4 s for S2, S3 and S4, corresponding to the pulse period of 3.3 s for a Mw7 earthquake [START_REF] Shahi | An Empirically Calibrated Framework for Including the Effects of Near-Fault Directivity in Probabilistic Seismic Hazard Analysis[END_REF]. The GMPEs spectra, on the other hand, do not reflect the directivity effects.

Following the American code UBC97, the seismic design response spectra of structures on hard rock (Vs30 = 800 m/s) in Beirut and other regions in Lebanon are computed and also plotted in Figure 6-14, and compared to the resultant simulations of the Mw7 rupture on the northern part of the Yammouneh Fault. On the rock site, the seismic design in Lebanon is underestimating the response spectra for S1. It perfectly fits the other stations in terms of median and standard deviation. Nonetheless, there is a large variability of the response spectra due to the different scenarios of the rupture, and to the aleatory behavior of the ground motion reflecting the natural variability in the source and the propagation process. Therefore, the design spectra in Lebanon would be underestimating the ground motions, especially at periods close to the pulse periods.

R

Sub-shear Super-shear

Conclusion

The GMPEs result from empirical analysis of collected ground motion records worldwide for various conditions. The peak ground acceleration and the acceleration spectrum can be computed using GMPEs based on the earthquake magnitude, distance and soil condition that is expected to occur. These conventional response spectra are frequently used in engineering practice, and in particular for the standard design spectrum of Lebanon, however, these GMPEs are not ready yet to reflect the near fault directivity effects. On the other hand, considering the earthquake rupture scenario for the area is a necessity, especially in the near-fault vicinity. A new hybrid approach to compute a broadband nearfault ground motion is presented. The simulation of the ground motion at low frequency (up to 1 Hz) takes into account the physics of the rupture and of the propagation of the seismic waves. Accordingly, the broadband ground motion is generated stochastically (up to 10 Hz) while respecting the characteristics of the low-frequency ground motion. This correlation between the low-frequency and the high-frequency ground motion is guaranteed by the stochastic model that is empirically calibrated by recordings of past events. The major advantage is that it is a combination of seismological approach (pseudodynamic source model) reflecting the physics of the rupture and accounting for the pulselike behavior of near-fault ground motion due to the directivity effect, on the one hand, and stochastic approach, which allows to covering a broadband frequency range.

In this work, subshear and supershear rupture scenarios on the Yammouneh Fault are generated to estimate the ground motion near the fault rupture for Mw7, on a segment of 70 km in the northern part of the Yammouneh Fault. Sites located at 5 km from the fault have a median pga values of 0.25 -0.43 g (depending on the location of the site with respect to the rupture); a pga up to 0.67 g is very probable. Farther sites located at 25 km from the fault (e.g. Beirut) have median pga values of 0.13 -0.26 g; a pga up to 0.35 g is very probable. The response spectra computed at the different sites showed an amplification of the spectral acceleration at periods around the period of the pulse (between 2 s and 5 s) issued from the directivity. This amplification severely affects highrise buildings with more than 30 stories [START_REF] Salameh | Seismic Response of Beirut (Lebanon) Buildings: Instrumental Results from Ambient Vibrations[END_REF].

This study can be improved in the future by taking into account the 3D Vs propagation crust medium under Lebanon, as well as the Vs30 spatial variability of the soil at the surface. Ground motion at stations covering the whole region could be computed to better assess the seismic hazard in Lebanon, along with incorporating other ruptures on the different faults in Lebanon, especially ruptures up to a Mw = 7.5 that are highly probable in Lebanon, and for which larger pga and response spectra are expected.

Conclusions and Perspectives

The earthquake engineering community is interested in the damaging effects of near-fault strong ground motions on different structural systems, and there is an increasing effort in the seismology community to incorporate the source rupture and the near-fault effects (extended source effects like the rupture directivity) within the probabilistic seismic hazard framework. The link between the seismological parameters (source properties, propagation medium and site characteristics) and the ground motion (peak ground velocity and acceleration and response spectra) remains a hot subject of seismological studies and the drive for research projects.

During my PhD, I investigated the role of some major seismological source parameters (rupture speed, stress drop, asperity size, correlation coefficient between source parameters, hypocenter location) at different sites to understand and to quantify the variations in ground-motion demands (PGA and PGV) (chapter 2). I concluded that PGA is mainly generated by abrupt changes of the rupture propagation, that is, stopping phases at the fault boundaries or strong heterogeneities of rupture speed along the rupture. PGA is mostly controlled by the location of the hypocenter, the average rupture speed, and to a lesser extent by the average stress drop (in the far-field) and the amplitude of the rupture speed heterogeneities. Interestingly, the correlation between the source parameters and the spatial correlation length (characteristic size of source heterogeneities) do not significantly affect the average PGA. Nevertheless, they play a significant role in the ground motion variability. This sensitivity analysis merits to be performed for a complete coverage of stations around the rupture fault length. In addition, the sensitivity analysis could be extended to the variability of the peak values.

Despite continuous expansion of the database of recorded earthquake ground motions, the empirical GMPEs are still somewhat unsatisfying for 1) estimating the strong ground motion in the near fault due to insufficient such recordings and 2) reproducing the directivity effects that generate a pulse with large amplitude, mainly observed in the near fault zone. When the rupture propagates towards a given site, the forward directivity effect is reflected by a peak in the response spectrum near the period of the pulse. In the literature, the period of the pulse was essentially estimated as function of the magnitude. In chapter 3, I considered the spatial variability of the pulse period recorded in the past events and introduced a simple equation to better predict its value. The numerical ground motion simulations based on more realistic homogeneous ruptures (chapter 2) could provide a complement to capture additional physical features contributing to the pulse generation.

In order to extend ground motion simulations to a broader frequency range useful for the earthquake engineering community, I presented in chapter 6, a new hybrid approach to compute a near-fault ground motion up to a frequency of 10 Hz, but still reflecting the directivity effects on the response spectrum at low frequency. The low frequency ground motion (up to 1 Hz) takes into account the physics of the rupture (as described in chapter 2) and of the propagation of the seismic waves in a 1D elastic medium. The broadband ground motion is then generated stochastically (up to 10 Hz) conditioned by the characteristics of the low-frequency ground motion. This correlation between the low-frequency and the high-frequency ground motion is guaranteed by the stochastic model that is empirically calibrated by near fault strong motion recordings of past events. As such, the method incorporates some important physical processes of the rupture (pseudo-dynamic source model), accounts for the pulse-like behavior of near-fault ground motion due to the directivity effect (reflected by an amplification of the response spectra at periods close to the period of the pulse) and allows to cover a broadband frequency range due to the use of stochastic approaches. This model provides suites of ground motion time histories with empirically calibrated duration, peak values and frequency content, and could be used to study the nonlinear response of structures.

I cannot conclude without commenting on Lebanon seismic hazard assessment. There is still work to be done, there is no doubt. Near fault strong ground motion simulations are indispensable for Lebanon seismic hazard. This small country is crossed by a large fault that makes all its cities in the near fault zone and prone to near fault effects. Lebanon is classified as moderate to high seismicity, nevertheless, no strong motion have been recorded. A direct application of the above-mentioned hybrid model was applied in chapter 6 to estimate the ground motion (pga and response spectra) for sites near the Yammouneh Fault, for a hypothetical Mw7 rupture. The 1D velocity structure was derived from ambient noise analysis (chapter 5). These simulations would help to better assess the seismic hazard in Lebanon, once combined with the known seismicity background.

Nevertheless, earthquakes larger than to Mw 7.2 could occur on the YF [START_REF] Sadek | Updated Seismic Hazard for Lebanon and Implications on Micro-Zonation of the Greater Beirut Area[END_REF][START_REF] Lefevre | Slip Deficit and Temporal Clustering along the Dead Sea Fault from Paleoseismological Investigations[END_REF]. The observation of speleothems in Jeita and Kanaan caves (~20 km from the YF), coupled with a statistical approach developed by [START_REF] Lacave | What Can Be Concluded about Seismic History from Broken and Unbroken Speleothems[END_REF], confirmed that broken speleothems indicate strong earthquakes with acceleration between 0.2 g and 0.6 g occurring in the past (Libris Report 2012). However, these accelerations are not necessarily attributed to the Yammouneh Fault ruptures. Ruptures up to a M7.5 should then be simulated on the Yammouneh Fault, as well as on other faults in Lebanon, mainly the Mount Lebanon thrust Fault that generated destructive M7.5 and large tsunamis in the past. And since the seismic hazard evaluation of Lebanon is mainly influenced by the seismic activity within the country but also from the neighboring regions, there is an interest to consider the rupture of the Levant Fault in different countries, especially since the frontiers between these regions are generally defined politically rather than on the basis of any seismo-tectonic boundaries [START_REF] Sadek | Updated Seismic Hazard for Lebanon and Implications on Micro-Zonation of the Greater Beirut Area[END_REF].

In chapter 5, I presented a 3D shear wave tomography for Lebanon between 1 km and 40 km using the technique of ambient noise cross-correlation.

The 3D VS pseudo model would help to, first, better comprehend the geodynamics in the region, and second, to better model wave propagations for seismic hazard assessment. Important interfaces are depicted at 2-3 km, 8-15 km, and larger than 20 km. The traces of the Yammouneh fault and the coastal basins are depicted at subsurface up to ~4 km. Homogeneous structures throughout the country are outlined deeper up to ~ 20 km. The deepest interface is probably the Moho, dipping from the west to the east between 20 and 40 km. A 1D velocity structure was deduced and used to perform the hybrid near-fault ground motion simulations for a Mw7 rupture on the Yammouneh fault (chapter 6). Nonetheless, the pseudo-model could be improved in the future. A higher number of receivers would significantly scale down the level of uncertainty of the inverted model. Small-scale arrays of receivers would allow resolving shallower depths. Body waves tomography would allow depicting more precisely the Moho depth and to image better the interaction at depth of the various faults.

Next to last, continuous measurements of ambient seismic noise over few years could be used for seismic monitoring. Actually, seismic velocity changes are related to co-seismic damage in the shallow layers and to deep co-seismic stress change and postseismic stress relaxation within the fault zone. The temporal evolution of the crust (variations of seismic velocities within the crust) can be tracked by computing cross-correlation functions at different dates for the same receiver pair and measuring the changes between the correlation functions [START_REF] Brenguier | Postseismic Relaxation along the San Andreas Fault at Parkfield from Continuous Seismological Observations[END_REF]. Seismic monitoring could help predicting the earthquake events in Lebanon.

Lastly, interdisciplinary collaborations between seismologists, geologists, earth scientists, speleologists, and historians of catastrophes and civil engineers would help better assess the natural geo-hazards in Lebanon, and mainly the seismic hazards.

  of source rupture and near-fault ground motion 1.1.1 General introduction about earthquakes 1.1.1.1 Why do we have earthquakes?

Figure 1 - 1 :

 11 Figure 1-1: Stress evolution and energy balance during faulting process at a given fault point. Shaded area represent the amount of radiated seismic energy (𝚬 𝒔 ), fracture energy (𝚬 𝒈 ) and friction energy (𝚬 𝒇 ), after (Ide 2003). Linear slip weakening friction law. 1.1.1.2 Can we be resilient to earthquakes?

Figure 1 - 2 :

 12 Figure 1-2: An illustration of the earthquake from the source rupture to the ground motion and the response of the structure.

Figure 1 - 3 :

 13 Figure 1-3: Idealized sinusoidal pulses in the far field, after (Kalkan and Kunnath 2006).

Figure 1 - 4 :

 14 Figure 1-4: Between-event and within-event components of ground-motion variability (L. Al Atik et al. 2010).

Figure 1 - 5 :

 15 Figure 1-5: The FAS of the 𝝎-squared model for displacement. 1.1.4.2 Finite-source rupture models

  a factor depending on the fault geometry). The stress drop (∆τ) is consequently expressed by: ∆τ = Gℇ ≈ G D ̅ L

Figure 1 - 6 :

 16 Figure 1-6: The barrier model and asperity model, respectively, for the aftershock and foreshock processes, after (Aki 1984).

Figure 1 - 7 :

 17 Figure 1-7: Relationship between normalized seismic moment M 0,n and fault roughness σ RMS . (a) Ruptures with near-center strength asperities increase M 0,n relative to the homogeneous strength case. Near-edge strength asperities have the opposite effect. (b) We observe an inverse relationship between σ RMS and M 0,n : rough faults release less seismic moment per tectonic loading than smooth ones, after (Zielke, Galis, and Mai 2017).

Figure

  Figure 1-8: Slip distribution for the 1940 Imperial fault rupture, after (Rockwell and Klinger 2013).

Figure

  Figure 1-10.

Figure 1

 1 Figure 1-9: A-section of partly eroded large slip fault surface at Mirrors locality on Dixie Valley fault, Nevada. B-Light detection and ranging (Lidar) fault surface topography as color-scale map, after (Sagy, Brodsky, and Axen 2007).

Figure

  Figure 1-10: Surface topography of the Castro Area fault, after (Candela, Renard, Schmittbuhl, et al. 2011).

Figure 1 -

 1 Figure 1-11: Selected finite-source rupture models for the 1999 Mw 7.6 Izmit earthquake (Turkey), obtained using different inversion strategies and different datasets. Black stars mark the hypocenter. Colors indicate fault slip (in meters), after (P.M. Mai et al. 2016).

Figure 1 -

 1 Figure 1-12: Slip distributions from finite-source studies illustrate the spatial variability of slip on the rupture plane, after (Mai and Beroza 2002).

Figure 1 -

 1 Figure 1-13: Corner frequency versus seismic moment (lower scale) and moment magnitude (upper scale). The dashed lines show constant stress drop of 0.1, 1, 10 and 100 MPa, after (Allmann and Shearer 2009).

Figure 1 -

 1 Figure 1-14: Joint distributions of stress drop and rupture speed values inferred from a database of apparent source time functions for ~100 shallow earthquakes, after (Chounet et al. 2017).

  .

Figure 1 -

 1 Figure 1-15: Distributions of kinematic source motions derived by dynamic rupture modeling, for two different input distributions of fracture energy. The

Figure 1 -

 1 Figure 1-16: Between-event residuals of PGA resulting from the nonparametric regression approach versus stress drop. Data points are color-coded following their region of origin in Japan. The correlation coefficients between stress drop and between event residual is 0.8, after (Oth, Miyake, and Bindi 2017).

Figure 1 -

 1 Figure 1-17: Stress parameters versus moment magnitude, after (Edwards and Fäh 2013).

Figure 1 -

 1 Figure 1-18: The acceleration ground motions resulting from different spatial distribution of the slip along the rupture area for different corner wavenumber of the slip, after (Causse et al. 2009). (Causse et al. 2009) used slip distributions following a k-2 spectral decay for different corner wavenumber values and showed examples of simulated accelerograms, assuming a constant rupture velocity. The case K=0.5 leads to large values of the correlation length of the slip distribution and hence smooth slip distributions, and results in low-amplitude ground acceleration (Figure 1-18).(Archuleta and Crempien 2015) simulated pseud-dynamic model and performed a sensitivity analysis of the ground motion to the slip correlation lengths using kinematic rupture simulations: they simulated synthetic white noise scaled to the seismic moment over the fault plane, and filtered it in the wavenumber domain with a Von Karman power spectrum based on the values of the slip correlation lengths. They found that as the asperity size decreases, the ground motion amplitude and its

Figure 1 -

 1 Figure 1-19: The acceleration ground motions resulting from different spatial distribution of the slip along the rupture area for different correlation lengths, after (Archuleta and Crempien 2015).

Figure 2 - 1 :

 21 Figure 2-1: Determination of the frequency controlling acceleration and velocity for the 1999 M=7.1 Ducze earthquake recorded at station 'Lamont 531'(Rrup =8 km, Vs30 = 640 m/s), expressed in terms of: (1) frequency associated with the maximum of Stockwell-transform amplitude at the time of PGA (resp. PGV) (fmaxST); (2) central frequency (Fcental) at the time of PGA (resp. PGV). The Central frequency is the 1 st order moment of the frequency distribution at the time of PGA (resp. PGV). The displayed component is the one associated with the

  1 and 5 s. The acceleration phase time 𝑇 𝑎𝑐𝑐 of the slip rate function is fixed to 10% of the rise time. The fault area is embedded at 0.5 km below the surface. The spatial distribution of final slip and rupture speed are tapered so as to avoid stress singularities at the fault boundaries. The values of the source parameters decrease with a quarter circular taper as they reach 20 % of the fault dimension at each side. Finally, we fix the hypocenter position to 20 % of the rupture length along the strike and 80 % along the dip, such that the rupture propagates unilaterally. This is consistent with[START_REF] Mai | Hypocenter Locations in Finite-Source Rupture Models[END_REF], who show that rupture initiate preferentially at depth, and with (McGuire, Zhao, and Jordan 2002), who obtained that 80 % of earthquake rupture are unilateral.

Figure 2 - 2 :

 22 Figure 2-2: Illustration of a rupture realization on the vertical fault plane, corresponding to a M7 event, and location of stations S1, S2, S3, S4 and S5. The white dotted lines represent the rupture front.

Figure 2

 2 Figure 2-3: a) to d). Spatial distribution of the kinematic rupture parameters (final slip D, rupture speed Vr, rise time Trise and peak slip velocity psv) for a homogeneous rupture of a M7 event, with L = 70 km, Vr = 2.71 km/s, Trise = 4.06 s, and psv = 106 cm/s tapered at the boundaries. The white contour lines represent the rupture propagation (each second) and the green contouring shows the location of the rupture front at the PGA time. e-1). Acceleration calculated in the far-field approximation using Equation (2-4), corresponding to the scenario described in a) to d), for X = 100 km and Qs = 220. This acceleration is referred to as the reference case. e-2). Accelerogram obtained using a wider taper at the fault boundaries (30 % tapering instead of 20 %) while the mean values of the kinematic parameters remain unchanged. e-3). Accelerogram due to a smaller rupture length while keeping M = 7 (55 km instead of 70 km); the stress drop, calculated using Equation (2-3) is increased by a factor of 1.29. e-4). Accelerogram computed for a faster rupture speed Vr =3.25 km/s. e-5).

Figure 2 -

 2 Figure 2-4a and 4b show two different realizations of heterogeneous ruptures with the same statistical properties of the source parameters (Table 2-2 -simulation B). These source models are associated with two different mechanisms of PGA generation. The PGA on Figure2-4a is induced by the rupture stopping at the top fault boundary, as observed for a homogeneous rupture. The PGA is however higher (0.041 g instead of 0.0147 g) because the rupture speed is heterogeneous and gives rise to a stronger rupture speed drop in this case. Furthermore, the PGA on Figure2-4b is not controlled by the same process because the rupture speed at the top fault above the nucleation is lower. It is generated by the large rupture speed patch located at the right edge of the rupture, resulting in an abrupt change of rupture velocity. Thus, the position of the high rupture speed patches with respect to the rupture nucleation and their interactions with the fault boundaries play a fundamental role in the PGA generation. Figure2-5 shows the distribution of the PGA time among the various rupture realizations. We consider only rupture models with Vr = 3 km/s. Interestingly, the probability density is maximum at about 3.5 s, which is the average time needed for the rupture to reach the top fault. This implies that the main mechanism and the main fault area implied in the PGA generation remain the same as for homogeneous ruptures.

Figure 2 - 4 :

 24 Figure 2-4: Two realizations of heterogeneous ruptures of a M7 event, with L = 70 km, 𝝁 𝑽𝒓 = 3 km/s, and 𝝁 𝑻𝒓𝒊𝒔𝒆 = 4 s, tapered at the boundaries and the corresponding acceleration computed using Equation (2-4). The statistical parameters are 𝝈 𝑫 = 𝟎. 𝟓 𝝁 𝑫 , 𝝈 𝑽𝒓 = 𝟎. 𝟓 𝝁 𝑽𝒓 , ax = 4 km, 𝝆 𝑫-𝑽𝒓 = 0.

  -6 and Figure2-7, we show the sensitivity of PGA and PGV to each source parameter for the five stations of Figure2-2 as well as in the case of the far-field approximation. The statistical percentiles (16 th , 50 th and 84 th ) of the PGA and the PGV sensitivity computed over the 150 rupture realizations are represented as colored patches.

Figure 2 - 6 :

 26 Figure 2-6: Sensitivity of the PGA to the different kinematic rupture parameters at stations S1 to S5, as well as for the far-field approximation. The 16 th , 50 th and 84 th percentiles of the PGA sensitivity computed for 150 simulations (50 realizations for each of the three specified hypocenters) are indicated for each couple of station and source parameter as colored patches. The median value of PGA sensitivity larger than 0.2 is shown explicitly.

Figure 2 - 7 :

 27 Figure 2-7: Same as Figure 2-6 for the PGV.

Figure 2

 2 Figure 2-8: (a) Coefficient of correlation between the values of PGA for all the combinations of the hypocenter locations at 20 % 'L', 50 % 'C', and 80 % 'R', for all the cases and for all the stations. (b) PGA values at stations S1 and S2 for case 2 used to compute the correlation coefficient in (a) marked by the blue rectangles, The variability of the PGA values is due to the spatial aleatory distribution of the source parameters along the fault, for the hypocenter located at 'L', 'C', and 'R' from the fault length.

  -9 and Figure 2-10.

Figure 2

 2 Figure 2-9: PGA variability for the different cases defined in Table 2-2 and for the different stations. The grey box delimits the values of the between-event variability obtained in some recent GMPEs (Causse and Song 2015). The parameter of interest in each case is: 𝝁 𝑽𝒓 (case 2, reference case 1), 𝝈 𝑫 (case 3, reference case 1), 𝝈 𝑽𝒓 (case 4, reference case 1), 𝒂 𝐱 (case 5, reference case 1), D-Vr + (case 6, reference case 4), D-Vr -(case 7, reference case 4) and ∆𝝉 (case 8, reference case 1).

Figure 2 -

 2 Figure 2-10: PGV variability for the different cases defined in Table 2-2 and for the different stations. The grey box delimits the values of the between-event variability obtained in some recent GMPEs (L. Al Atik et al. 2010). The parameter of interest in each case is: 𝝁 𝑽𝒓 (case 2, reference case 1), 𝝈 𝑫 (case 3, reference case 1), 𝝈 𝑽𝒓 (case 4, reference case 1), 𝒂 𝐱 (case 5, reference case 1), D-Vr + (case

Appendix 2- 3 :

 3 The convergence of the variability of PGA for 50 simulations. The Blue curves represent the evolution of the 𝝈 𝒍𝒏(𝑷𝑮𝑨) for a hypocentre located at 20% from the fault length, and the Red curves represent the evolution of 𝝈 𝒍𝒏(𝑷𝑮𝑨) for the hypocentre located at 80% form the fault length. Each column represents the station that makes symmetrical cases with the hypocenter positions. Each row represents the source scenarios defined in Table 2convergence of the variability of PGV for 50 simulations. The Blue curves represent the evolution of the 𝝈 𝒍𝒏(𝑷𝑮𝑽) for a hypocentre located at 20% from the fault length, and the Red curves represent the evolution of 𝝈 𝒍𝒏(𝑷𝑮𝑽) for the hypocentre located at 80% form the fault length. Each column represents the station that makes symmetrical cases with the hypocenter positions. Each row represents the source scenarios defined in Table2

Figure 3 - 1 :

 31 Figure 3-1: Distribution of the stations with respect to the top-fault rupture projection for the synthetic dataset of strike-slip rupture simulations. The black lines represent iso-values of the probability to observe a pulse, as given by (Shahi and Baker 2011).

Figure 3

 3 Figure 3-2: (a) Description of the rupture propagation. Slip is supposed to be constant on the fault plane (equal to 𝟏 𝒎) or heterogeneous. The displayed slip distribution has a correlation length of 𝟐𝟎 𝒌𝒎 along strike and 𝟕 𝒌𝒎 along dip, and has a standard deviation of 𝟎. 𝟕𝟓 𝒎. The white star represents the rupture initiation. The rupture speed is constant. (b-left) Simulated velocity at station 4 (fault-normal component) for homogeneous or heterogeneous rupture, or considering the rupture area represented by the shaded rectangle. (b-right) Extracted pulses and values of the pulse periods using Baker (2007) algorithm at station 4. The color version of this figure is available only in the electronic edition.

Figure 3 - 3 :

 33 Figure 3-3: Illustrative scheme describing the parameters used in Equation (2-1). 𝑯𝒚𝒑𝑫 denotes the hypocentral distance, 𝒄𝒍𝒔𝑫 denotes the closest distance to the rupture area, 𝑫 denotes the length of the rupture area that breaks toward the site, and 𝑽𝒓 is the rupture velocity. The color version of this figure is available only in the electronic edition.

  Figure 3-4: Values of the synthetic velocity pulse periods at the 12 stations for three different rupture scenarios with uniform slip. The circles represent the extracted pulse periods using the wavelet approach (Baker 2007), denoted by Tp (sim) and the filled squares represent the calculated pulse periods using Equation (2-1), denoted by Tp (eq1).

Figure 3 - 5 :

 35 Figure 3-5: Simulated velocity time series (black line) (fault-normal component) and extracted pulses (gray line) using Baker (2007) algorithm. The slip is assumed to be constant over the fault plane. Tpmeas. denotes the period of the extracted pulse, while Tpred. denotes the pulse period predicted from Equation (2-1).

Figure 3 - 6 :

 36 Figure 3-6: Database of pulse periods from 16 different earthquakes at 110 stations. The circles represent the extracted pulse periods using the wavelet approach of (Baker 2007) as listed in the NGA-West2 database, denoted by Tp (NGA), and the filled squares represent the calculated pulse periods using Equation (2-1), denoted by Tp (eq1).

Figure 3 -

 3 Figure 3-7: (a) Value of the pulse period Tp (NGA) as a function of moment magnitude. (b) D as a function of moment magnitude for the NGA-West2 dataset considered in this study (see Appendix 3-1).

Figure 3

 3 Figure 3-8: (a) Spatial distribution of the mean pulse period value and (b) standard deviation of the natural logarithm pulse periods for 𝑽𝒓/𝑽𝒔 and 𝑻𝒓𝒊𝒔𝒆 following a uniform distribution in the range [𝟎. 𝟔 -𝟎. 𝟗] and [𝟎. 𝟏 -𝟏] , respectively. The thin black lines in figure (a) represent the iso-values of the probability to observe a velocity pulse, as defined by Shahi & Baker (2011). The star represents the hypocenter and the thick black lines, the fault surface projection. The color version of this figure is available only in the electronic edition.

Figure 4 - 1 :

 41 Figure 4-1: Schematic map of the Arabia-Africa-Eurasia zone of plate interaction, after (Reilinger et al. 2006). Double lines are extensional plate boundaries, plain lines are strike-slip boundaries, and lines with triangles are thrust faults. Dark numbers are slip rates (mm/yr) on block-bounding faults (number in parentheses are dip slip and those without are strike slip). White arrows and corresponding numbers are plate velocities (mm/yr) relative to Eurasia. Curved arrows show sense of block rotation relative to Eurasia. AE: Aegean plate, AN: Anatolian plate, CAU: Caucasus plate.

Figure 4 -

 4 Figure 4-2: a) The Dead Sea Fault system running from the south from the Red Sea to the north to Turkey. b) The main fault branches of the Dead Sea Fault crossing Lebanon and the historical large earthquakes (M > 6.7), modified after (Daëron et al. 2007).

Figure 4 - 3 :

 43 Figure 4-3: Local Seismicity map for Lebanon (2006-2013), by Lebanese CNRS (http://www.cnrs.edu.lb/).

Figure 4 - 4 :

 44 Figure 4-4: Contour map of PGA with a 10% probability of exceedance in 50 years, after (Huijer, Harajli, and Sadek 2016).

Figure

  Figure 4-5: major geological units and faults (bold) of Lebanon by (Grant, Wartman, and Abou-Jaoude 2016) Mount-Lebanon Range (1) and its foothills are bound by the Yammouneh fault to the East, Mount-Lebanon Thrust (MLT) fault to the North, and the Roum fault to the South. The Bekaa Valley (2) divides the high Mount-Lebanon Range from the lower Anti-Lebanon Range (3). Rolling southern alluvial valleys and lowlands (4) extend across Lebanon south of the Roum fault.

Figure 4 - 6 :

 46 Figure 4-6: Main wadis and springs in the Lebanon mountain area and the Bekaa. After (Khair et al. 1994).

Figure 4 - 7 :

 47 Figure 4-7: Geological cross-sections along the 2D seismic profiles: A. Batroun -Ainata profile crossing the Qartaba Structure; and B. Aley -Barr Elias profile crossing the Levant Fracture and the Bekaa valley. The inset map shows the location of the sections and the seismic profiles by (Nader, Browning-Stamp, and Lecomte 2016).

Figure 5 - 1 :

 51 Figure 5-1: The Dead Sea Fault system running from the south from the Red Sea to the north to Turkey. b) The main fault branches of the Dead Sea Fault crossing Lebanon, after (Daëron et al. 2007).

Figure 5 - 2 :

 52 Figure 5-2: Major geological units and faults of Lebanon, and the distribution of the 21 seismic stations recording ambient noise in 2011-2012, modified after (Grant, Wartman, and Abou-Jaoude 2016). Mount-Lebanon Range (1) and its foothills are bound by the Yammouneh fault to the East, Mount-Lebanon Thrust fault to the North, and the Roum fault to the South. The Beqaa Valley (2) divides the high Mount-Lebanon Range from the lower Anti-Lebanon Range (3). Rolling southern alluvial valleys and lowlands (4) extend across Lebanon south of the Roum fault.

Figure 5 - 3 :

 53 Figure 5-3: Stacked cross-correlation of signals recorded at LU12 and LU08 seismological stations band-pass filtered between various period bands (top panel) and between 1 and 25 s (bottom panel). The red box delimitates the time range corresponding to velocities between 1 and 5 km/s, the blue box encompasses the time band of the envelope of the signal and the yellow box represents the time band of the noise.

Figure

  Figure 5-4: Rayleigh wave group velocities measured between selected pairs of stations at a) 1 s, b) 4 s, c) 7 s, d) 12 s, e) 18 s and f) 22 s of periods. Black lines indicate the borders of Lebanon.

Figure

  Figure 5-5: Local Rayleigh wave group velocity map at a) 1 s, b) 4 s, c) 7 s, d) 12 s, e) 18 s and f) 22 s od periods. Black lines indicate the borders of Lebanon.

Figure 5 - 6 :

 56 Figure 5-6: Average Shear-wave velocity VS profiles (red curve) with the corresponding standard deviation of VS with depth. The solid black curves represent the VS profiles of the surrounding regions from (Pinsky et al. 2013) and the dashed lines refer to the

  Figure 5-7 (a) and indicates three main seismic contrasts at around 2-3 km, 10-15 km and 20 km throughout Lebanon. The average VS profile from the individual VS estimates and the corresponding standard deviation are plotted in Figure 5-6. The inverted VS profiles fit with the range of VS values found in the surrounding regions. Horizontal sections of the VS tomography at different depths is presented in Figure 5-8 and discussed in details later. 5.4.4.2 The uncertainty on inverted VS profiles

  . Sections of the 3D COV at different depth are presented in Figure5-9. Large values of COV are observed around 2-3 km, related to the uncertainty on the inverted seismic depth interface. For depths larger than 20 km, large COV values are also observed as a consequence of the decrease surface waves resolution for largest wavelengths.

Figure 5 - 7 :

 57 Figure 5-7:The inverted 3-D VS model under Lebanon area and the corresponding resolution in terms of coefficient of variation of the shear wave velocity.

Figure 5 - 8 :

 58 Figure 5-8: Shear wave velocity maps at different depths in km/s. The black lines denote Lebanon borders and the fault within the region. The black dots denote the stations.

Figure 5 - 9 :

 59 Figure 5-9: The resolution in terms of coefficient of variation of the shear wave velocity maps at different depths. The black lines denote Lebanon borders and the fault within the region. The black dots denote the stations.

5. 5 . 1

 51 Spatial variation of VS at different depths First, the spatial variation of VS at various depths is shown in Figure5-8, while Figure5-9 indicates the corresponding COV. At z =1 km, a Vs ~2.5 km/s zone outspread at the Mount-Lebanon range (Figure5-1) surrounded by a Vs ~2.3 km/s zone. Low shear-velocity patches (VS ~2.1 km/s) trace in the north-west as well as in the south-west. The presence of such low VS matches the features of the surface geology of Lebanon. The landscape of Lebanon is principally dominated by the two mountain chains striking N-NE, and in between the Beqaa Valley (associated to the Vs ~2.5 km/s zone). On the northern coast of Lebanon, the Aakar plain, Tripoli coastal plain and foothills of the Lebanon mountains, and Chekka narrow coastal plain extend (Figure5-1) in front of the mountain range (associated to the Vs ~2.3 km/s zone).

  Figure 5-10 (b) and (c) reveal the surface topography for each of the 2 profiles; (d) and (e) illustrate the results of our VS tomography up to 8 km, with the corresponding resolution in Figure 5-11. Profile 1shows a top layer of VS ~2.5 km/s, and about 2 km thick, overlaying a 1.5 km average thickness unit with VS ~2.65 km/s, and then a layer of VS ~3 km/s. Profile 2 shows the same stratigraphy of the layers, however, the Yammouneh Fault trace is well defined, separating the west side from the east side under which a thick layer of ~3 km of VS ~2.8 km/s is clearly detected.

Figure 5 -

 5 Figure 5-10: a) the geological map of Lebanon showing the location of the profiles A and B in white, and profile 1 and profile 2 in magenta. b) and c) the tomography of profiles 1 and 2. d) and e) The vertical cross-section of the profiles 1 and 2 resulting from this work. f) and g) The geological cross-sections of profiles 'A' and 'B' after (Nader, Browning-Stamp, and Lecomte 2016). The interface between the Jurassic and Triassic denoted in black dots is reproduced on the tomography plots (d) and (e) as dashed lines. The lower dashed line represents the interface corrected by the topography.

Figure

  Figure 5-11: The corresponding resolution of the tomography shown in Figure 5-10.

  assessment of group velocity maps using multi-scale checkerboard tests. a(1-2): input models with velocity anomalies (dv/v) of size 0.2 • and 0.4 • ; b(1-2): outputs for period of 2 s; c(1-2): outputs for period 8 s; d(1-2): outputs for period 20 s. The black lines represent the borders of Lebanon. The black dots represent the fault lines.

Figure 6 - 1 :

 61 Figure 6-1: Map showing the population density distribution in Lebanon measured as the number of persons per square kilometer of land area, after the online Earth-data pool EOSDIS. A hypothetical rupture length of 70 km considered in this study on the Yammouneh Fault is indicated as black line, the stars representing the hypocenters in study (L, C and R). The colored dots represent the location of the stations where the ground motion is synthetized.

Figure 6 - 2 :

 62 Figure 6-2: Flowchart of the hybrid methodology applied for the broadband ground-motion simulation. F: type of faulting, Mw: earthquake magnitude, ZTor: depth to top of rupture plane, Rrup: source-to-site distance, Vs30: site characteristics, sord and ΘorΦ: directivity parameters, Vp: pulse amplitude, Tp: pulse period, ϒ: parameter characterizing the number of oscillations within the pulse, Φ: phase angle, and tmax,p: time of the peak of the modulating envelope. 6.5.1.2 1D shear-wave velocity medium in Lebanon

  at two rock sites: Bhannes and the campus of the American University of Beirut (AUB) downtown Beirut. Phase velocities were available from 4 to 10 Hz (AUB) and from 40 to 75 Hz at Bhannes. From suites of best inverted Vs profiles shown in Figure6-3, we then adopt a representative Vs profile for the surface to the first km depth for rock sites with Vs-30 = 800 m/s (rock-like site classified as soil type B in the UBC97 code). The 1D propagation medium properties finally used in the ground motion simulation are shown in Table6-2.

Figure 6 - 3 :

 63 Figure 6-3: Ensemble of Vs profiles inverted from small aperture seismic ambient noise arrays at two sites: Bhannes and the campus of the American University of Beirut (AUB).

Figure 6 - 4 :

 64 Figure 6-4: Pulse detection and extraction of its characteristic parameters. The blue solid lines show an example of a generated low frequency ground motion at station S3 for a subshear rupture, and blue dashed lines show an example of a generated low frequency ground motion at station S3 for a supershear rupture. The dashed curves represent the fitted pulses respectively. The pulses are obtained from a wavelet decomposition (see (Baker 2007) for details).

Figure 6 - 5 :

 65 Figure 6-5: On the left, source rupture simulation showing an example of the spatial distribution of the final slip, the rupture speed and the rise time along the fault. The white dotted plots on each illustrate the rupture propagation each 1 s. On the right, two examples of the 2 components of horizontal acceleration, showing two examples of the low-frequency simulated pulse-like ground motion recorded at station 3, for subshear rupture, in red, with their corresponding broadband stochastic ground motions in black. The response spectra issued from 50 stochastic simulations for this scenario are also plotted. The red dashed vertical line represents the period of the pulse detected.

Figure 6 - 6 :

 66 Figure 6-6: On the left, source rupture simulation showing an example of the spatial distribution of the final slip, the rupture speed and the rise time along the fault. The white dotted plots on each illustrate the rupture propagation. On the right, two examples of the 2 components of horizontal acceleration, showing two examples of the low-frequency simulated pulse-like ground motion recorded at station 3, for supershear rupture, in red, with their corresponding broadband stochastic ground motions in black. The response spectra issued from 50 stochastic simulations for this scenario are also plotted. The red dashed vertical line represents the period of the pulse detected.

Figure 6 - 7 :

 67 Figure 6-7: On the left, source rupture simulation showing an example of the spatial distribution of the final slip, the rupture speed and the rise time along the fault. The white dotted plots on each illustrate the rupture propagation. On the right, two examples of the 2 components of horizontal acceleration, showing two examples of the low-frequency simulated non pulse-like ground motion recorded at station 3, for subshear rupture, in grey, with their corresponding broadband stochastic ground motions in black. The response spectra issued from 50 stochastic simulations for this scenario are also plotted.

Figure 6 -

 6 Figure 6-8: (Top panel) in blue, the GMRoTD50 pga (g) values for the 50 ground motions generated at low frequency GM1Hz (up to 1 Hz, deterministic method coupling rupture modeling and a 1D velocity structure) for the 6 different rupture case scenarios L, C, R, Ls, Cs, Rs computed at station S1, and (bottom panel) the pulse periods associated to their acceleration time series. In red, the GMRoTD50 pga (g) values for the 50*50*6 ground motion up to 10 Hz (GM10Hz). The median values and the one standard deviation of the pga (considering a lognormal distribution) issued from all the broadband simulations are shown as solid and dashed horizontal lines.

Figure 6 -

 6 Figure 6-12: PGA one-point statistics (mean and standard deviation) computed at stations 1 through 4, for the different scenarios: L, C, R, Ls, Cs and Rs. On top, the PGAs are computed from the low-frequency deterministic simulations. At bottom, PGAs are computed from the hybrid broadband stochastic model up to 10 Hz.

Figure 6 -

 6 Figure 6-13: GMPE Response spectra for stations located at 5 km and 25 km away from the fault. The response spectra are computed for class B soil. HM is for hybrid method.

Figure 6 -

 6 Figure 6-14: Response spectra computed at each station, issued from the 15000 acceleration time histories (6x50x50), using the Duhamel integration. The red curve and the black curves represent the median and the standard deviation respectively. The design response spectra for Lebanon (for Beirut and all other regions) are also plotted in blue for comparison.

Table 1

 1

	-1.

Table 1 -1: Values of the Between-event variability of 𝑷𝑮𝑨 reported by some recent ground motion prediction equations for crustal events based on the Next Generation Attenuation-West2 database, except for Akkar and Bommer (2010), who used a European database, after (Causse and Song 2015).

 1

Table 1 -

 1 2: Compilation of correlation relationships between local-scale and large-scale source and fault parameters.

	Fault geometry	Stress and strength (friction parameters)	Kinematic rupture parameters	Data used to infer correlations	Correlati-on	References
			local-scale Parameters
			D -Vr	Finite source inversions	+	(Song et al., 2009)
				database of
			D-Vr	dynamic	0	(Schmedes et al., 2010)
				rupture simulations
				database of
			D-Vr	dynamic	0	(Mai et al., 2017)
				rupture simulations
			D -psv	Finite source inversions	+	(Song et al., 2009)
			D-risT	Finite source inversions	+	(Song et al., 2009)
				database of
			D-risT	dynamic	+	(Schmedes et al., 2010)
				rupture simulations
				database of
			D-risT	dynamic	+	(Mai et al., 2017)
				rupture simulations
			Vr -psv	theory	+	(Ida, 1973)
				database of
			Vr -psv	dynamic	+	(Schmedes et al., 2010)
				rupture simulations

Table 2 -1: Characteristics of the 15 accelerograms available in the PEER database (https://ngawest2.berkeley.edu/, last access June 27th, 2018) for 6.7 < M < 7.3, Vs30 > 600 m/s and Rrup < 50 km. The frequencies controlling the PGA and PGV values are expressed in terms of: (1) frequency corresponding to the maximum amplitude of the Stockwell-transform (Stockwell, Mansinha, and Lowe 1996) at the time of PGA (resp. PGV); (2) Central frequency (that is, first order moment of the frequency distribution at the time of PGA (resp

 2

. PGV)). Error! Reference source not found. displays

an example of Stockwell-transform of elected near-fault accelerograms. The symbols in the table are defined in the following: M: the moment magnitude, Rrup: the distance to the rupture, fc-a: the estimated corner frequency, fcentral: represents the central frequency of acceleration, fm-a: the frequency of maximum amplitude of acceleration S-transform, fc-v: the estimated corner frequency, fm-v: the frequency of maximum amplitude of velocity S-transform.

	Landers 06/28/1992 (SilentvValley -Poppet Flat)	7.3	51	660	0.2	9.4	8	5	2.7
	Landers 06/28/1992 (Twentynine Palms) 2.3 Mechanism of PGA generation in kinematic source models 7.3 41 635 0.2	6.7	4.5	5.5	4
	Manjil 06/20/1990 (Abbar) 2.3.1 Earthquake source model	7.3	12	725	0.4	10	11	4	0.5
	Strong motion data	M	Rrup (km)	Vs30 (m/s)	fc-a (Hz)	fcentral (Hz)	fm-a (Hz)	fc-v (Hz)	fm-v (Hz)
	Kobe 01/16/1995 (Chihaya)	6.9	50	610	0.1	6.1	3.7	4	3.5
	Kobe 01/16/1995 (MZH)	6.9	70	610	0.15	3.5	1.7	2.4	1.7
	Kobe 01/16/1995 (Nishi-Akashi)	6.9	7	610	0.2	3.8	2.2	2.1	2.1
	Duzce 11/12/1999 (Lamont 1060)	7.1	26	780	0.15	5.6	5.4	1.3	0.5
	Duzce 11/12/1999 (Lamont 531)	7.1	8	640	0.1	4.6	4	2.9	1.7
	Hector Mine 10/16/1999 (Heart Bar State Park)	7.1	61	625	0.1	5.8	6.2	3.6	3.2
	Hector Mine 10/16/1999 (Twentynine Palms)	7.1	42	635	0.2	8.9	5.5	3.8	1.8
	Hector Mine 10/16/1999 (Hector)	7.1	12	730	0.15	4	1.9	1.3	0.7
	Hector Mine 10/16/1999 (Joshua Tree N.M. -Keys View)	7.1	50	690	0.1	4.3	3.7	3.5	3.7
	Darfield 09/03/2010 (RPZ)	7.0	57	640	0.15	5.3	2.4	1.8	1.3
	Darfield 09/03/2010 (CSHS)	7.0	43	640	0.5	3	1.9	1.4	0.6
	Darfield 09/03/2010 (LPCC)	7.0	26	650	0.1	7	5.1	3.6	1.7

Table 2 -2: The source parameters

 2

: M, L, W, ∆𝝉, D,

Vr, ax, az, psv, and Trise stand for the magnitude, length, width, stress drop, slip, rupture speed, spatial correlation lengths along the strike and along the dip directions, peak slip velocity and rise time, respectively. 𝝁 represents the mean value, 𝝈 is the standard deviation and 𝝆 is the coefficient of correlation.

	2.3.5 Mechanism of PGA generation for heterogeneous ruptures
	Simulation id	M	L (km)	W (km)	μD (cm)	Δτ (MPa)	σD/μD	μVr (km/s)	σVr/μVr	ρD-Vr	ax (km)	az (km)	μpsv (cm/s)	σpsv	ρD-psv	risT (s)	risT min	risT max	circular taper
	A B	7 70	14	120 2.87	0 0.5	3	0 0.5	-0	-4	-2	106	0 50	0.8	f(D, psv)	0.1 5	L/5 & W/5
	1						0.5	3	0.25	0	16	5
	2						0.5	2.5	0.25	0	16	5
	3						1	3	0.25	0	16	5
	4		70		120 2.87	0.5	3	0.5	0	16	5
	5	7		14			0.5	3	0.25	0	4	2	115	50 0.8	f(D, psv)	0.1 5	L/5 & W/5
	6						0.5	3	0.5	1	16	5
	7						0.5	3	0.5	-1	16	5
	8		55		153	3.63	0.5	3	0.25	0	16	5

Appendix 2-1: 1-D velocity model used in our simulations.

	H (km)	Vp (km/s)	Vs (km/s)	Rho (g/cm 3 )	Qs
	0.00	6.20	3.58	2.70	220
	14.00	6.80	3.93	2.86	220
	34.00	8.05	4.65	3.28	220
	50.00	8.25	4.76	3.29	220
	80.00	8.50	4.91	3.50

Table 3 -1: List of earthquakes considered in the present study. Earthquake Year Mw 1 L 2 Trise 3 Vr 4 Vr/Vs 5 Stations

 3

	Coyote Lake	1979 5,74	6,6	0,43 2,68 0,77	4	SS 6
	Parkfield	2004 6,00	40	0,88 2,7	0,8	11	SS
	Morgan Hill	1984 6,19	27	0,43 2,58 0,80	2	SS
	Imperial Valley	1979 6,53	50	0,87 2,70 0,87	12	SS
	Superst. Hills	1987 6,54	20	0,62 2,47 0,77	2	SS
	Bam, Iran	2003 6,60	15	1,40 2,80 0,92	1	SS
	Kobe, Japan	1995 6,90	60	1,24 2,68 0,80	5	SS
	Duzce, Turkey	1999 7,14 46,8 1,90 2,80 0,80	2	SS
	1 Earthquake Magnitude
	2 Fault length in km
	3 Rise time in s
	4 Rupture speed in km/s
	5 Ratio between the rupture speed and the shear wave speed
	6 Strike Slip

Table 3 -2: Rupture speed values for the 1995 Kobe earthquake according to various published source models. Source model reference for the Kobe earthquake Vr 8

 3

	(Zeng and Anderson 2000)	2,8
	(Yoshida et al. 1996)	2,5
	(Wald 1996)	2,8
	(Koketsu, et al., 1998)	2,5
	(Ide et al., 1996)	3
	(Horikawa et al., 1996)	3
	(Cho and Nakanishi 2000)	3,4
	(Yujia, Koketsu, and Ohno 2013)	3,1
	(Sekiguchi et al., 2000)	3,1

Appendix 3-1: Dataset of pulse periods considered in this study, build from the NGA-West2 database.

  Overview of the Dead Sea Fault in the Levant regionThe Dead Sea Transform (DSTF), also called the Levant Fault, is one of the active continental transform faults, like San Andreas Fault of California, Alpine Fault of New Zealand and North Anatolian Fault of Turkey. The DSTF began about 25 million years ago when Arabia was still part of the African continent. As rifting began to open the Red Sea, the Arabian plate split from Africa and began to move northward.

		Chapter 3: Spatial variability of the directivity pulse periods observed during an earthquake Chapter 3: Spatial variability of the directivity pulse periods observed during an earthquake Chapter 3: Spatial variability of the directivity pulse periods observed during an earthquake
		Imperial Valley Northridge Chi-Chi, Taiwan	Agrarias Pacoima Kagel Canyon 26,04 10,30 WGK 32,95	0,65 7,26 9,94	10,28 2,338 19,88 0,728 32,08 5,74
		Imperial Valley Northridge Chi-Chi, Taiwan	El Centro -Me Geot Pacoima Dam (up. left) 26,85 21,84 TCU031 80,49	0,07 7,01 30,17	21,84 3,423 19,88 0,84 58,92 5,929
		Imperial Valley Northridge Chi-Chi, Taiwan	El Centro Array #6 LA -Sepulveda VA H TCU040	29,22 19,45 69,51	1,35 8,44 22,06	29,19 3,773 17,54 0,931 50,77 6,433
		Imperial Valley Northridge Chi-Chi, Taiwan	El Centro Array #5 Pardee -SCE TCU063	29,53 31,05 36,35	3,95 7,46 9,78	29,26 4,13 24,08 1,232 28,21 6,552
		Earthquake Imperial Valley Northridge Chi-Chi, Taiwan	Station Name El Centro Array #7 Rinaldi Receiving Sta CHY024	HypD 9 ClstD 10 D 11 29,38 0,56 29,34 4,375 Tp(NGA) 12 20,62 6,5 19,55 1,246 25,39 9,62 17,92 6,65
		Coyote Lake Imperial Valley Northridge Chi-Chi, Taiwan	Gilroy Array #3 Brawley Airport Newhall -Fire Sta TCU104	12,49 44,29 26,78 49,93	7,42 10,42 5,92 12,87	6,46 39,05 4,396 1,155 21,50 1,372 45,87 7,189
		Coyote Lake Imperial Valley Northridge Chi-Chi, Taiwan	Gilroy Array #6 EC County Center FF LA Dam TCU059	9,12 30,73 21,10 53,97	3,11 7,31 5,92 17,11	6,46 29,41 4,417 1,232 19,50 1,617 45,87 7,784
	Coyote Lake Coyote Lake Imperial Valley Imperial Valley Northridge Northridge Chi-Chi, Taiwan Chi-Chi, Taiwan 4 STATE OF ART Gilroy Array #4 Gilroy Array #2 El Centro Array #3 El Centro Array #10 Sylmar -OV Med FF Newhall -W P C Rd. NST TCU046	11,08 13,55 30,33 30,46 24,24 27,76 89,20 69,35	5,7 9,02 12,85 8,6 5,3 5,48 38,42 16,74	6,46 6,46 27,47 4,501 1,351 1,463 28,70 4,515 19,55 2,436 24,08 2,982 57,86 7,875 58,92 8,043
		Parkfield Imperial Valley Northridge Chi-Chi, Taiwan	Cholame 3E El Centro Array #4 Sylmar -Converter Sta TCU082	14,37 28,90 21,87 37,07	5,55 7,05 5,35 5,16	11,46 0,518 28,03 4,788 19,54 2,982 36,31 8,099
		Parkfield Imperial Valley Northridge Chi-Chi, Taiwan	Stone Corral 1E Holtville Post Office Jensen Filter Plant Adm 21,78 10,82 22,16 TCU026 106,51	3,79 7,5 5,43 56,12	8,61 20,85 4,823 0,574 19,58 3,157 57,50 8,372
		Parkfield Imperial Valley Northridge Chi-Chi, Taiwan	Cholame 4W El Centro Diff. Array Sylmar -Conv Sta East TCU064	14,74 29,00 22,16 59,68	4,23 5,09 5,19 16,59	11,46 0,7 28,23 6,265 19,51 3,528 49,82 8,456
	4.1	Parkfield Superstition Hills Kornbloom Road Slack Canyon Northridge Jensen Filter Plant Gen 21,80 32,55 21,27 Chi-Chi, Taiwan TCU103 53,04	2,99 18,48 5,43 6,08	30,52 0,854 10,54 2,128 19,60 3,535 50,77 8,687
		Parkfield Superstition Hills Parachute Test Site Cholame 3W Loma Prieta Los Gatos -Lex. Dam Chi-Chi, Taiwan TCU034	14,62 18,35 26,83 88,24	3,63 0,95 5,02 35,68	11,46 1,022 18,33 2,394 23,60 1,568 58,92 8,869
		Parkfield Bam, Iran Loma Prieta Chi-Chi, Taiwan	Cholame 2WA Bam Gilroy -Historic Bldg. TCU136	14,10 13,94 33,10 49,40	3,01 1,7 10,97 8,27	11,46 1,078 13,39 2,023 24,71 1,638 45,87 8,8816
		Parkfield Kobe, Japan Loma Prieta Chi-Chi, Taiwan	Fault Zone 9 KJMA Gilroy Array #2 TCU056	12,86 25,58 34,52 40,53	2,85 0,96 11,07 10,48	11,38 1,134 25,36 1,092 24,71 1,729 37,48 8,939
		Parkfield Kobe, Japan Loma Prieta Chi-Chi, Taiwan	Fault Zone 1 Takatori Gilroy Array #3 TCU033	11,67 22,19 35,94 93,41	2,51 1,47 12,82 40,88	10,09 1,19 21,88 1,554 24,71 2,639 58,92 8,974
		Parkfield Kobe, Japan Loma Prieta Chi-Chi, Taiwan	Fault Zone 12 Takarazuka Saratoga -Aloha Ave TCU128	13,66 42,55 32,35 63,80	2,65 0,27 8,5 13,13	12,30 1,19 42,47 1,806 24,74 4,571 58,92 9,023
		Parkfield Kobe, Japan Loma Prieta Chi-Chi, Taiwan	EADES Port Island (83 m) Saratoga -W V Coll. NSY	12,83 26,28 32,20 63,80	2,85 3,31 9,31 13,13	11,32 1,218 26 2,534 24,74 5,649 58,92 9,163
		Parkfield Kobe, Japan Cape Mendocino Centerville Beach, N F Cholame 1E Port Island (0 m) Chi-Chi, Taiwan TCU039	14,02 26,28 29,58 71,95	3 3,31 18,31 19,89	11,46 1,33 25,74 2,828 13 1,967 58,92 9,331
		Morgan Hill Duzce, Turkey Cape Mendocino Petrolia Coyote Lake Dam -SW 25,98 Bolu 43,58 10,52 Chi-Chi, Taiwan TCU045 77,91	0,53 12,04 8,18 26	25,82 1,071 32,61 0,882 6,67 2,996 58,92 9,338
		Morgan Hill Duzce, Turkey Cape Mendocino Bunker Hill FAA Gilroy Array #6 IRIGM 487 Chi-Chi, Taiwan TCU038	37,32 26,72 21,86 73,55	9,87 2,65 12,24 25,42	27,69 1,232 24,67 10,052 12 5,362 50,77 9,576
		Landers Tabas, Iran Chi-Chi, Taiwan	Lucerne Tabas TCU102	44,58 55,54 46,26	2,19 2,05 1,49	48,06 5,124 55,08 6,188 46,78 9,632
		Landers Chi-Chi, Taiwan Chi-Chi, Taiwan	Yermo Fire Station CHY006 TCU049	86,28 41,26 39,73	23,62 9,76 3,76	66,18 7,504 33,54 2,5704 39,77 10,22
		Landers Chi-Chi, Taiwan Chi-Chi, Taiwan	Barstow TCU TCU101	95,02 37,07 45,75	34,86 5,16 2,11	66,18 9,128 36,31 4,508 45,87 10,318
		Kocaeli, Turkey Chi-Chi, Taiwan Chi-Chi, Taiwan	Yarimca TCU076 TCU051	25,07 17,91 39,35	4,83 2,74 7,64	24,35 4,949 15,73 4,732 37,68 10,381
		Kocaeli, Turkey Chi-Chi, Taiwan Chi-Chi, Taiwan	Izmit TCU075 TCU087	16,86 22,16 56,21	7,21 0,89 6,98	15,41 5,369 21,68 4,998 58,92 10,395
		Kocaeli, Turkey Chi-Chi, Taiwan Chi-Chi, Taiwan	Gebze TCU029 TCU052	49,68 79,60 40,38	10,92 28,04 0,66	48,40 5,992 58,92 5,285 41,70 11,956
	Hypocenter Distance in km Kocaeli, Turkey Chi-Chi, Taiwan Chi-Chi, Taiwan	Arcelik CHY101 TCU068	56,02 32,95 48,52	13,49 9,94 0,32	49,16 7,791 32,08 5,341 50,96 12,285
	Closest Distance in km Distance between the hypocenter and the closest distance from the station to the fault rupture area in San Fernando Pacoima Dam (up. left) 17,60 1,81 13,44 1,638 Chi-Chi, Taiwan TCU036 68,28 19,83 50,77 5,383 Chi-Chi, Taiwan TCU047 86,76 35 58,92 12,313
	km Pulse period in s Northridge Chi-Chi, Taiwan Chi-Chi, Taiwan	Pacoima Dam (downstr) 26,85 TCU065 27,85 TCU053 41,97	7,01 0,57 5,95	19,88 0,588 27,93 5,74 41,34 13,118
			92 93 94

  [START_REF] Rezaeian | A Stochastic Ground Motion Model with Separable Temporal and Spectral Nonstationarities[END_REF] proposed a stochastic model that adequately represents the nonstationary characteristics of real earthquake ground motions both in time and frequency domains. Temporal and spectral nonstationarities are achieved through modulation in time (by multiplying the stochastic process with a deterministic function that varies over time) and by varying the filter parameters over time. The physically based parameters (Ia, D5-95, tmid, fmid, ξf, f') completely define the time modulation and the evolutionary frequency content of the nonstationary ground motion model: Ia represents the expected Arias intensity of the acceleration process, D5-95 represents the effective duration of the motion, tmid time at which 45% level of the expected Arias intensity is reached, fmid represents the filter frequency at tmid, f' represents the rate of change of the filter frequency with time, and ξf represents the filter damping ratio.[START_REF] Rezaeian | A Stochastic Ground Motion Model with Separable Temporal and Spectral Nonstationarities[END_REF] identified the model parameter values by studying recorded ground motion in the NGA_West2 database. Based on this database, empirical predictive equations for the model parameters are

Table 6 -2: Geophysical properties of the propagation medium

 6

	H (km)	Vp (km/s)	Vs (km/s)	ρ (g/cm 3 )
	0.00	1.83	0.8	1.8
	0.04	3.66	1.6	2
	0.08	5.45	2.38	2.20
	1.00	5.55	2.60	2.20
	2.00	5.62	2.90	2.40
	3.00	5.81	3.00	2.60
	8.00	6.10	3.15	2.80
	15.00	6.49	3.35	3.00
	25.00	8.72	4.5	3.30
	6.5.2 Broad-band ground motion (~0.1-10 Hz)

: Depth (H), Compressional wave speed (Vp), Shear wave speed (Vs), and the density (ρ).

Table 6 -3: The hybrid model (HM) resultant median and standard deviation values of the PGA simulated at the stations S1 through S4, compared to the GMPEs output at 5 km and 25 km from the rupture.

 6

		5 km				25 km
	PGAHM	S1		S2		S3		S4
	μ(g)	0.43		0.25		0.26		0.13
	μ± σ	0.28	0.67	0.17	0.35	0.16	0.35	0.08	0.20
	PGAGMPE	ASK14	BSSA14	CB14	ASK14	BSSA14	CB14
	μ(g)	0.36	0.34		0.67	0.37	0.4		0.47
	μ± σ	0.20 0.67	0.18	0.63 0.38	1.20 0.21 0.64	0.21	0.78 0.30	0.74

Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, 38000 Grenoble, France

Earthquake Research Center, Korea Institute of Geoscience and Mineral Resources (KIGAM), 124 Gwahang-no, Yuseong-gu, Daejon

34132, Korea.

Remerciements

Quelques mots de mon coeur, après trois années riches et intenses.

Acknowledgments

We thank an anonymous reviewer and Paul Spudish for valuable comments and insight that improved the article. This study has been partially supported by Université Grenoble Alpes (AGIR-POLE-PAGE project), by Institut National des Sciences de l'Univers (CNRS-Tellus-ALEAS project), by the STREST project (Seventh Framework Programme (FP7/2007-2013) under Grant Agreement Number 603389), and by the JEAI SAMMOVA (2015-2017) funded by Institut de Recherche pour le Développement (IRD).

, à ma mère, à mon père, à Anthony et à Joyce, pour leurs prières et leurs encouragements, pour leur soutien, ------------------------------------------------------------

Data and resources

The Next Generation Attenuation-West2 Project (NGA-West2) database was searched using http://peer.berkeley.edu/ngawest2/databases/ (last accessed March 2016). The rupture speed values reported by several finite-source inversion models were searched using http://equake-rc.info/SRCMOD/ (last accessed March 2016). The pulse periods of the synthetic velocity time series were computed using Shahi & Baker (2007) algorithm (MATLAB computer code available at https://github.com/shreyshahi/PulseClassification, last accessed March 2016).

 0 1 2 3 4 5 6 7 8 9 10 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 2 0 2 1 2 2 2 3 2