Écouter les séismes : tracer les rayons d'une planète

Une onde S qui atteint la frontière noyau-manteau est finie. Le noyau externe est liquide, un liquide n'a aucune résistance au cisaillement, et une onde de cisaillement a besoin de quelque chose à cisailler — alors 2 891,5 km plus bas, la phase s'arrête. Ce que cela laisse de l'autre côté d'un grand séisme, c'est une zone d'ombre sismique de la taille d'un continent, et c'est l'une des observations qui nous ont appris que le milieu de la planète est fondu, des décennies avant que quiconque puisse le modéliser.

Cette ombre est à l'écran dans Tremor chaque fois que vous sélectionnez un événement important, et la façon dont elle y arrive est tout le principe en miniature. Aucune ligne ne dessine d'ombre. Aucun seuil sur la distance angulaire n'existe où que ce soit dans le moteur de rendu. Nous traçons des rayons dans un modèle de vitesse, nous demandons à chacun où il se trouve à l'instant de rejeu courant, et nous écrivons les réponses dans un tampon préalloué.

Les rayons morts au noyau écrivent NaN. Le rendu de lignes interrompt sa course sur un NaN et en commence une autre. Les trous sont les zones d'ombre.

Le modèle de temps de trajet ak135, en 22 nœuds

Tout l'aval repose sur une table. ak135 est un modèle de référence unidimensionnel standard de la Terre : pour chaque profondeur, une vitesse de compression et une vitesse de cisaillement moyennées sur la planète entière — assez approximatif pour tenir dans un paragraphe, assez bon pour que les sismologues localisent encore de vrais séismes avec. Notre version fait 22 lignes de [depthKm, vp, vs], interpolées linéairement par morceaux en rayon plutôt qu'en profondeur, à environ 0,1 km/s du modèle publié. Le code source parle de qualité visuelle et le pense : c'est un outil de dessin, pas une table de temps de trajet ak135 contre laquelle on pointerait une phase.

L'astuce de la table tient à la façon dont les discontinuités y sont encodées : une profondeur qui apparaît deux fois est une discontinuité. Le constructeur parcourt la liste de nœuds, la découpe en segments continus partout où une profondeur se répète, et enregistre le saut. Rien d'autre dans le moteur ne connaît les frontières par leur nom — le Moho, la 410, la 660, la frontière noyau-manteau et la frontière du noyau interne sont une seule structure de données, et en ajouter une autre reviendrait à ajouter deux lignes.

Le noyau externe liquide n'est pas un cas particulier dans le moteur. C'est un zéro dans une table.

Ce zéro abat un travail considérable. Quand un rayon S atteint la frontière noyau-manteau, l'étape de réfraction va chercher la vitesse de cisaillement de l'autre côté, ne trouve rien dans quoi se réfracter, et met fin au rayon — le traceur n'apprend jamais qu'il a heurté un noyau. Sous la frontière du noyau interne, la vitesse de cisaillement revient, parce que le noyau interne est solide ; le traceur s'y propagerait volontiers si une onde de cisaillement pouvait y parvenir. Aucune n'y parvient, et une assertion vérifie qu'aucune ne l'a jamais fait.

Courber 360 rayons à travers une planète

Une Terre à symétrie sphérique vous offre une dimension : un rayon partant d'une source ne quitte jamais le plan qui contient la source et le centre de la planète. Le traceur est donc bidimensionnel : il fait marcher chaque rayon dans ce plan de grand cercle et le stocke sous forme d'un rayon, d'une distance angulaire et d'un temps accumulé — trois flottants par échantillon conservé, dans un Float32Array plat, un échantillon gardé tous les six pas.

L'intégration relève de la théorie classique des rayons : avec un indice de réfraction n = 1/v, l'équation des rayons d/ds(n·d̂) = ∇n dit qu'un rayon se courbe vers l'indice le plus élevé — vers le milieu le plus lent. Dans un milieu à symétrie radiale, le gradient pointe le long du rayon vecteur, si bien que toute la courbure se réduit à une douzaine de lignes d'arithmétique scalaire, exécutées une fois par pas de huit kilomètres.

if (r > 0.5) {
  const dndr = -velocityGradientAt(r, phase) / (v * v);
  const k = (dndr * DS_KM) / r;
  const n = 1 / v;
  const wx = n * dx + k * x;
  const wz = n * dz + k * z;
  const wl = Math.hypot(wx, wz);
  if (wl > 0) {
    dx = wx / wl;
    dz = wz / wl;
  }
}

