Voici un petit exemple de comment tracer des rais sismiques (pour impressionner vos collègues pendant des dîners mondains, évidemment). Je vous présente ici un petit tutoriel pour utiliser la librairie TauP disponible comme package Python. Cela fait partie du package ObsPy, une boite à outils (toolbox) pour les sismologues et les observations sismiques.

Si vous avez l’habitude des notebooks Jupyter, vous pouvez aller voir sur mon Github et télécharger le notebook de là bas. Et sinon, suivez le guide ici, en ouvrant un terminal et votre éditeur de texte préféré (comme le mien, c’est VIM, je vous laisse choisir)e

Python (3) et packages

Pour pouvoir me suivre, il faut Python 3, et installer quelques packages: numpy, matplotlib (les deux packages habituels pour tous les projets scientifiques), et obspy.taup. Pour l’installation, j’utilise Brew (sur Mac) et pip… Mais vous pouvez aussi faire une installation toute propre avec Conda.

# import statements
import numpy as np 
import matplotlib.pyplot as plt #for figures
from obspy.taup import TauPyModel #TauP

Initialisation

Commençons par initialiser le modèle… On va choisir un modèle de Terre qui s’appelle IASP91. Si vous allez voir la base de données de TauP, vous pouvez en choisir d’autres! C’est ce modèle qui vous donne les densités et vitesses sismiques en fonction du rayon.

model = TauPyModel(model="iasp91")

Ensuite, on va pouvoir déterminer trois choses:

  • la profondeur du séisme (tout sera fait « par rapport » à son emplacement, donc son emplacement est par défaut à la distance 0)
  • la position de la station sismique (par définition, à la surface, et on choisit la distance en degrés au séisme)
  • les phases qui nous intéressent (P, S, PKPKP, etc.)

Et on peut demander soit les temps d’arrivées, soit les trajectoires totales:

parameters = {"source_depth_in_km":55, "distance_in_degree":67, "phase_list":["P", "PSPSPS"]}

model.get_travel_times(**parameters)
model.get_get_ray_paths(**parameters)

Remarque sur Python: parameters est un dictionnaire, et l’utilisation des ** permet de « déplier » le dictionnaire et d’utiliser les valeurs du dictionnaire comme variables de la fonction.

Et pour visualiser les trajectoires, on demande simplement un plot():

arrivals = model.get_get_ray_paths(**parameters)
arrivals.plot()

Screen Shot 2018-02-19 at 19.59.48

Jouons un peu…

Pour les noms des phases, il faut se rappeler le nom de chaque phase dans chaque couche de la Terre:

  • P et S : ondes P et S dans le manteau
  • K : ondes P dans le noyau liquide (pas d’onde S)
  • I et J : ondes P et S dans la graine

Et les noms des réflections:

  • sans rien : réflections à la surface (typiquement: phase PP, PPP, etc.)
  • minuscule p et s : reflection à la surface, si l’onde partait « vers le haut »
  • c : à l’interface noyau/manteau
  • i : à l’interface graine/noyau liquide

En l’absence de phases données, le défaut est une série de phases pré-entrées. Voici par exemples toutes les phases que l’on obtiendrait avec une distance épicentrale de 70 degrés, et un tremblement de terre à 50km de profondeur:

arrivals = model.get_ray_paths(50, 70)
arrivals.plot()

Screen Shot 2018-02-19 at 20.10.09.png

Un exemple : la graine (PKIKP, PKiKP, PKP)

Et finalement, voilà un petit code pour faire une jolie figure sur l’observation de la graine terrestre:

fig = plt.figure(figsize=(20, 10)) 
ax1 = plt.subplot(121) 
ax2 = plt.subplot(122, projection='polar') 
ax2.set_theta_zero_location('N') 
ax2.set_theta_direction(-1)


distances = np.linspace(10, 180, 100)


phases = ["P", "PKP", "PKiKP", "PKIKP"]
colors = {"P":"b", "PKP":"r", "S":"g", "PKiKP":"k", "PKIKP":"k"}

for i, value in enumerate(distances):
   arrivals = model.get_travel_times(source_depth_in_km=source_depth,
                                     distance_in_degree=value,
                                     phase_list=phases)
   if len(arrivals)==0: continue
   for arrival in arrivals:
      ax1.scatter(value,arrival.time, c=colors[arrival.phase.name], s=10, linewidth=0)
 
ax1.set_ylabel("time (s)")
ax1.set_xlabel("distance (degree)")

for i, value in enumerate(distances[::10]):
   arrivals = model.get_ray_paths(source_depth_in_km=source_depth,
                                  distance_in_degree=value,
                                  phase_list=phases)
   if len(arrivals)==0: continue
      arrivals.plot(ax=ax2,plot_type='spherical', legend=False, label_arrivals=True,
                    show=False
inner_core.png
Trajectoires des ondes sismiques qui échantillonnent la graine, telles que possiblement observables entre 10 et 180 degrés, tous les 10 degrés.

Ici, vous remarquerez que j’ai fait fait deux visualisations différentes, dans deux subplots (les premières lignes du code). J’ai d’un coté les temps d’arrivées (le nombre de secondes entre l’évènement et l’arrivée de l’onde à la station) pour une centaine de stations, et de l’autre les trajets des ondes dans la Terre pour une dizaine de ces stations. On voit que en fonction de la distance épicentrale (en degrés), les ondes qui arrivent sont différentes. On peut donc y voir la zone d’ombre du noyau (la courbe bleue des ondes P ne rejoint pas celle rouge des PKP), et les réflexions à la surface de la graine (la courbe noire). Si vous regardiez la zone autour de 150 degrés d’un peu plus près, vous verriez qu’il y a superposition de plusieurs PKP et de PKiKP/PKIKP.

Conclusions:

Vous avez donc vu ici comment vous servir du package TauP pour impressionner vos amis, et passer pour un vrai sismologue! N’hésitez pas à aller voir les autres outils de ObsPy, ou même d’autres informations disponibles dans les fonctions de TauP.

Sur la structure de la Terre et les observations, on en parlait ici. Et pour en savoir plus,  sur le package obspy.TauP !

Et n’oubliez pas que ce sont des trajectoires pour des Terres qui ont une structure parfaite (pas de variations latérales!), et que surtout, ça ne vous dit pas l’énergie de l’onde qui arrive (pour les PKIKP, par exemple, pas grand chose), ou si votre onde sismique n’est pas « cachée » par une autre arrivant en même temps.

Publicités