Un pas rectiligne convient jusqu'à ce qu'un rayon rencontre une discontinuité, où l'approximation se voit immédiatement : franchissez la frontière noyau-manteau quelques kilomètres trop tard et l'angle de réfraction est faux, et le rayon atterrit là où il n'a rien à faire. Les franchissements sont donc résolus exactement. stepToRadius() résout la quadratique qui place le rayon précisément sur la frontière, le pas s'y tronque, et la loi de Snell-Descartes est appliquée dans la base tangentielle/normale locale : sin θ₂ = sin θ₁ · v₂/v₁. Si cela dépasse un, il n'y a pas de rayon transmis et la direction se réfléchit — réflexion totale, et c'est ainsi que les phases réfléchies sur le noyau apparaissent sans qu'on les demande. Sinon, la direction réfractée est reconstruite et décalée de 0,02 km, pour que le pas suivant ne puisse pas redétecter le même franchissement.

Trois cent soixante rayons s'ouvrent en éventail depuis l'hypocentre, entre 0,25° et 179,75° d'angle de départ, la profondeur de la source étant bornée à 1–700 km. Un éventail complet coûte 10–20 ms par phase — sur le thread principal, sans Web Worker, ce qui n'est défendable que parce que cela survient une fois par sélection et non une fois par image.

Les phases que nous n'avons jamais écrites

La sismologie a un riche alphabet de phases — P, S, PKP, PKIKP, PcP, ScS — et la manière tentante de le rendre consiste à traiter chacune comme un cas particulier : dessiner cet arc-ci pour une phase de noyau, masquer la P directe au-delà de tel angle. Tremor n'a rien de tout cela. Le traceur connaît les vitesses, les gradients, les frontières et la loi de Snell-Descartes ; toutes les phases de l'image découlent de ces quatre choses.

L'éventail S s'arrête au noyau parce que la vitesse de cisaillement y est nulle. L'éventail P fait quelque chose de plus intéressant : la vitesse de compression chute de 13,66 à 8,0 km/s en entrant dans le noyau externe, si bien que les rayons se réfractent brutalement vers le bas, que les arrivées directes en surface s'épuisent au-delà de 104°, et que l'énergie ayant traversé le noyau ressurgit sous forme de PKP bien plus loin. Les deux bords de l'ombre P sortent de la réfraction, pas d'un if — et c'est précisément pour cela qu'ils demandent vérification.

Le moteur est donc livré avec un autocontrôle exécutable, lancé en développement, dans le même esprit que l'assertion ISS d'Orrery : un moteur physique que personne ne peut juger à l'œil a besoin d'un test qui, lui, le peut.

• Le temps de trajet croît de façon monotone le long de chaque rayon. (Un rayon qui remonte le temps signifie que l'intégrateur a cassé.)
• La P directe à 60° arrive vers 600 s, à ±60 s près ; à 90°, vers 780 s à ±80 s près.
• La S directe à 60° arrive vers 1 090 s, à ±110 s près.
• La dernière arrivée de P directe en surface tombe entre 90° et 118° — assez large pour tolérer un modèle de qualité visuelle, assez étroit pour attraper une réfraction cassée.
• Des arrivées PKP existent au-delà de 130°.
• Aucun rayon S n'est jamais trouvé à plus de 5 km sous la frontière noyau-manteau.

Un isochrone par image, zéro allocation

Un rayon tracé est un chemin dans l'espace et dans le temps ; un front d'onde en est une tranche. sampleIsochron() cherche par dichotomie, dans les échantillons de chaque rayon, l'encadrement contenant l'instant de rejeu courant, et en interpole une position — un couple rayon/distance par rayon, dans un tableau préalloué. Les rayons qui n'ont jamais atteint cet instant, ou qui sont morts au noyau, écrivent (NaN, NaN).

Les fronts sont dessinés en LineSegments plutôt qu'en polyligne, et ce choix porte tout. Une polyligne relierait le dernier échantillon vivant avant un trou au premier qui le suit, soudant une corde en travers de l'ombre. Les segments ne le font pas : un NaN termine une course, le couple valide suivant en ouvre une autre, et l'ombre se creuse d'elle-même dans la géométrie. Le ruban lumineux est bâti de quads indépendants pour la même raison — une bande partagée ressouderait le trou que les lignes viennent de laisser.

Rien ici n'est reconstruit à chaque image : 720 emplacements d'échantillons, 2 880 sommets de lignes et 1 440 quads de ruban sont alloués une fois, et setDrawRange fait le travail. Chaque front est dessiné deux fois, en miroir à ±delta, parce qu'un front d'onde dans une Terre à symétrie sphérique est une surface de révolution autour de l'axe de la source et que la coupe en expose deux demi-plans. Les fronts tiennent 120 s de temps simulé et s'estompent en vieillissant — la P du presque-blanc vers la braise, la S du vert d'eau pâle vers le sarcelle — pour qu'un rejeu chargé se lise comme une séquence et non comme un flou.

Une Terre en coupe faite de géométrie, pas de découpe

Pour regarder un front d'onde traverser l'intérieur, il faut ouvrir la planète, et la voie évidente — les plans de découpe — est le mauvais outil pour les coques : du travail de fragment sur chaque pixel de chaque couche, et des faces creuses au bout du compte. Les coques sont donc construites incomplètes. Une SphereGeometry avec phiStart = π et phiLength = 1.5π est une coque de 270° amputée d'un quadrant, orientée vers la caméra par défaut. Le trou est dans le maillage : il ne coûte rien.

Une couche ne peut pas user de cette astuce. Les littoraux et le graticule doivent tourner — c'est tout leur intérêt — et ils sont donc la seule chose de la scène véritablement découpée par le GPU, par deux plans avec clipIntersection activé. Les données de littoral sont décodées depuis le module empaqueté pour Orrery, puis ré-interpolées sphériquement le long de grands cercles à environ 2° par subdivision, pour que les longues cordes ne s'affaissent jamais sous la surface et ne disparaissent pas dans la sphère.

La scène tourne ensuite sur deux référentiels, et c'est la pièce dont nous sommes discrètement le plus fiers. Le référentiel de coupe ne bouge jamais : coin, strates, filets de frontière et axe polaire sont statiques, si bien que les plans de découpe en espace monde sont des constantes et que le décor des faces de coupe se cuit en un seul maillage indexé et un seul appel de dessin. Le référentiel géographique — littoraux, graticule, tout le catalogue — est amené par interpolation sphérique pour que l'hypocentre sélectionné chevauche l'axe polaire +Y, l'arête commune aux deux faces exposées. Tout plan passant par l'axe de la source est une section de grand cercle valide d'une Terre à symétrie sphérique, et il se trouve qu'ici deux d'entre eux sont des murs visibles. La coupe est honnête quel que soit l'événement choisi — et instantanée, sans la moindre animation, pour les visiteurs qui ont demandé une réduction des animations.

Quatre mille séismes dans un seul maillage

Le catalogue arrive de deux flux de résumé de l'USGS — la journée écoulée à toute magnitude, plus M4,5 et au-dessus sur le dernier mois — via un proxy qui réduit chaque événement à un tuple de position et gère son propre cache. Une revalidation au niveau de la route mettrait en cache un 503 aussi volontiers que des données, ce qui, sur Orrery, a figé un jour un état hors ligne six heures durant. Ici, le TTL de 90 secondes est ancré à l'horodatage du plus ancien des deux flux, si bien qu'un flux en échec continue d'être retenté pendant que le flux sain continue de servir. Périmé vaut mieux que vide.

À l'écran, tout le catalogue — jusqu'à 4 000 événements — est un unique InstancedMesh d'octaèdres. Chaque instance se place au rayon (6371 − depth) / 1000 en unités de scène : à sa profondeur hypocentrale réelle plutôt qu'épinglée en surface, si bien qu'une zone de subduction se lit comme une dalle plongeant dans le manteau. L'échelle est logarithmique en magnitude, 0.012 + 0.02 · 1.9^(mag − 4), bornée des deux côtés — plafonnée pour qu'un M8 reste une perle et non une lune, avec un plancher pour que les micro-séismes existent encore. La couleur, c'est l'âge : de l'orange vif à la braise sur les 24 premières heures, puis de la braise à l'ardoise au fil d'une traîne de 30 jours. L'histoire de la semaine écoulée est visible avant même le premier clic.

La sonification des séismes, par un décodeur miniSEED écrit à la main

Les sismogrammes sont publics. N'importe qui peut tirer du mouvement du sol brut des services FDSN dataselect, et les stations IU du Global Seismographic Network enregistrent en continu depuis des décennies. La dernière chose que Tremor fait d'un événement, c'est donc de proposer de vous le jouer — de la vraie sonification de séisme, pas un synthétiseur qui l'imite.

Le navigateur ne peut pas le faire lui-même, pour deux raisons indépendantes. Les données arrivent en miniSEED compressé en Steim, un format binaire à codage différentiel dont aucun navigateur n'a jamais entendu parler, et le service n'envoie aucun en-tête CORS, si bien qu'une page ne pourrait même pas récupérer les octets quand bien même elle saurait les lire. Les deux moitiés vivent sur notre serveur : un proxy récupère, un décodeur écrit à la main déballe, et un tampon Int16 normalisé parvient au client.

Choisir une station est un petit problème de classement doublé d'une contrainte physique. Tremor connaît 31 stations large bande IU et veut quarante minutes d'enregistrement sur la composante verticale, mais la plus proche n'est pas la meilleure : en deçà de 2°, l'instrument sature sur un gros événement ; au-delà de 90°, le noyau commence à faire de l'ombre à la P. La bande téléséismique mène donc le classement, la distance départage, trois candidates sont essayées dans l'ordre, et un enregistrement de moins de 45 minutes est refusé d'emblée — les archives ont besoin de temps, et le client anticipe ce refus plutôt que de vous le laisser découvrir.

Le décodeur miniSEED est à peu près ce que le code web fait de plus bas niveau : un DataView, des lectures gros-boutistes, et un en-tête blockette-1000 qui indique l'encodage en présence. Steim empaquette les échantillons dans des trames de 64 octets faites de seize mots de 32 bits, dont le premier est une carte de quartets — seize codes de deux bits décrivant l'empaquetage des quinze autres. Steim-1 offre une, deux ou quatre différences par mot ; Steim-2 ajoute des sous-codes logeant cinq différences de 6 bits, six de 5 bits ou sept de 4 bits dans un seul mot, plus un emplacement de 30 bits pour les valeurs aberrantes. Les mots un et deux de la trame zéro ne sont pas des différences du tout, mais les constantes d'intégration X0 et Xn, ce qui rend le flux auto-vérifiant. Les deux encodages sont implémentés ; là où la synchronisation des enregistrements dérive de plus de 1,5 échantillon, nous comptons le trou et concaténons quand même, parce qu'en sonification la continuité l'emporte sur l'exactitude.

Ce qui nous amène à la partie sans code. L'énergie sismique vit autour de 0,05–10 Hz ; l'audition commence vers 20 Hz. Le remède du manuel est un décaleur de hauteur, et nous n'en avons écrit aucun, parce que l'API Web Audio permet de déclarer la fréquence d'échantillonnage d'un tampon. Quarante échantillons par seconde, déclarés comme 40 × 150 = 6 000 Hz, se lisent 150 fois plus vite et 150 fois plus haut — dans la fenêtre légale des fréquences de tampon, de 3 000 à 768 000 Hz — ce qui pose la bande sismique à 7,5–1500 Hz. Quarante minutes deviennent seize secondes. Ce n'est pas une transformation des données : ce sont les données, lues vite.

La dynamique est le problème le plus dur. Une arrivée P téléséismique peut valoir cent fois la coda qui la suit, et jouée linéairement, le visiteur entend un tic puis le silence. Le client applique donc une compression de dynamique d'exposant proche de 0,6 — dans le navigateur, pas sur le serveur, parce que le travail du serveur est de compter honnêtement et que le façonnage appartient là où le visiteur peut l'entendre s'appliquer. Avant transport, l'enregistrement est démoyenné (SNZO, par exemple, se tient à quelque 1 500 comptes du zéro) et normalisé au pic. La sparkline trace le pic par casier plutôt que la valeur efficace, pour la même raison : l'arrivée est un pic, et la moyenne est justement l'opération qui la masquerait.

Le lecteur lui-même est délibérément petit : un GainNode, des rampes de 30 ms planifiées une à l'avance, un contexte audio créé à l'intérieur du gestionnaire de clic pour que la politique d'autoplay n'ait rien à redire, une lecture qui s'arrête quand l'onglet s'en va. La sourdine n'appartient pas à Tremor : il lit la sourdine du synthétiseur, partagée à l'échelle du site avec tous les jeux du Labo. Coupez le son d'un jeu, et la planète se tait aussi.

Une pastille dans le coin de l'écran annonce temps de trajet de qualité visuelle, pour la même raison que celle d'Orrery annonce propagation de qualité visuelle. Un modèle 1-D sans corrections de station et sans structure tridimensionnelle place un front d'onde là où l'œil l'attend, pas là où un instrument l'enregistrerait. Fidèle à l'œil, pas au mètre ; ne vous en servez pour rien d'autre que l'émerveillement. Le dire franchement ne coûte rien et achète au reste de l'image toute sa crédibilité.

Les manques, en toute honnêteté : le tracé de rayons tourne sur le thread principal, et la page n'a pas encore de couverture de tests de bout en bout — seulement les assertions de développement, qui vérifient la physique et non les boutons. Les deux sont sur la liste.