Ssu-Ting Lai
email: lai@ipgp.fr

Nobuaki Fuji

Ikuo Katayama

Luis Fabian Bonilla

Yann Capdeville

S.-T Lai

Rock Deformation Monitoring Using Monte Carlo Waveform Inversion

Keywords: modélisation numérique, propagation des ondes, mécanique des roches, expérience en laboratoire, nucléation sismique, supervision, méthode des éléments spectraux, inversion de forme numerical modelling, wave propagation, rock mechanics, laboratory experiment, earthquake nucleation, monitoring, spectral element method, waveform inversion, rock deformation numerical modelling, wave propagation, rock mechanics, laboratory experiment, earthquake nucleation, monitoring, spectral element method, waveform inversion, rock deformation

d'onde, déformation des roches Time-lapse rock deformation monitoring using waveform inversion

Acknowledgement

First and foremost, I would like to thank my advisor, Nobuaki Fuji, for providing all kinds of support, encouragement and good ideas during these three and a half years.

I am very thankful for the freedom he has given me to pursue independent work and to explore on my own. I thank my co-supervisor, Fabian Bonilla. He always gave me good direction, support and plenty of advice. He guided me to appreciate the importance of finding and thinking about scientific questions. I am very grateful to my two supervisors. They gave me all the guidance when I was lost in my academic career or lost confidence. Their patience and warm words supported me in completing my Ph.D. Alexandre Schubnel's insightful comments and constructive criticisms were thought-provoking, and he helped me discover and verify some ideas on rock experiments. I wish to thank Yann Capdeville for sharing his theoretical knowledge during these years. The discussions helped me understand and solve many theoretical and technical issues in my work. I am very much grateful to Guilhem Barruol, Céline Gelis and Stéphanie Durand, for their help and guidance in my comité de thèse. I would also like to thank Ikuo Katayama, Harsha Bhat, and Xiaoping Jia for their encouragement, good suggestions and discussions. A thought of gratitude goes to École Doctoral de l'IPGP and the Make Our Planet Great Again (MOPGA) project, for giving to me the opportunity of this PhD degree.

I would like to thank the members of my PhD defense committee, Alexandrine Gesret, Giuliana Rossi, Pascal Bernard, Stéphanie Durand and Lucile Bruha for their comments and suggestions on this dissertation. Many thanks to the colleagues and friends of the Seismology group. All of them helped me a lot and they made these three years and a half of research very enjoyable. I love our tradition of daily coffee break and chocolate break so much. Thanks to Sylvie Contamina for all her assistance along the way and her help for the numerous missions. Thanks to the technical staff of the IPGP Seismology group. And, special thanks to Geneviève Moguilny and Ludovic Aubry for their precious help during these years. Every Friday night gathering is what keeps me going. Furthermore, I would like to thank my family and my friends, who are now scattered around the world. My friends and family have given me endless love and strength to get this far. I am especially thankful to Zhikai, Pauline, Nassima, Shipra, Marina, Aline, Pablo, Sung-Ping, Kuan-Yu, Natalie, Pei-Shi, Kuang-Yu, Hwa-Ting, Hsiang-Hui, Cynthia, Yu-Pei, Yu-Ping and my TaiGE partners for their constant support and for being always present during the countless times when I needed it. Thanks, Egon Charles for being always ready to listen to my complaints and relieve my stress during my hardest time before the PhD defense.

I will never be able to find the proper words to express my gratitude to my dearest mother, Shu-Chuan Liang. Thank you for providing me with the best education, opportunities and talents. Thanks for teaching me that the world is large and worth exploring. You opened my horizons and always encouraged me not to be afraid to step out of my comfort zone and to follow my dreams. In my toughest times, you always accompany me to talk and unravel the problems in my heart, no matter where you are. Thanks for being my first wonderful mentor and my root, always inspiring, encouraging and supporting me. I must share this achievement with you.

Imagerie time-lapse de la déformation d'une roche à l'aide de l'inversion des formes d'ondes

Résumé

La supervision microscopique et macroscopique de l'évolution de la rupture est un élément clé pour nous aider à comprendre l'état des contraintes et les caractéristiques sismiques. À l'échelle microscopique, la mesure des propriétés des roches sous déformation en laboratoire nous fournit une vue globale des caractéristiques dynamiques et statistiques d'une roche. À l'aide ainsi de l'observation macroscopique, les expériences de déformation des roches à l'échelle du laboratoire nous ont apporté les connaissances essentielles pour comprendre la physique de la nucléation des séismes. La compréhension de l'évolution des fissures est également importante dans l'industrie pour surveiller le stockage des hydrocarbures et du CO 2 , puisque leur injection souterraine dans les réservoirs déclenche des ruptures de failles sismiques. Par conséquent, l'imagerie des scénarios de stockage time-lapse, l'étude de la sensibilité de la réponse sismique et l'investigation du processus de rupture devienneront indispensables dans l'industrie et la physique des séismes. Cette thèse propose une méthodologie d'imagerie sismique de l'état de l'art utilisant l'équation des ondes pour caractériser les changements structurels pendant la déformation de la roche. Nous avons donc développé une stratégie d'imagerie 0D et 3D time-lapse (4D) des données des expériences sismiques-actives à l'échelle du laboratoire, en proposant une méthode d'inversion pour imager efficacement les changements structurels élastiques/anélastiques pendant la déformation. Nous combinons d'abord les données observées pendant des expériences statiques et le données synthétiques pour calibrer les signaux enregistrés. Puis, nous déterminons les propriétés des ondes qui se propagent dans les échantillons. Ensuite, nous obtenons des informations sur l'hétérogénéité et l'anisotropie 3D des roches en comparant les formes d'onde observées et les synthétiques. Enfin, les méthodes d'inversion de formes d'onde et d'imagerie de migration pour les sources actives sont mises en oeuvre sur les données réelles afin de reconstruire les changements de la structure localisée. La méthode d'inversion que nous proposons a un potentiel plus élevé que les méthodes classiques de laboratoire qui utilisent uniquement les temps d'arrivée.

Abstract

Microscopic and macroscopic monitoring of crack evolution is a key to helping us understand the stress statement and seismic features of the failures. Microscopically, the measurement of rock properties under deformation in the laboratory provides us with a comprehensive view of a rock's dynamic and statistical characteristics. Together with macroscopic observation, laboratory-scale experiments of rock deformation have brought us the essential knowledge to understand the physics of earthquake nucleation. Understanding crack evolution is also critical in the industry for monitoring hydrocarbon and CO 2 storage since their underground fluid injection in the reservoirs is known to trigger seismic fault ruptures. Therefore, imaging the timelapse storage scenarios, studying the sensitivity of the varying seismic response and investigating the rupture process is becoming indispensable and extremely important in industry and earthquake physics. This thesis proposes a state-of-the-art time-lapse wave-equation-based seismic imaging methodology to infer structural changes during rock deformation. We thus developed 0D time-lapse and 3D time-lapse (4D)

active-source seismic imaging on rock samples in the laboratory, proposing an inversion method to efficiently image the elastic/anelastic structural changes during the deformation.

We first combine laboratory data of static experiments and numerical modelling to calibrate the recorded signals and determine the properties of waves propagating in the samples. Second, we obtain 3D heterogeneity and anisotropy information on rocks by comparing the observed and synthetic waveforms. Finally, waveform inversion and migration imaging methods for active sources are implemented on waveforms recorded in rock deformation experiments to reconstruct localised structural changes. Our proposed inversion method has a higher potential than the classic laboratory methods using only travel times.

Résumé

La supervision microscopique et macroscopique de l'évolution de la rupture est un élément clé pour nous aider à comprendre l'état des contraintes et les caractéristiques sismiques. À l'échelle microscopique, la mesure des propriétés des roches sous déformation en laboratoire nous fournit une vue globale des caractéristiques dynamiques et statistiques d'une roche. À l'aide ainsi de l'observation macroscopique, les expériences de déformation des roches à l'échelle du laboratoire nous ont apporté les connaissances essentielles pour comprendre la physique de la nucléation des séismes. La compréhension de l'évolution des fissures est également importante dans l'industrie pour surveiller le stockage des hydrocarbures et du CO 2 , puisque leur injection souterraine dans les réservoirs déclenche des ruptures de failles sismiques.

Par conséquent, l'imagerie des scénarios de stockage time-lapse, l'étude de la sensibilité de la réponse sismique et l'investigation du processus de rupture devienneront indispensables dans l'industrie et la physique des séismes. Cette thèse propose une méthodologie d'imagerie sismique de l'état de l'art utilisant l'équation des ondes pour caractériser les changements structurels pendant la déformation de la roche.

Nous avons donc développé une stratégie d'imagerie 0D et 3D time-lapse (4D) des données des expériences sismiques-actives à l'échelle du laboratoire, en proposant une méthode d'inversion pour imager efficacement les changements structurels élastiques/anélastiques pendant la déformation.

Résumé Long

Cette thèse vise à introduire une nouvelle technique d'imagerie dans la mécanique des roches pour observer précisement les expériences de déformation des roches à l'échelle laboratoire afin de mieux comprendre la nucléation des séismes. Pourquoi les tremblements de terre se produisent-ils ? C'est un mystère des sciences de la Terre depuis des siècles. Malgré de nombreux efforts d'exploration globale et locale de l'intérieur de la Terre, statique ou dynamique, les tremblements de terre sont toujours un mystère. [START_REF] Sengupta | Three dimensional model of seismic velocity variation in the Earth's mantle[END_REF]; [START_REF] Kreemer | A geodetic plate motion and Global Strain Rate Model[END_REF]; [START_REF] French | Broad plumes rooted at the base of the Earth's mantle beneath major hotspots[END_REF]; [START_REF] Beroza | Machine learning improves forecasts of aftershock locations[END_REF]; [START_REF] Noda | Energy-Based Scenarios for Great Thrust-Type Earthquakes in the Nankai Trough Subduction Zone, Southwest Japan, Using an Interseismic Slip-Deficit Model[END_REF], quels que soient les efforts in situ et numériques déployés pour comprendre la dynamique des ruptures sismiques (e.g. [START_REF] Bhat | Role of fault branches in earthquake rupture dynamics[END_REF][START_REF] Passelègue | Dynamic rupture processes inferred from laboratory microearthquakes[END_REF], la prédiction des séismes est intrinsèquement impossible (e.g. [START_REF] Geller | Shake-up time for Japanese seismology[END_REF]. Pourtant, une meilleure compréhension de la physique des tremblements de terre permettra de réduire les pertes de vies humaines et de biens que les tremblements de terre causent.

La sismologie des tremblements de terre a été maintenue à jour grâce à des études continues. Il convient de noter que les scientifiques ont eu un débat sur la poule et l'oeuf à propos d'un tremblement de terre et de ses failles sismiques jusqu'à la fin du XIXe siècle. Par exemple, on pensait qu'une dislocation soudaine des failles était à l'origine des grands tremblements de terre à une époque antérieure (e.g. [START_REF] Mallet | On the dynamics of earthquakes; being an attempt to reduce their observed phenomena to the known laws of wave motion in solids and fluids[END_REF][START_REF] Koto | On the case of the great earthquake in central Japan, 1891[END_REF], alors qu'une théorie quantitative de séisme a été introduite après le tremblement de terre de San Francisco de 1906 [START_REF] Lawson | The California Earthquake of April 18[END_REF]Reid, 1910, e.g.). En comparant les observations géologiques et géodésiques de ce séisme avec les autres grands séismes sur la dislocation des failles sous contrainte, [START_REF] Reid | The mechanics of the earthquake[END_REF] a mis en avant la théorie du rebond élastique et a conduit aux fondements généralement acceptés des mécanismes de faille sur l'origine des séismes : i) la nucléation des tremblements de terre est un processus d'accumulation, où la déformation aux discontinuités de la roche augmente jusqu'à ce que la roche ne puisse plus la supporter, provoquant finalement une rupture ; ii) la rupture commence à la surface de la discontinuité et s'étend vers le bas, le rebond élastique relatif de la roche déformée se produit sous ses contraintes jusqu'à ce que la déformation soit soulagée ; iii) les déplacements finis pendant l'événement sont confinés à une zone de quelques kilomètres de failles et se propagent à une vitesse inférieure à la vitesse des ondes de compression de la roche. Ces affirmations sont la base de la compréhension de la relation entre les failles et les séismes et ont été utilisées par la suite dans toutes les études sur l'origine des séismes. [START_REF] Nakano | Note on the nature of forces which give rise to the earthquake motions[END_REF] a étudié théoriquement la propagation des ondes élastiques qui sont générées par différents types de forces exercées en un point d'un solide élastique infini, ce qui fournit une base théorique pour déterminer la position initiale de la dislocation de la faille. [START_REF] Steketee | Some Geophysical Applications of the Elasticity Theory of Dislocations[END_REF] a d'abord suggéré que les dislocations dans un demi-espace élastique pouvaient être utilisées pour modéliser la déformation due à la formation de failles. Puis, il a réalisé des applications des noyaux de déformation en utilisant la formule de Volterra, qui donne le champ de déplacement comme une intégrale sur la surface de la faille impliquant la nucléation de la déformation [START_REF] Mansinha | The displacement fields of inclined faults[END_REF]. Des sources cinématiques, y compris dans un modèle de faille rectangulaire simple, ont été proposées successivement par [START_REF] Haskell | Total energy and energy spectral density of elastic wave radiation from propagating faults[END_REF][START_REF] Haskell | Total energy and energy spectral density of elastic wave radiation from propagating faults. Part II. A statistical source model[END_REF] pour décrire le double couple équivalent à la dislocation de la faille. Les modèles de propagation des fissures, qui décrivent les caractéristiques essentielles des fissures au sein d'une faille circulaire [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], d'une faille dip-glissement [START_REF] Starr | Slip in a crystal and rupture in a solid due to shear[END_REF] et d'une faille à glissement latéral [START_REF] Knopoff | Energy Release in Earthquakes[END_REF], du mécanisme de fracture, ont été développés en introduisant une analogie entre les séismes et les fissures. Le mécanisme des roches a ensuite été introduit dans la mécanique sismique par le biais d'expériences sur la friction des roches. [START_REF] Brace | Stick-Slip as a Mechanism for Earthquakes[END_REF]; [START_REF] Byerlee | The mechanics of stick-slip[END_REF] ont mesuré la force de friction des roches, suggérant que la contrainte de cisaillement sur les failles était élevée. Ils ont réalisé que les phénomènes de friction par glissement observés en laboratoire pouvaient être utilisés pour comprendre les ruptures survenant lors de grands séismes.

Les expériences qu'ils ont menées à l'aide d'échantillons de roches découpés à la scie, soumis à une contrainte aux deux extrémités et à une pression de confinement, ont permis de constater que le glissement entre les deux morceaux de roches n'était pas lisse, mais qu'il s'effectuait plutôt selon un cycle de stick-slip, c'est-à-dire que les surfaces se bloquent momentanément, puis se relâchent soudainement et glissent.

Il s'agissait de la première théorie bien acceptée des séismes de type stick-slip sur des failles préexistantes. La mesure des propriétés des roches en laboratoire nous fournit des vues microscopiques et statistiques sur les caractéristiques dynamiques et/ou statistiques d'une roche (e.g. Lockner et al., 1977;Toksöz et al., 1979). Des expériences en laboratoire sur la propagation des ondes dans des échantillons de roche sont généralement réalisées pour déduire la relation entre les changements des conditions externes et les caractéristiques sismiques de la roche, telles que la vitesse, l'atténuation et l'anisotropie, ce qui nous aide à mieux comprendre le processus d'affaiblissement et de fracturation de la roche. L'évolution de la tectonique des plaques à l'échelle mondiale a également permis de mieux comprendre la mécanique des failles. La dislocation des failles à la marge des plaques, qui provoque les tremblements de terre, a été expliquée par la théorie de la tectonique des plaques qui a également découvert l'origine des déformations élastiques de la croûte.

La mécanique de la rupture a été étudiée sous différents aspects et méthodes. La structure des zones de faille (figure 1.1) est fortement liée au processus de rupture des séismes et peut nous donner un aperçu des conséquences des séismes. Elle dépend de la profondeur de la formation, du protolithe, de l'environnement tectonique (par exemple glissement, extension ou compression), de l'ampleur du déplacement et de l'écoulement des fluides [START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF]. Le modèle conceptuel de la structure des zones de faille a évolué et s'est développé au cours des dernières décennies, passant d'un noyau de faille unique [START_REF] Chester | Implications for mechanical properties of brittle faults from observations of the Punchbowl fault zone, California[END_REF]) à un noyau de faille avec plusieurs branches compliquées, anastomosées et reliées, entraînant des blocs ou des lentilles de protolithe fracturé entre les couches [START_REF] Faulkner | On the internal structure and mechanics of large strike-slip fault zones: Field observations of the Carboneras fault in southeastern Spain[END_REF]. Les chercheurs ont mis en oeuvre des paramètres plus simples en laboratoire que dans la réalité pour comprendre le processus de fracture sous pression afin d'imiter les conditions réelles des tremblements de terre d'une manière contrôlable et de relier d'autres expériences de déformation des roches et la déformation de la croûte (e.g. [START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF]). On considère que les roches fragiles subissent une certaine initiation et propagation des fractures avant leur rupture (e.g. [START_REF] Hoek | Brittle fracture propagation in rock under compression[END_REF] Van der [START_REF] Van Der Pluijm | Earth structure. An introduction to structural geology and tectonics[END_REF]. (figure 1.2 et 1.3), où la microfissure est généralement considérée comme ayant un rôle clé dans les noyaux de fracture. Le processus de fracture peut être divisé en quatre étapes dans une expérience de compression confinée avec un tracé de contrainte-déformation correspondant (Van der [START_REF] Van Der Pluijm | Earth structure. An introduction to structural geology and tectonics[END_REF]. Un plan de faille se produira en raison de la connexion de microfissures en un plan continu, alors que le plan se formera généralement à environ 30°de la direction de la contrainte principale. Puisque les vitesses moyennes des ondes élastiques d'un échantillon de roche chargé peuvent refléter la nucléation et/ou la croissance localisée des fissures, de nombreuses études ultérieures ont porté sur les changements temporels d'une roche expérimentale pendant sa déformation afin d'obtenir des informations sur le processus macroscopique de génération des tremblements de terre (Scholz et al., 1973).

Grâce à ces avancées scientifiques dans l'ensemble de la communauté des sciences de la Terre, à des observations précises et denses et à des capacités numériques, les points de vue modernes sur les tremblements de terre se sont enrichis et diversifiés.

Afin de comprendre les comportements de la rupture dynamique de la croûte terrestre et du manteau lors de la nucléation des séismes, des recherches macroscopiques et microscopiques ont été menées dans le cadre de la modélisation numérique (e.

g. Fliss et al., 2005;Gabriel et al., 2012), des expériences à l'échelle du laboratoire (e.g. Guéguen and Schubnel, 2003;Benson et al., 2006) et de l'observation de séismes naturels (e.g. Di Carli et al., 2010;Lee et al., 2011). Cette thèse tente de proposer une nouvelle méthodologie sur des mesures à l'échelle du laboratoire en utilisant des techniques d'imagerie sismique de pointe appliquées à la tomographie globale et régionale autrement. Les études classiques des expériences de déforma-tion des roches ont été menées en utilisant une simple sélection manuelle des temps d'arrivée des sources actives (Lockner et al., 1977;Watanabe and Higuchi, 2015;Zaima and Katayama, 2018, e.g.). En revanche, la technique de l'imagerie sismique a été beaucoup plus développée depuis la première génération de tomographie globale (e.g. [START_REF] French | Broad plumes rooted at the base of the Earth's mantle beneath major hotspots[END_REF][START_REF] Konishi | Waveform inversion for localized three-dimensional seismic velocity structure in the lowermost mantle beneath the Western Pacific[END_REF][START_REF] Yuan | High-frequency localized elastic fullwaveform inversion for time-lapse seismic surveys[END_REF]Fuji et al., 2021). Nous voulons appliquer ces méthodologies d'imagerie sismique de pointe à un volume réduit de données dans un petit échantillon de roche. Avant cette étude, plusieurs efforts ont été faits pour imager sismologiquement l'intérieur d'échantillons de roche. Möllhoff et al. (2010) a simulé la fracture artificielle préexistante dans la roche en utilisant la méthode du réseau élastique discret 2D, tandis que Lai et al. (2019) a étudié les changements du modèle de vitesse pendant la déformation de la roche en utilisant la méthode des différences finies 2D. Cependant, la modélisation numérique 2D est insuffisante pour décrire la propagation des ondes dans une roche 3D, même si la fonction de transfert de Lai et al. (2019) ou une approximation 2D de source linéaire (Igel et al., 2002) semble fonctionner dans une certaine mesure. Pour dépasser les limites, des efforts prometteurs ont été faits récemment en matière de modélisation numérique 3D des formes d'ondes sismiques à l'échelle du laboratoire, malgré certaines difficultés. Yoshimitsu et al. (2016a), par exemple, a modélisé l'effet de la géométrie sur l'aluminium cylindrique en utilisant la méthode des 3D différences finies pour un cylindre de 50 mm de diamètre et 100 mm de hauteur avec une fréquence dominante entre 200 kHz-800 kHz. Leurs comparaisons entre les données synthétiques et les données observées, limitées aux expériences sur les métaux et plastiques homogènes, ont montré la faisabilité de la reproductibilité numérique des données sismiques actives en laboratoire. En plus de la simulation de la propagation des ondes, Brantut (2018a) a suivi l'émission acoustique sur du grès à petite échelle en utilisant une méthode 3D de marche rapide avec la méthode des différences finies et a vérifié concomitamment l'évolution des changements de vitesse de la roche pendant sa déformation. [START_REF] Sarout | Stressdependent permeability and wave dispersion in tight cracked rocks: Experimental validation of simple effective medium models[END_REF] ont estimé l'évolution des paramètres des fissures en inversant les ondes élastiques et les données de perméabil-ité. Zotz-Wilson et al. (2020) ont imagé le développement des fractures en utilisant la migration des données sismiques actives, qu'ils ont appelée inversion de décorrélation des ondes coda. Les études présentées sont encourageantes, mais certains des traitements numériques du problème complexe actuel devraient être favorisés.

La condition aux limites employée dans la méthode des différences finies est irréaliste (Komatitsch and Vilotte, 1998) détecter la limite noyau-manteau de la Terre [START_REF] Gutenberg | Über erdbebenwellen. vii a. beobachtungen an registrierungen von fernbeben in göttingen und folgerung über die konstitution des erdkörpers (mit tafel)[END_REF]. La limite noyaumanteau de Mars, en revanche, n'a été détectée que deux ans après le déploiement d'une seule station InSight [START_REF] Stähler | Seismic detection of the martian core[END_REF]. Malgré la quantité limitée de données, les sismologues se sont efforcés de mettre à jour et d'affiner la structure intérieure de la Terre globale et locale et des autres planètes en faisant progresser les efforts théoriques et numériques pour extraire le plus d'informations possible.

La sismologie globale et d'exploration a d'abord utilisé les temps de parcours des ondes de corps ou la dispersion de la vitesse de phase des ondes de surface pour effectuer la tomographie sismique ; mais plus tard, des techniques d'inversion de forme d'onde ont été établies pour traiter les caractéristiques du train d'onde entier (e.g. [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF][START_REF] Geller | Two efficient algorithms for iterative linearized inversion of seismic waveform data[END_REF]. En profitant de l'abondance et de la qualité croissantes des réseaux sismiques globaux/régionaux/locaux ainsi que des développements théoriques, l'inversion de forme d'onde révèle des caractéris-tiques détaillées de l'intérieur de la Terre à diverses résolutions (e.g. [START_REF] French | Broad plumes rooted at the base of the Earth's mantle beneath major hotspots[END_REF][START_REF] Konishi | Waveform inversion for localized three-dimensional seismic velocity structure in the lowermost mantle beneath the Western Pacific[END_REF]Fuji et al., 2021). Grâce à ces pratiques, les sismologues comprennent maintenant l'équilibre entre le rang des informations indépendantes au sein des données observées, et le degré et le type d'informations que la communauté doit acquérir, ce qui encourage les travaux interdisciplinaires qui se couplent avec la géodynamique et la géochimie (e.g. [START_REF] Schuberth | Traveltime dispersion in an isotropic elastic mantle: strong lower-mantle signal in differential-frequency residuals[END_REF][START_REF] Franken | Seismic wave-based constraints on geodynamical processes: An application to partial melting beneath the réunion island[END_REF][START_REF] Fuji | Small-data seismology-from a centimetric rock sample to Mars[END_REF].

La surveillance time-lapse vise à estimer les changements des propriétés du milieu en fonction du temps, où la différence entre les données sismiques mesurées à différents pas de temps devrait être causée par le changement des propriétés du milieu, ce qui se traduit par la structure sismique. L'imagerie time-lapse est également cruciale dans l'industrie et appliquée à la surveillance du stockage des hydrocarbures et du CO 2 , à la gestion des réservoirs et à la surveillance à long terme des changements souterrains d'un état de contrainte et des caractéristiques sismiques (e.g. [START_REF] Arts | Monitoring of CO2 injected at Sleipner using time-lapse seismic data[END_REF][START_REF] Picotti | Seismic modeling to monitor CO$_{2}$ geological storage: The Atzbach-Schwanenstadt gas field[END_REF][START_REF] Riahi | Time-lapse analysis of ambient surface wave anisotropy: A three-component array study above an underground gas storage[END_REF][START_REF] Planès | Time-lapse monitoring of internal erosion in earthen dams and levees using ambient seismic noise[END_REF][START_REF] Stork | Assessing the potential to use repeated ambient noise seismic tomography to detect CO 2 leaks: Application to the Aquistore storage site[END_REF][START_REF] Cotton | Time-lapse velocity analysis -Application to onshore continuous reservoir monitoring[END_REF][START_REF] Solazzi | Seismic Signatures of Fractured Porous Rocks: The Partially Saturated Case[END_REF][START_REF] Dupuy | Combined geophysical and rock physics workflow for quantitative CO2 monitoring[END_REF][START_REF] Yuan | High-frequency localized elastic fullwaveform inversion for time-lapse seismic surveys[END_REF]. Arts Outre la compréhension de la physique des séismes, les effets physiques des roches sont largement reconnus comme jouant un rôle clé dans l'interprétation quantitative des données sismiques dans le domaine du stockage du CO 2 . Par conséquent, de nombreuses études liées à la pétrophysique ont essayé de se connecter avec la recherche sur la surveillance sismique et vice versa [START_REF] Vanorio | The rock physics basis for 4D seismic monitoring of CO2 fate: Are we there yet?[END_REF][START_REF] Picotti | Seismic modeling to monitor CO$_{2}$ geological storage: The Atzbach-Schwanenstadt gas field[END_REF][START_REF] Njiekak | CO2 rock physics as part of the Weyburn-Midale geological storage project[END_REF][START_REF] Prasad | CO2 messes with rock physics[END_REF]. En considérant différents modèles de roche et en analysant la sensibilité de la réponse sismique avec des saturations variables de CO 2 , [START_REF] Picotti | Seismic modeling to monitor CO$_{2}$ geological storage: The Atzbach-Schwanenstadt gas field[END_REF] a construit un modèle pétro-élastique pour détecter les fuites potentielles dans le réservoir. [START_REF] Njiekak | CO2 rock physics as part of the Weyburn-Midale geological storage project[END_REF] a modélisé et analysé les formes d'ondes ultrasoniques enregistrées lors d'expériences de saturation en CO 2 sur des roches. La chute des vitesses des ondes P et S s'est avérée marquer la phase CO 2 . [START_REF] Prasad | CO2 messes with rock physics[END_REF] a construit un système pour les propriétés acoustiques de l'interaction CO 2 -eau et a appliqué des micrographies time-lapse pour surveiller les effets de l'interaction eau de mer-grès et l'interaction entre l'eau riche en CO 2 et le carbonate. En combinant les connaissances de la physique des roches, des fluides et des propriétés des ondes sismiques, ils ont étudié les changements des ondes sismiques se propageant à l'intérieur des roches en raison de l'injection de CO 2 , qui s'est avérée être loin d'être un pur problème de substitution de fluide.

La relation temporelle entre les propriétés physiques des roches et les ondes sismiques permet de révéler le processus de nucléation des séismes et est ensuite largement utilisée dans l'industrie. Il est donc nécessaire de développer des modèles numériques et de mener des expériences en laboratoire pour déterminer avec pré-cision les caractéristiques des roches dans des conditions de stress. Les mesures à l'échelle du laboratoire à l'aide de méthodes sismiques, ou la surveillance sismique en temps réel elle-même, sont souvent difficiles. Le manque d'études menant des expériences numériques et de laboratoire sur des spécimens relativement petits dans des conditions bien contrôlées peut entraîner des divergences dans l'interprétation des données expérimentales. Même si certaines études considèrent, par exemple, l'impact de la taille des transducteurs ultrasoniques (Kovalyshen et al., 2020), l'effet du transducteur, l'effet de la géométrie [START_REF] Yoshimitsu | Geometric effect on a laboratory-scale wavefield inferred from a three-dimensional numerical simulation[END_REF], qui sont nécessaires pour obtenir des informations adéquates dans les formes d'onde, il y a eu peu d'études qui ont analysé quantitativement les formes d'onde sismiques actives timelapse pendant les expériences de laboratoire (e.g. [START_REF] Lai | Rock deformation monitoring using Monte Carlo waveform inversion[END_REF]. Par conséquent, dans cette étude, nous voulons combiner systématiquement ensemble dans ce projet les expériences de laboratoire, la modélisation numérique et la surveillance sismique time-lapse.

La roche est, par nature, un mélange complexe de minuscules grains qui s'alignent de manière aléatoire ou systématique. Les ondes sismiques sont sensibles aux différentes moyennes spatiales du milieu qu'elles traversent. Dans les milieux hétérogènes tels que les roches, la dispersion des ondes sismiques peut conduire à des caractéristiques significatives des ondes de Coda. L'interférométrie des ondes est couramment utilisée pour surveiller les changements mineurs des ondes de Coda, qui sont des ondes multidiffusées enregistrées après l'arrivée des ondes de corps [START_REF] Singh | Coda Wave Interferometry for Accurate Simultaneous Monitoring of Velocity and Acoustic Source Locations in Experimental Rock Physics[END_REF][START_REF] Obermann | 4-D Imaging of Subsurface Changes with Coda Waves: Numerical Studies of 3-D Combined Sensitivity Kernels and Applications to the $$M_{\mathrm{w}}$$7.9, 2008 Wenchuan Earthquake[END_REF]Zotz-Wilson et al., 2020). Cependant, la plupart des études de surveillance des faibles changements dans la subsurface terrestre à l'aide d'ondes multi-dispersées/coda se sont limitées à l'identification des faibles changements plutôt qu'à la compréhension de la physique de ces changements. Par rapport à la méthode standard qui consiste à mesurer les changements de vitesse en choisissant les temps de première arrivée ou l'analyse des ondes de Coda, nous déployons une analyse de forme d'onde complète pour extraire les informations physiques sur les changements structurels. étudiées (e.g. [START_REF] Sivaji | A physicalmodel study of the statistics of seismic waveform fluctuations in random heterogeneous media[END_REF]. En appliquant l'approximation de Born, les problèmes de diffusion sismique causés par une hétérogénéité à petite échelle et aléatoire peuvent être résolus [START_REF] Hudson | The use of the Born approximation in seismic scattering problems[END_REF][START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition[END_REF]. Cependant, aucune recherche ne considère et n'analyse systématiquement le cadre expérimental du laboratoire, l'entrée de l'ondelette source, la propagation des ondes, les phases de l'onde, l'hétérogénéité aléatoire en même temps. [START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF]: (a) a model with a single high-strain core surrounded by a fractured damage zone [START_REF] Chester | Implications for mechanical properties of brittle faults from observations of the Punchbowl fault zone, California[END_REF] and (b) a multiple-core model, where many strands of high-strain material enclose lenses of fractured protolith [START_REF] Faulkner | On the internal structure and mechanics of large strike-slip fault zones: Field observations of the Carboneras fault in southeastern Spain[END_REF] The reference is the waveform at 0 MPa compression stress state. . .

3.

12 Spectral ratio of observed S-wave waveforms of the dry experiments.

The reference is the waveform at 0 MPa compression stress state. . .

3.13

Cross-correlation tests for P-wave waveforms of the dry experiments.

3.14 Cross-correlation tests for S-wave waveforms of the dry experiments.

3.15 Cross-correlation tests time-lag matrix. . . . . . . . . . . . . . . . .

3.16

The snapshot of the P-wave wavefield at different time steps from (a)

2 to (f) 7 µs, showing P-wave propagation process inside the medium.

3.17 The snapshot of the S-wave wavefield at different time steps from (a)

2 to (f) 7 µs, showing S-wave propagation process inside the medium.

Q -1 (Q -1 /Q -1 0 ) of

Chapter1 Introduction

This thesis aims to introduce a novel imaging technology to rock mechanics to accurately observe the laboratory-scale rock deformation experiments to understand earthquake nucleation further. Why do earthquakes happen? It has been a major mystery in Earth sciences throughout centuries. Albeit numerous efforts in global and local exploration of Earth's interior, static or dynamic (e.g., [START_REF] Sengupta | Three dimensional model of seismic velocity variation in the Earth's mantle[END_REF][START_REF] Kreemer | A geodetic plate motion and Global Strain Rate Model[END_REF][START_REF] French | Broad plumes rooted at the base of the Earth's mantle beneath major hotspots[END_REF][START_REF] Beroza | Machine learning improves forecasts of aftershock locations[END_REF][START_REF] Noda | Energy-Based Scenarios for Great Thrust-Type Earthquakes in the Nankai Trough Subduction Zone, Southwest Japan, Using an Interseismic Slip-Deficit Model[END_REF], regardless of in situ and numerical efforts in understanding the seismic rupture dynamics (e.g., [START_REF] Bhat | Role of fault branches in earthquake rupture dynamics[END_REF][START_REF] Passelègue | Dynamic rupture processes inferred from laboratory microearthquakes[END_REF], earthquake prediction is inherently impossible (e.g., [START_REF] Geller | Shake-up time for Japanese seismology[END_REF]). Yet, a better understanding of the physics of earthquakes will reduce the loss of human life and property that earthquakes could cause.

Earthquake seismology has been kept updated with continuous studies. It is worth noting that scientists had had a chicken-and-egg debate on an earthquake and its seismic faults until the late nineteenth century. For instance, a sudden dislocation of the faults was thought to cause large earthquakes in an earlier time (e.g., [START_REF] Mallet | On the dynamics of earthquakes; being an attempt to reduce their observed phenomena to the known laws of wave motion in solids and fluids[END_REF][START_REF] Koto | On the case of the great earthquake in central Japan, 1891[END_REF] whereas a quantitative theory of earthquake genesis was introduced after the 1906 San Francisco earthquake [START_REF] Lawson | The California Earthquake of April 18[END_REF][START_REF] Reid | The mechanics of the earthquake[END_REF]. By comparing geological and geodetic observations of this earthquake with the other great earthquakes on fault dislocation under in situ stress, [START_REF] Reid | The mechanics of the earthquake[END_REF] put forward the theory of elastic rebound and led to generally accepted foundations of fault mechanisms on earthquakes' origin: i) earthquakes nucleation is an accumulative process, where the strain at rock discontinuities grows until the rock cannot endure it, eventually causing a rupture; ii) the rupture begins at the discontinuity surface and extends in the volume, the relative elastic rebound of the strained rock occurs under stresses until the strain is relieved; iii) finite displacements during the event are confined to a zone of a few kilometres of faults and spread at a rate less than the compressional-wave speed of the rock. These statements lay the basic understanding of the relationship between faults and earthquakes and have later been used in all the studies of earthquakes origin. [START_REF] Nakano | Note on the nature of forces which give rise to the earthquake motions[END_REF] theoretically studied the propagation of elastic waves, which are generated by various force types exerted at a point in an infinite elastic solid, which provides a theoretical basis for determining the initial fault dislocation position. [START_REF] Steketee | Some Geophysical Applications of the Elasticity Theory of Dislocations[END_REF] first suggested that dislocations in an elastic half-space could be used to model deformation due to faulting and conducted applications of strain nuclei using Volterra's formula, which gives the displacement field as an integral over the fault surface involving nucleation of strain [START_REF] Mansinha | The displacement fields of inclined faults[END_REF]. Kinematic sources, including in a simple rectangular fault model, have been proposed successively by [START_REF] Haskell | Total energy and energy spectral density of elastic wave radiation from propagating faults[END_REF][START_REF] Haskell | Total energy and energy spectral density of elastic wave radiation from propagating faults. Part II. A statistical source model[END_REF] to describe double couple equivalent to fault dislocation. Crack propagation models, which describe essential features of cracks within a circular fault [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], a dip-slip fault [START_REF] Starr | Slip in a crystal and rupture in a solid due to shear[END_REF] and a strike-slip fault [START_REF] Knopoff | Energy Release in Earthquakes[END_REF], of fracture mechanism, were developed by introducing an analogy between earthquakes and cracks. Rock mechanics was later brought into earthquake mechanics via experiments on rock friction. [START_REF] Brace | Stick-Slip as a Mechanism for Earthquakes[END_REF] and [START_REF] Byerlee | The mechanics of stick-slip[END_REF] measured the frictional strength of rocks, suggesting that shear stress on faults is high. They realised that the stick-slip friction events observed in the laboratory could be used to understand ruptures occurring during large earthquakes. Their experiments using saw-cut rock samples with stress at both ends and confining pressure revealed that the sliding between the two pieces of rocks was not smooth but rather in a cycle of stick-slips, which represent the surfaces momentarily lock and then suddenly release and slide. This was the first well-accepted [START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF]: (a) a model with a single high-strain core surrounded by a fractured damage zone [START_REF] Chester | Implications for mechanical properties of brittle faults from observations of the Punchbowl fault zone, California[END_REF] and (b) a multiple-core model, where many strands of high-strain material enclose lenses of fractured protolith [START_REF] Faulkner | On the internal structure and mechanics of large strike-slip fault zones: Field observations of the Carboneras fault in southeastern Spain[END_REF] . theory of stick-slip earthquakes on pre-existing faults. Measurement of rock properties at laboratory provides us with microscopic and statistic views on the dynamic and/or statistic characteristics of a rock (e.g., Lockner et al., 1977;Toksöz et al., 1979). Laboratory experiments of wave propagation in rock samples are usually carried out to infer the relationship between changes in external conditions and seismic characteristics of the rock, such as velocity, attenuation and anisotropy, which helps us foster our understanding of the weakening and fracturing process of the rock.

The evolution in global plate tectonics has also brought some insights into fault mechanics. Dislocation of faults at plate margin, which causes the earthquakes, was explained by the theory of plate tectonics that also solved the origin of elastic strains in the crust [START_REF] Le Pichon | Plate tectonics[END_REF].

Rupture mechanics has been studied using various aspects and methods. Fault zone structure (figure 1.1) is strongly related to the earthquake rupture process and can give us insight as to the evidence of earthquake consequence. It depends on the σ1 σ3 Figure 1.2: Schematic diagram of the role of micro cracks in fracture process: A fault plane will occur due to the connection of micro cracks into one continual plane, whereas the plane will generally form at around 30°to the direction of principle stress.

depth of formation, protolith, tectonic environment (e.g., strike-slip, extension or compressional), the magnitude of displacement and fluid flow [START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF].

The conceptual model for fault zone structure evolved and developed over the past decades from a single fault core [START_REF] Chester | Implications for mechanical properties of brittle faults from observations of the Punchbowl fault zone, California[END_REF] to fault core with plural complicated branches, anastomose and links, entraining blocks or lenses of fractured protolith between the layers [START_REF] Faulkner | On the internal structure and mechanics of large strike-slip fault zones: Field observations of the Carboneras fault in southeastern Spain[END_REF]. Researchers implemented settings in the laboratory simpler than the reality to understand the fracture process under pressure to mimic real earthquake conditions in a controllable way and to connect further rock deformation experiments and deformation of the crust (e.g., [START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF]. Brittle rocks are considered to undergo a certain fracture initiation and propagation before rock failure (e.g., [START_REF] Hoek | Brittle fracture propagation in rock under compression[END_REF][START_REF] Van Der Pluijm | Earth structure. An introduction to structural geology and tectonics[END_REF]) (figure 1.2 and 1.3), where micro crack is commonly considered as a key role in the fracture nuclei. The fracture process can be divided into four stages in a confined compression experiment with accompanying stress-strain plot (Van der [START_REF] Van Der Pluijm | Earth structure. An introduction to structural geology and tectonics[END_REF]. A fault plane is generated by the connection of micro-cracks into one continual plane, whereas the plane will generally form at around 30°to the direction of the principal stress. Since averaged elastic wave velocities of a loaded rock sample could reflect the localised nucleation and/or growth of cracks, many subsequent studies worked on time-lapse changes of an experimental rock during its deformation to get insights on the macroscopic earthquake generation process (Scholz et al., 1973).

Through these scientific advances in the whole Earth science community and precise and dense observations and numerical capacities, modern views on earthquakes have enriched and diversified. In order to understand the behaviours of dynamic rupture of the Earth's crust and mantle during the earthquake nucleation, both macroscopic and microscopic investigations have been made in numerical modelling (e.g., Fliss et al., 2005;Gabriel et al., 2012), laboratory-scaled experiments (e.g., Guéguen and Schubnel, 2003;Benson et al., 2006) and natural earthquakes observation (e.g., Di Carli et al., 2010;Lee et al., 2011). This thesis attempts to propose a novel methodol- ogy on laboratory-scaled measurements using state-of-the-art seismic imaging techniques applied to global and regional tomography. Classical studies of rock deformation experiments had been conducted using simple manual picking of arrival times of active sources (e.g., Lockner et al., 1977;Watanabe and Higuchi, 2015;Zaima and Katayama, 2018). In contrast, seismic imaging technology has been much further developed since the first generation of global tomography (e.g., [START_REF] French | Broad plumes rooted at the base of the Earth's mantle beneath major hotspots[END_REF][START_REF] Konishi | Waveform inversion for localized three-dimensional seismic velocity structure in the lowermost mantle beneath the Western Pacific[END_REF][START_REF] Yuan | High-frequency localized elastic fullwaveform inversion for time-lapse seismic surveys[END_REF]Fuji et al., 2021). We want to apply these state-of-the-art seismic imaging methodologies to a small-data regime in a small rock sample. Prior to this study, there have been several efforts to image the interior of rock samples seismologically. Möllhoff et al. (2010) simulated the preexisting artificial fracture in the rock using 2D discrete elastic lattice method, whereas Lai et al. (2019) studied the velocity model changes during rock deformation using the 2D finite differences method. However, 2D numerical modelling is insufficient to describe wave propagation in a 3D rock, even though the transfer function of Lai et al. (2019) or a 2D line source approximation (Igel et al., 2002) appear to be effective in some cases. To overcome the limits, there have been some promising recent efforts of 3D numerical modelling of seismic waveforms on a laboratory scale, albeit some difficulties. Yoshimitsu et al. (2016a), for instance, modelled the geometry effect on cylindrical aluminium using the 3D finite difference method for a cylinder with 50 mm in diameter and 100 mm in height with a dominant frequency between 200 kHz-800 kHz. Their comparisons between synthetics and observed data, limited to homogeneous metal and plastic experiments, have shown the feasibility of the numerical reproducibility of active seismic data in the laboratory. In addition to wave propagation simulation, Brantut (2018a) monitored acoustic emission on small-scaled sandstone using a 3D fast-marching method with the finite-difference method and concomitantly constructed the evolution of velocity changes of the rock during its deformation. [START_REF] Sarout | Stressdependent permeability and wave dispersion in tight cracked rocks: Experimental validation of simple effective medium models[END_REF] estimated the evolution of crack parameters by inverting elastic waves and permeability data. Zotz-Wilson et al. (2020) imaged the development of fractures using migration of active seismic data, which they called coda wave decorrelation inversion. The presented studies are encouraging, but some of the numerical treatments of the present complex problem should be fostered. The boundary condition employed in the finite difference method is unrealistic (Komatitsch and Vilotte, 1998), and the kernels used in the migration are over-simplified using generalised ray theory. The present work aims to extend the idea, combine with 3D numerical modelling on wave propagation, and carefully treat the entire waveforms to extract time-lapse information on elastic and anelastic parameters changes (figure 1.4).

Since we would like to introduce state-of-the-art seismic imaging techniques in the laboratory of rock mechanics, here we present a (very) brief overview of global and exploration seismology. Only twenty years after the first record of a teleseismic event (Tokyo-Potsdam) in 1889 (von Rebeur-Paschwitz, 1889), Gutenberg and Geiger used less than 100 traveltimes to detect the core-mantle boundary of the Earth [START_REF] Gutenberg | Über erdbebenwellen. vii a. beobachtungen an registrierungen von fernbeben in göttingen und folgerung über die konstitution des erdkörpers (mit tafel)[END_REF]. The core-mantle boundary of Mars, on the other hand, was detected only two years after the deployment of a single station InSight [START_REF] Stähler | Seismic detection of the martian core[END_REF]. Albeit the limited quantity of data, seismologists have struggled to update and refine the interior structure of the global and local Earth and other planets by advancing the theoretical and numerical efforts to extract as much information as possible. Global and exploration seismology has first used traveltimes of body waves or phase velocity dispersion of surface waves to perform seismic tomography (e.g., [START_REF] Sengupta | Three dimensional model of seismic velocity variation in the Earth's mantle[END_REF][START_REF] Montagner | Global anisotropy in the upper mantle inferred from the regionalization of phase velocities[END_REF]) but later waveform inversion techniques were established to handle the whole wavetrain features (e.g., [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF][START_REF] Geller | Two efficient algorithms for iterative linearized inversion of seismic waveform data[END_REF]. By profiting from the increasing abundance and quality of global/regional/local seismic networks along with the theoretical developments, waveform inversion reveals detailed features of the Earth's interior on various resolutions (e.g., [START_REF] French | Broad plumes rooted at the base of the Earth's mantle beneath major hotspots[END_REF][START_REF] Konishi | Waveform inversion for localized three-dimensional seismic velocity structure in the lowermost mantle beneath the Western Pacific[END_REF]Fuji et al., 2021). Throughout these practices, seismologists now understand the balance between the rank of independent information within the observed data and the degree and type of information that the community needs to acquire, which encourages The schematic schema of our ultimate time-lapse imaging strategy during rock deformation experiments. The grey rectangle represents the cross-section of the sample. t = 1∆t, 2∆t and 3∆t indicate the deformation stages that the higher ∆t represents that higher deformed condition. Our idea is to build a system that 1) process a huge database of numerical models and their corresponding synthetic waveforms; 2) match the experimental waveforms and synthetic waveforms and obtain the degree (%) of agreement between the experimental data and the simulated waveform; 3) obtain the most potential structural changes among the first few bestmatched models.

interdisciplinary works that couple with geodynamics and geochemistry, for instance (e.g., [START_REF] Schuberth | Traveltime dispersion in an isotropic elastic mantle: strong lower-mantle signal in differential-frequency residuals[END_REF][START_REF] Franken | Seismic wave-based constraints on geodynamical processes: An application to partial melting beneath the réunion island[END_REF][START_REF] Fuji | Small-data seismology-from a centimetric rock sample to Mars[END_REF]. Time-lapse monitoring aims to estimate the changes of medium properties as a function of time, where the difference in seismic data measured at different time steps should be caused by the change in the properties of the medium, which is in terms of seismic structure. "Time-lapse imaging" is also crucial in industry and applied to monitoring hydrocarbon and CO 2 storage, reservoir management and the long term monitoring of underground changes in a stress state and seismic features (e.g., [START_REF] Arts | Monitoring of CO2 injected at Sleipner using time-lapse seismic data[END_REF][START_REF] Picotti | Seismic modeling to monitor CO$_{2}$ geological storage: The Atzbach-Schwanenstadt gas field[END_REF][START_REF] Riahi | Time-lapse analysis of ambient surface wave anisotropy: A three-component array study above an underground gas storage[END_REF][START_REF] Planès | Time-lapse monitoring of internal erosion in earthen dams and levees using ambient seismic noise[END_REF][START_REF] Stork | Assessing the potential to use repeated ambient noise seismic tomography to detect CO 2 leaks: Application to the Aquistore storage site[END_REF][START_REF] Cotton | Time-lapse velocity analysis -Application to onshore continuous reservoir monitoring[END_REF][START_REF] Solazzi | Seismic Signatures of Fractured Porous Rocks: The Partially Saturated Case[END_REF][START_REF] Dupuy | Combined geophysical and rock physics workflow for quantitative CO2 monitoring[END_REF][START_REF] Yuan | High-frequency localized elastic fullwaveform inversion for time-lapse seismic surveys[END_REF]. [START_REF] Arts | Monitoring of CO2 injected at Sleipner using time-lapse seismic data[END_REF] compared the measured seismic data to estimate and predict the time-lapse migration of the CO 2 after injecting CO 2 into the underground reservoir. [START_REF] Riahi | Time-lapse analysis of ambient surface wave anisotropy: A three-component array study above an underground gas storage[END_REF] obtained the phase velocities and back azimuth angles of Rayleigh and Love waves processing by three-component beamforming algorithm and used them to perform a time-lapse analysis of an underground gas storage facilities. [START_REF] Planès | Time-lapse monitoring of internal erosion in earthen dams and levees using ambient seismic noise[END_REF] applied seismic interferometry to monitor the temporal changes in earthen dams and levees. The surface wave velocity reduction has been estimated in a laboratory-scale experiment and observed on a field embankment. [START_REF] Cotton | Time-lapse velocity analysis -Application to onshore continuous reservoir monitoring[END_REF] developed a 4D prestack time migration to detect velocity changes obtained after cross correlation of the base and monitor surveys. This method is applied to heavy oil high resolution data for revealing the time-lapse steam propagation within the reservoir. Time-lapse inversion is then developed by inverting the data measured at different times to obtain the changes of model parameters.

The subtraction of seismic data and synthetics gives the information of structural variations in the medium. The migration imaging technique is used in this study on experimental data to obtain the structural variation during deformation experiments.

In addition to understanding the physics of earthquakes, rock physical effects are widely recognised as a key role in the quantitative interpretation of seismic data in the field of CO 2 storage. Therefore, numerous studies related to petrophysics have looked for connections with seismic monitoring research and vice versa [START_REF] Vanorio | The rock physics basis for 4D seismic monitoring of CO2 fate: Are we there yet?[END_REF][START_REF] Picotti | Seismic modeling to monitor CO$_{2}$ geological storage: The Atzbach-Schwanenstadt gas field[END_REF][START_REF] Njiekak | CO2 rock physics as part of the Weyburn-Midale geological storage project[END_REF][START_REF] Prasad | CO2 messes with rock physics[END_REF]. Considering different rock models and analysing the sensitivity of the seismic response with varying saturations of CO 2 , [START_REF] Picotti | Seismic modeling to monitor CO$_{2}$ geological storage: The Atzbach-Schwanenstadt gas field[END_REF] built a petro-elastic model to detect the potential leakages in the reservoir. [START_REF] Njiekak | CO2 rock physics as part of the Weyburn-Midale geological storage project[END_REF] modelled and analysed the ultrasonic waveforms recorded during CO 2 saturated experiments on rocks.

The drop of P and S-wave velocities were found to mark the CO 2 phase. [START_REF] Prasad | CO2 messes with rock physics[END_REF] built a system for the acoustic properties of CO 2 -water interaction and applied time-lapse micrographs to monitor the effects of the seawater-sandstone interaction and the interaction between CO 2 -rich-water and carbonate. Combining the knowledge of rock physics, fluids and seismic waves properties, they studied the changes in the seismic waves propagating inside the rocks owing to CO 2 injection, which was discovered to be far from being a pure fluid-substitution problem.

The time-lapse relationship between seismic and rock physical properties helps reveal the earthquake nucleation process and is further widely used in industry. Therefore, it is necessary to develop numerical models and conduct laboratory experiments to determine the characteristics of rocks under stressed circumstances accurately.

Laboratory-scale measurements using seismic methods, or time-lapse seismic monitoring itself, are often challenging. The lack of studies conducting laboratory and numerical experiments on relatively small specimens under well-controlled conditions can cause discrepancies in interpreting the experimental data. Even though there are studies considering, for instance, the impact of the size of ultrasonic transducers (Kovalyshen et al., 2020), transducer effect, geometry effect [START_REF] Yoshimitsu | Geometric effect on a laboratory-scale wavefield inferred from a three-dimensional numerical simulation[END_REF], which are necessary for obtaining adequate information in waveforms, there have been few studies that quantitatively analysed time-lapse active seismic waveforms during laboratory experiments (e.g., [START_REF] Lai | Rock deformation monitoring using Monte Carlo waveform inversion[END_REF]. Therefore, in this study, we want to systematically combine laboratory experiments, numerical modelling and time-lapse seismic monitoring.

Rock is, by nature, a complex mixture of tiny grains that align either randomly or systematically. Seismic waves are sensitive to different spatial averages of the physical properties of the medium that they traverse through. In heterogeneous media such as rocks, scatters of seismic waves can lead to significant features in the coda waves. Wave interferometry is commonly used to monitor the minor changes in coda waves, which are multiply-scattered waves that are recorded after the arrival of the body waves [START_REF] Singh | Coda Wave Interferometry for Accurate Simultaneous Monitoring of Velocity and Acoustic Source Locations in Experimental Rock Physics[END_REF][START_REF] Obermann | 4-D Imaging of Subsurface Changes with Coda Waves: Numerical Studies of 3-D Combined Sensitivity Kernels and Applications to the $$M_{\mathrm{w}}$$7.9, 2008 Wenchuan Earthquake[END_REF]Zotz-Wilson et al., 2020). However, most of the studies for monitoring weak changes within the earth's subsurface using multiply scattered/coda waves have been limited to identifying weak changes rather than understanding the physics of these changes. Compared to the standard method of measuring the changes in velocity by picking the first-arrival times or coda wave analysis, we use full-waveform analysis to extract the physical information on the structural changes.

One of the sophisticated ways to simulate the scatters in rocks is by introducing random heterogeneity to numerical models using the Born series. Random heterogeneity generates incoherent scattered waves when seismic waves propagate through a heterogeneous medium. The scattered waves are overlapped with the source signal, and complex waveforms are observed at the stations. We may regard the scattered waves as noise in the observed waveforms. However, the scattered waves are not the same as the noise, randomly appearing in the time series. Elastic waves propagating through different laboratory-scaled rock samples with varying scales of heterogeneity have been studied (e.g., [START_REF] Sivaji | A physicalmodel study of the statistics of seismic waveform fluctuations in random heterogeneous media[END_REF]. By applying Born approximation, seismic scattering problems caused by small-scale and random heterogeneity can be solved [START_REF] Hudson | The use of the Born approximation in seismic scattering problems[END_REF][START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition[END_REF]. However, there is no research systematically considering and analysing at the same time the experimental laboratory setting, source wavelet input, waves propagation, waveform phases, random heterogeneity. This thesis wants to fill this gap.

This study first introduces our numerical modelling method throughout the whole re-search. We then present the Monte Carlo waveform inversion for estimating elastic and anelastic parameters and their evolution during laboratory rock deformation experiments. The transducer-to transducer one-source one-station active seismic data of dry and water-saturated samples are obtained from Zaima and Katayama (2018).

Based on the comparisons, we recovered the time-lapse variations of seismic velocities and attenuation of rock samples during deformation experiments, which we interpreted as crack developments. The systematic analysis by comparing laboratory and numerical experiments in metal and rock examples are reported later in this thesis. Numerical modelling was conducted to distinguish the phases in the waveforms and the effect of the instruments. We explored the phases of the experimental waveforms by computing the seismic wavefield modelling using the spectral element method (SEM) and calculating the sensitivity kernel. Random heterogeneity was introduced to the model to include scatters in the synthetic waveforms and approximate the fundamental physics of rock. This study shows that heterogeneity and grain sizes can be obtained by matching the experimental waveforms to the synthetic waveforms simulated from our models. Anisotropy of the rock samples can also be calculated effectively and through our method. Our study can be powerful since we can predict the physical properties of rocks by sending active-source seismic signals and matching the experimental waveforms with the synthetic waveforms generated by our models. The method is highly likely to be further used to predict the rocks' type effectively and further serve the structural imaging during the deformation and fracturing of the rocks. Chapter 7 introduces the migration imaging method using laboratory seismic data and their perspectives for exploring the time-lapse changes during rock deformation. Chapter 3 and chapters 4-6 are two independent scientific papers. The publications at national and international conferences are also available:

• Lai, S., Fuji, N., Katayama, I. ( 2019

Chapter2

Numerical modelling of wave propagation in a centimetric cylinder

Introduction

To calculate elastic wavefield in a medium with curve discontinuities or free surface boundaries, we can use a tensorial formulation of the wave equation and curvilinear coordinate transformations to make the computational uniform grid coincide with the curved interface in the physical domain (e.g., Komatitsch and Tromp, 1999). First, we used the optimally accurate finite difference method (OPT-FD) on the collocated Cartesian introducing zig-zag free-surface discontinuities (Lai et al., 2019). We then used ABAQUS software (Dassault Systeme), which numerically resolve the wave equation with the finite element method (FEM), to build the mesh of the model and perform the simulations since FEM implements the traction-free condition on the boundary and it is simple to adapt the mesh to a curved surface [START_REF] Igel | Computational seismology: a practical introduction[END_REF]. After comparison using 2D numerical modelling, we confirmed that the response function estimation could roughly rebuild the waveform, except for the attenuation influences by considering all the effects, such as 3D effect, instrumental response, attenuation.

However, the large discrepancy between observed and synthetic waveforms before convolution with the response function implied the necessity to consider a more real-istic configuration and design the laboratory experiments by changing the size of the sample and experimental geometry. By extending the study to the 3D dimensions, we can understand the weaknesses and ruptures in realistic conditions. Therefore, we implemented 3D numerical modelling using spectral element method (SEM) in most of the studies. Note that this study does not focus on the numerical modelling scheme itself: we thus review the linear seismic wave propagation straightforwardly in the following section; we then discuss the advantages and disadvantages of the methods tested.

Brief review of linearised seismic wave propagation

In order to simulate synthetic waveforms for an active source in a centimetric cylinder, we use linearised elastic equation and Born approximation for 3D random heterogeneities. The expressions below, i.e. equation of motion and reciprocity theorem, are common in seismology (e.g., [START_REF] Aki | Quantitative Seismology[END_REF] and general terms of partial derivatives (Fréchet derivatives) of waveforms can be found in many pieces of literature [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF][START_REF] Geller | Two efficient algorithms for iterative linearized inversion of seismic waveform data[END_REF][START_REF] Dahlen | Theoretical Global Seismology[END_REF]Takeuchi and Geller, 2000;[START_REF] Fichtner | Full seismic waveform modelling and inversion[END_REF][START_REF] Zhao | An efficient and flexible approach to the calculation of three-dimensional full-wave Fréchet kernels for seismic tomography-II. Numerical results[END_REF][START_REF] Fuji | Finite-frequency structural sensitivities of short-period compressional body waves[END_REF]. Here, we would like to briefly review the seismic wave propagation phenomenon and the reconstruction of partial derivatives from a technical point of view of dimensions and units, which is indispensable in our following random heterogeneity experiments. This practice aims to clarify units of all the quantities used in Fréchet derivative computation and avoid frequent dimension errors. For this reason, we describe all the concrete units in the international system of quantities (SI), which can be more helpful than a mere dimension description. We write a linearised elastic equation of motion in Cartesian coordinates:

ρ ∂ 2 ∂t 2 u i -Λ i j pq u p,q , j = f i for (t,x x x) in [0, T ) × Ω (2.1) with ρ [kg m -3 ] the density, Λ i j pq [N m -2 ] the elastic moduli, f i [N m -3 ] the i-th external force density, t [s] the time and u i [m] the i-th component of displace- ment. T [s]
is the duration of wave propagation and Ω denotes the whole volume of medium. The subscript , j denotes the partial derivative with respect to j-th coordinate and Einstein's summation convention is applied unless specified. Boundary conditions should be added to the left-hand side for [0, T ) × ∂ Ω. ρ is time invariant but Λ i jpq can be time variant when dealing with anelasticity as we discuss in the following sections.

Green's function, double-couple source, reciprocity

We first define Green's functions as follows:

ρ ∂ 2 ∂t 2 G il -Λ i j pq G pl,q , j = δ il δ (x x x-s s s)δ (t) for (t,x x x; 0,s s s) in [0, T )×Ω×Ω (2.2) with G il [m N -1 s -1 ]
is the i-th component of Green's function for a l-th direction impulse 1 [N s] exerted at the source point s s s at time t = 0. δ denotes either Kronecker's delta or Dirac's delta function. Note that the dimension of Dirac's delta function is the dimension of the inverse of its argument. For an arbitrary pair of x x x and s s s these Green's functions can be prepared. For a time-marching schemes, it is hardly possible to directly compute Green's functions, and one usually needs to use sinc or Ricker wavelets to generate pseudo-Green's functions. We thus describe the explicit form of displacement for a point force with an finite source-time function:

u n (t,x x x; 0,s s s) = t 0 dτS l (t -τ)G nl (τ,x x x; 0,s s s) for (t,x x x) in [0, T ) × Ω. (2.3)
where S l (t) [N] is a source-time function excited at s s s. This equation can be further simplified if we can assume a constant radiation pattern throughout an event S l (t) = S l S (t) with a demensionless normalised rate function S (t) as follows:

u n (t,r r r; 0,s s s) = t 0 dτS (t -τ)S l G nl (τ,r r r; 0,s s s) for (t,r r r) in [0, T ) × Ω. (2.4)
This equation can be used to apply the reciprocity theorem for three-component displacements and a double-couple source, as it will be shown later in this section.

When we discuss displacements due to a point double-couple source with a moment tensor rate function M i j (t) [N m] at s s s with a Dirac source time function, we can convolve the source time function with synthetics and integrate over the source volume when it is necessary. We replace f i in equation 2.1 with double-couple source:

f i (t,x x x) = -δ (x x x -s s s)M i j, j (t).
(2.5)

The representation theorem for such a source can be given by the volume integral over the whole rock and the aid of Gauss' theorem as follows:

u n (t,r r r; 0,s s s) = t 0 dτM i j (t -τ)G ni, j (τ,r r r; 0,s s s) for (t,r r r) in [0, T ) × Ω (2.6)
with the local Cartesian partial derivative operator , j which can be applied regardless of the choice of global coordinates in the system (in a global Cartesian system, we can use the same operator on any point). This equation can be further simplified if we can assume a constant radiation pattern throughout an event:

M i j (t) = M i j S (t)
with a demensionless normalised rate function S (t) as follows:

u n (t,r r r; 0,s s s) = t 0 dτS (t -τ)M i j G ni, j (τ,r r r; 0,s s s) for (t,r r r) in [0, T ) × Ω. (2.7)
This equation is simpler to handle together with equation 2.4 since we can use the same dimensionless "source-time function" S (t) for the reciprocity theorem below.

When we propose the equations of motion 2.2 for a set of two Green's function G in (t,x x x; 0,x x x A ) and G mi (t,x x x; 0,x x x B ), convolving them with G mi (t,x x x; 0,x x x B ) and G in (t,x x x; 0,x x x A ) respectively, taking a volulme integral over the whole rock, we obtain the reciprocity theorem with the aid of Gauss' theorem as:

G mn (t,x x x B ; 0,x x x A ) = G nm (t,x x x A ; 0,x x x B ).
(2.8)

Plugging this result into equation 2.7, we obtain:

u n (t,r r r; 0,s s s) = M i j t 0 dτS (t -τ)G in, j (τ,s s s; 0,r r r) for (t,r r r) in [0, T ) × Ω. (2.9)
In our application, we need to compute the component perpendicular to the rock's surface due to an explosion source (diagonals of M i j ) and a shear source (nondiagonals). We thus use this equation by simulating the seismic wavefield with a n-directed single force and S (t) located at the receiver r r r, computing the strain at the source point s s s so that we can either obtain explosion-related and shear-related seismograms.

Born approximation

To simulate synthetic waveforms for random 3D heterogeneity added to a homogeneous model, we propose to use Born approximation (the first-order approximation of Born series). We re-define homogeneous Green's functions as follows:

ρ 0 ∂ 2 ∂t 2 G il -Λ 0 i j pq G pl,q , j = δ il δ (x x x -s s s)δ (t) for (t,x x x; 0,s s s) in [0, T ) × Ω × Ω.
(2.10)

We now decompose the wavefield u i into the initial (homogeneous) model-related wavefield and the 3D perturbation to the initial model:

u i = u 0 i + δ u i + δ 2 u i + δ 3 u i + • • • for (t,x x x) in [0, T ) × Ω (2.11)
where the δ u i is the first-order scattering and δ 2 u i the second-order scattering, etc.

with elastic parameter perturbation:

Λ i j pq = Λ 0 i j pq + δ Λ i j pq .
(2.12)

In this thesis, we are interested in P-wave velocity perturbation. Thus, we omit the discussion of the density perturbation δ ρ, which requires displacement Green's functions instead of strain Green's functions. When substituting this decomposition into equation 2.1 for a source at s s s, we obtain Born series:

ρ 0 ∂ 2 ∂t 2 u 0 i -Λ 0 i j pq u 0 p,q , j = S l (t)δ il δ (x x x -s s s) for (t,x x x) in [0, T ) × Ω (2.13) ρ 0 ∂ 2 ∂t 2 δ u i -Λ 0 i j pq δ u p,q , j = -δ Λ i j pq u 0 p,q , j for (t,x x x) in [0, T ) × Ω (2.14) ρ 0 ∂ 2 ∂t 2 δ 2 u i -Λ 0 i j pq δ 2 u p,q , j = -δ Λ i j pq δ u p,q , j for (t,x x x) in [0, T ) × Ω (2.15) • • • Equation 2.
13 is the equation of motion for the initial model with the physical sources. The right-hand sides of all the equations are the same operator for the initial model, and the left-hand sides can be interpreted as physical/virtual sources.

The main interest of the Born series is that we can recycle the (simpler) operator for the initial models (e.g., Takeuchi and Geller, 2000). While the complete set of Born series including high-order scatterings provides an exact solution, we truncate this series at equation 2.14 to handle the first-order scattering only (Born approxi-

mation) when |Λ 0 Λ 0 Λ 0 | >> |δ Λ δ Λ δ Λ| and thus |u 0 u 0 u 0 | >> |δ u δ u δ u| >> |δ 2 u δ 2 u δ 2 u|.
The solution for the zeroth-order equation 2.13 is given by the convolution form:

u 0 i (t,x x x; 0,s s s) = t 0 dτS l (t -τ)G il (τ,x x x; 0,s s s) for (t,x x x) in [0, T ) × Ω. (2.16)
Note that time integral is limited from 0 to t if there is no energy before 0 in S l (t).

The solution for Born approximation equation 2.14 is thus given by: δ u i (t,x x x; 0,s s s)

(2.17)

= - Ω d 3 η η η t 0 dτ δ Λ l j pq (η η η)u 0 p,q (t -τ,η η η; 0,s s s) , j G il (τ,x x x; 0,η η η) (2.18) = - Ω d 3 η η η t 0 dτ τ 0 dτ δ Λ l j pq (η η η)S k (t -τ -τ )G pk,q (τ ,η η η; 0,s s s) , j G il (τ,x x x; 0,η η η) (2.19) = Ω d 3 η η η t 0 dτ τ 0 dτ δ Λ l j pq (η η η)S k (τ )G pk,q (t -τ -τ ,η η η; 0,s s s)G il, j (τ,x x x; 0,η η η) (2.20) = Ω d 3 η η η t 0 dτ τ 0 dτ δ Λ l j pq (η η η)S k (τ )G pk,q (t -τ -τ ,η η η; 0,s s s)G li, j (τ,η η η; 0,x x x) for (t,x x x) in [0, T ) × Ω.
(2.21) Equation 2.17 can be obtained from equation 2.16 plugging the virtual source term in equation 2.14. We then substitute u 0 by equation 2.16 in order to obtain equation 2.18. When taking integration by parts for a larger Ω than necessary and applying Gauss' theorem to make the surface integral vanish:

∂ Ω d 2 η η η t 0 dτ τ 0 dτ δ Λ l j pq (η η η)S k (τ )G pk,q (t -τ -τ ,η η η; 0,s s s)G il (τ,x x x; 0,η η η)n j = 0
with n j the j-th component of normal vector, we can obtain equation 2.19. This equation is in convolutional form which can be iteratively used for higher Born series. However, for first-order Born approximation, one can make use of reciprocity (that we do not prove here) to further obtain equation 2.20 which requires the computation of only two Green's function from the source s s s and from the receiver x x x in our example: as is well accepted in seismology, the Born approximation can be described as the cross-correlation of the two Green's functions.

In this study we are interested in P-wave velocity perturbation. We take an isotropic approximation:

Λ l j pq = λ δ li δ pq + µδ l p δ jq + µδ lq δ j p (2.22)
with Lamé's parameters λ and µ [N m -2 ] that relates the P-wave velocity α [m s -1 ]:

α = λ + 2µ ρ . (2.23)
Taking the total derivative of the above equation, we obtain:

dα = λ + 2µ 2α ln ρdρ + 1 2ρα dλ + 1 ρα dµ.
(2.24)

When setting λ to be the only non-zero perturbed parameter, we can simplify the elastic moduli perturbation as:

δ Λ l j pq = 2ρ 0 α 0 δ αδ l j δ pq .

(2.25)

We thus can obtain the Born approximation for δ α perturbation as:

δ u i (t,x x x; 0,s s s) = 2 Ω d 3 η η η t 0 dτ τ 0 dτ ρ 0 α 0 δ αS k (τ )G pk,p (t -τ -τ ,η η η; 0,s s s)G li,l (τ,η η η; 0,x x x).
(2.26)

We thus need to store the sum of diagonals of strain Green's functions (e.g., [START_REF] Fuji | Finite-frequency structural sensitivities of short-period compressional body waves[END_REF]. This equation can be directly used for synthetic seismogram computation for a 3D heterogeneous model, by truncating at the second term of the left-hand side of equation 2.11 as we will show in some examples later on.

Another application of this equation is to illustrate the waveform Fréchet derivative for the P-wave velocity at each point. If we apply a delta-function perturbation:

δ α(η η η) = 1 [m s -1 ] × δ (η η η -q q q) [m -3 ],
we obtain the waveform P-wave velocity Fréchet derivative:

∂ u i (t,x x x; 0,s s s) ∂ α(q q q) [s m -3 ] = 2ρ 0 (q q q)α 0 (q q q) t 0 dτ τ 0 dτ S k (τ )G pk,p (t -τ -τ ,q q q; 0,s s s)G li,l (τ,q q q; 0,x x x).

(2.27)

Mapping Fréchet derivatives helps us to understand the sensitivity of seismic parameters to the waveforms (or their attributes such as traveltimes and amplitudes).

Practically speaking, it could be helpful to stock Fréchet derivatives using equation 2.27 to realise Born approximation with the discretised version of equation 2.26 as follows:

δ u i (t,x x x; 0,s s s) = N p ∑ I=1
∂ u i (t,x x x; 0,s s s) ∂ α(q q q I ) ∆α(q q q I )∆V I (2.28) with ∆α(q q q I ) [m s -1 ] P-wave velocity perturbation averaged for a small portion of volume ∆V I [m 3 ] in the vicinity of I-th point q q q I and N p the number of perturbed voxels. We use these Fréchet derivatives to illustrate temporal sensitivity and perform migration in the later chapters. Note that all the discussions above can be extended to any elastic and anelastic parameters straightforwardly (e.g., [START_REF] Fuji | A methodology for inversion of broadband seismic waveforms for elastic and anelastic structure and its application to the mantle transition zone beneath the Northwestern Pacific[END_REF][START_REF] Zhao | An efficient and flexible approach to the calculation of three-dimensional full-wave Fréchet kernels for seismic tomography-II. Numerical results[END_REF].

Since the migration image is equivalent to the zeroth iteration of full-waveform inversion (or the gradient direction) (e.g., [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF][START_REF] Fuji | Small-data seismology-from a centimetric rock sample to Mars[END_REF] which can be expressed as the correlation between the waveform residual and Fréchet derivatives (e.g., [START_REF] Fuji | A methodology for inversion of broadband seismic waveforms for elastic and anelastic structure and its application to the mantle transition zone beneath the Northwestern Pacific[END_REF], we express the P-wave migration image M at q q q for data d i as follows:

M (q q q) ∝ ∑ i ∑ r r r ∑ s s s t 0 dt [d i (t,r r r; 0,s s s) -u i (t,r r r; 0,s s s)] ∂ u i (t,x x x; 0,s s s) ∂ α(q q q) (2.29) = 2ρ 0 (q q q)α 0 (q q q) t 0 dt ∑ s s s t 0 dτdτ ρ 0 α 0 δ αS k (τ )G pk,p (t -τ -τ ,η η η; 0,s s s) ∑ r r r G li,l (τ,η η η; 0,r r r) [d i (t,r r r; 0,s s s) -u i (t,r r r; 0,s s s)] .
(2.30)

Depending on the configuration, the migration image can be obtained by summing the cross-correlation between the forward wavefield and back propagated wavefield with the wavefield residuals from all the receivers r r r as in industry (equation 2.30).

In chapter 7, we illustrate the migration images for cylindrical samples but up to this stage we use equation 2.29, since we have a few receivers and sources but a large number of 3D points.

2D numerical modelling: OPT2D and ABAQUS

We used and modified the full-waveform-based elastic imaging methods (OPT2D)

for laboratory rock experiments [START_REF] Geller | Optimally accurate second-order time-domain finite difference scheme for the elastic equation of motion: one-dimensional cases[END_REF]Fuji et al., 2016;Lai et al., 2019). We used an optimally accurate finite-difference method on collocated

Cartesian coordinate with zig-zag free-surface discontinuities for the 2D circle forward modelling. We modified the boundary conditions for an optimally accurate collocated finite-difference method based on predictor-corrector scheme [START_REF] Takeuchi | Optimally accurate second order time-domain finite difference scheme for computing synthetic seismograms in 2-D and 3-D media[END_REF]Fuji et al., 2016) for 2D P-SV medium. The code can be found on Github https://github.com/seismobassoon/OptimallyAccurate2D with visualisation facilities. We performed the modelling in the same geometry and conditions as the experiments of Zaima and Katayama (2018) We take the velocity model from Zaima and Katayama (2018), comparing it with the observed waveform. It shows the shortcomings of using the hand-picking method with only first arrival to infer the velocity models. We use full-waveform for data fitting instead. However, the waveform misfit is too large to go beyond the first arrival fitting for the moment, and we thus further explored the way to better simulate the waveforms.

Many reasons cause the differences between synthetics and observed data. In our case of 2D circle disk forward modelling, 3D effect, instrumental response, attenuation, finite-surface source, and other unknowns may cause a considerable discrep-ancy between observed and synthetic data. Hence, we design the response function estimation, which involves all the impact factors, to rebuild the synthetics similar to the observed waveform. Since we have the relation:

F -1 F(u obs m1 )(t) F(u syn m1 )(t) • F(u syn m1 )(t) = u obs m1 (t), (2.31)
we can have the response function of the experiment by taking the relationship between the observed waveform and the synthetic waveform of the same model (m1).

Even though the time-lapse P-wave velocity can be slightly changed due to the first setting of P-wave velocity (m1), it will not be as slow as the V P based on hand-picked traveltimes, and we value more about the ∆V P but not the absolute V P in this study.

Then, we can obtain the response function R(ω) as follows:

R(ω) = F(u obs m1 )(t) • 1 F(u syn m1 )(t) (2.32)
We can rebuild the waveform similar to the experimental results by using only the numerical synthetics and convolving them with the response function:

R(t) * u syn m2 (t) ∼ = u obs m2 (t) (2.33)
We test several norms to automatically and accurately detect the time-lapse structure changes (figure 2.2). Therefore, we first comply with the eyeball-based trial-anderror method in order to find the best of the velocity models. We then introduce automated Monte-Carlo methods based on waveform misfit to construct a system to find the models fast and accurately. We processed a series of velocities with an interval of ∆v = 10 m/s, so that we have 191 velocity models for autofitting. The

Spectral element method and viscoelasticity

To understand the complicated phases of the experimental records, we simulate the waveforms in models with the approximate geometry used in experiments. The spectral element method (SEM3D software) simulates the waveform propagation to handle complicated geometry. We briefly review the spectral element method in this section, and the readers are referred to Komatitsch and Vilotte (1998); Chaljub et al. (2007); Peter et al. (2011) for more details.

We modified the strong form of the equation of motion 2.1 by multiplying it with an arbitrary test function w

(m)
i , taking a volume integral with the aid of integral by parts and Gauss' theorem to construct the weak form:

Ω dV w (m) i ρ ∂ 2 ∂t 2 u i + Ω w (m) i, j Λ i j pq u p,q = Ω w (m) i f i S f (t) + Ω w (m) i, j M i j S M (t) (2.34)
with S f and S M the normalised source-time function for a single force and a double-couple source, respectively. Note that we assume natural boundary condition (free surface) on ∂ Ω but any particular boundary condition can be added as an additional boundary condition. We expand the displacement u on the same set of basis functions as w i to construct Galerkin method as in Geller and Ohminato (1994); Komatitsch and Vilotte (1998) as:

u i (t,x x x) = ∑ n c n (t)w (n) i (x x x) (2.35)
with c n the coefficient for each basis function. The above equation can be reformulated in a matrix form as:

M M M∂ 2 t c c c(t) +K K Kc c c(t) = g g g(t), (2.36)
where M M M is the mass matrix:

M mn = Ω dV w (m) i ρw (n) i , (2.37)
and K K K is the stiffness matrix:

M mn = Ω dV w (m) i, j Λ i j pq w (n) p,q , (2.38)
and g g g is the excitation force vector:

g m (t) = Ω dV w (m) i f i S f (t) + Ω w (m) i, j M i j S M (t).
(2.39)

SEM introduces the combination of Lagrange polynomials as interpolants defined on

Gauss-Lobatto-Legendre (GLL) points for the elastic wave equation. SEM allows a very significant reduction in the complexity and cost of the algorithm since the mass matrix M M M is diagonal by construction.

Attenuation

Here, attenuation is described as viscoelasticity. We use three standard Zener linear solids (SLS: figure 2.8) (Liu et al., 1976) to approximate a nearly constant quality factor Q (Emmerich and Korn, 1987). SLS assumes that Q does not depend on the frequency and can be expressed as follows:

Q -1 (ω) = ∑ N r 1 ωθ l γ l θ 2 l +ω 2 1 + ∑ N r 1 ω 2 γ l θ 2 l +ω 2 , (2.40) θ l = ω 0 2π ω 1 ω 0 l-1 N r -1 , l = 1, ..., N r . (2.41)
We define a band of angular frequency ranging from ω 0 to ω 1 with N r points in logarithmic scale in order to minimize the Q -1 (ω) to the given Q -1 re f (ω). The values of γ in equation 2.40 is a combination of the stress and strain relaxation times associated to each SLS:

γ = 1 N τ ε τ σ -1 , (2.42)
where N is the number of SLS, γ is the modulus defect normalised by the relaxed modulus for each SLS and is determined following the assumption of a constant Q re f on the seismic band of interest. From a constant Q re f and a prior fixed τ σ chosen by [START_REF] Lombard | Numerical modeling of transient twodimensional viscoelastic waves[END_REF], τ ε can be determined by the classic least-square fit. As input parameters, we set the number of standard linear solids and the seismic periods of interest (the frequency range). The frequency is sampling for 1/(stress relaxation times) (1/τ σ ), which is the frequency sampling for Q chosen by [START_REF] Lombard | Numerical modeling of transient twodimensional viscoelastic waves[END_REF]. In this study, we use 3 SLS on a band of relaxation frequency ranging from 33.3 kHz to 3.3 MHz that the configuration of Q -1 to the frequency with

Q = 15 is shown in figure 2.8.
From Q re f (data) constant and τ σ (a priori fixed), gives τ ε for each standard linear solid for the time being: only a classical least-squares fit, that means linear operators in the inversion process. Now we can obtain values of τ σ and τ ε : the first is the inverse of ω τ s and the second can be given by:

κ = 1 N τ ε τ σ -1 (2.43) τ σ = 1 ω τ σ (2.44) τ ε = τ σ * (1 + N solid * κ).
(2.45) However, everything needed in the SEM attenuation process is caught with γ and

ω τ σ .
SEM3D, instead of Q P and Q S , requires the number of standard linear solid, N r , for attenuation, attenuation period band and model period. The greater N r , the smaller the error introduced by describing the medium with constant quality factor in terms of a finite number of relaxation mechanisms. At the same time, increasing N r greatly increases the computational cost, especially the memory requirements. In the practice, N r > 3 is rarely implemented in the literature, and N r = 1 is widely used, especially in the 3D context [START_REF] Robertsson | Viscoelastic finitedifference modeling[END_REF]. In this study, we use three solids for attenuation. The attenuation band is set on the period from 3.0×10 -7 s to 3.0×10 -5 s (Figure . 2.8).

Chapter3

0D monitoring of weakening process from granite deformation

Introduction

In this chapter, we want to develop a methodology that inverts the seismic waveforms for structural changes during laboratory experiments performed by Zaima and Katayama (2018). They performed a triaxial compression experiment on the Aji granite and carried out active seismic measurements during the rock deformation.

The idea is to robustly and quantitatively infer the subtle structural changes of elasticity and anelasticity, and three-dimensional structure in the near future. A small number of previous studies attempted "waveform inversion" of laboratory experiments. We use the spectral element method (SEM: see a brief review in chapter 2)

to model wave propagation in a rock sample during deformation. The spectral element method can accurately handle the boundary conditions in 3D complex models (Capdeville et al., 2003;Peter et al., 2011;De Basabe and Sen, 2014;Chaljub et al., 2015). However, due to its computational costs, there are only a few studies devoted to the comparison of SEM synthetics with laboratory experimental data (e.g., Pageot et al., 2017;Solymosi et al., 2018). Here, we perform SEM modelling on one-source and one-receiver experimental data to show the adaptability of SEM to the seismic data in the laboratory. Throughout 3D simulation, we systematically construct a database of waveforms for different sets of elastic and anelastic parameters: P-wave velocity (V P ), S-wave velocity (V S ), P-wave quality factor (Q P ) and S-wave quality factor (Q S ). We then perform a Monte Carlo waveform inversion using an objective function that combines three kinds of cost functions: L 1 -norm and L 2 -norm of waveform misfit and zero-lag cross-correlation. We would like to explore rock elastic properties and attenuation changes using numerical simulation and Monte Carlo inversion. Throughout our approach, we will show how numerical models help understand the observations and illustrate the seismic characteristics of rock deformation.

Main results and discussions can be found in a published article:

Methodology

This study conduct inversions to image elastic and anelastic structural changes during rock deformation. Hence, we first acquire the active seismic data during the deformation experiment. Meanwhile, we prepare a set of synthetic seismograms with the aid of the spectral element method. We then select the preferred models for the observed seismograms (compressional and shear strains) for each differential stress step to obtain the time-lapse trend of elastic and anelastic parameters. First, we revisit the experiment configurations briefly and provide known and unknown properties in the experiments. Second, we describe the numerical modelling scheme.

Finally, we present our data processing strategy based on the data and modelling conditions and detail the inverse problem implemented in this study.

Inversion problem

We use waveforms for the differential stress ∆σ :

d T (∆σ ) = d T PP (∆σ ) d T SS (∆σ ) (3.1)
where d PP and d SS are time series vectors of the compressional and shear strain waveforms recorded at the antipodal receiver (figure 3.1 a). The superscript T denotes transpose. We then generate synthetic seismograms u(m) with the same geometrical and source configurations for a set of seismological structural parameters m. We define an objective function to maximise:

S(m, ∆σ ) = α d T (∆σ )u(m) |d(∆σ )||u(m)| + β 1 |d(∆σ ) -u(m)| + γ 1 |d(∆σ ) -u(m)| 2 (3.2)
with α, β , γ, the weighting factors for zero-lag cross correlation, reciprocals of L 1norm and L 2 -norm, respectively. First, equation 3.2 is used to constrain the main phase with a large α. Second, equation 3.2 with large β and γ is used to investigate the amplitude of the waveform, and to further define the change in Q. The choice of (α, β , γ) depends on the dataset since β and γ are difficult to fix from one experiment to another. Here, we use our ad-hoc values as discussed in the results section, although this problem should and will be afforded in forthcoming studies.

As the Aji granite is fine-grained and nearly isotropic, we can assume a homogeneous medium with a wavelength of 2-3 mm. We use four parameters to represent the model vector m T = (V P ,V S , Q P , Q S ). This parameterisation is to homogenise 3D an-/elastic heterogeneities. The time-lapse change of these global parameters will be related to the microscopic short-wavelength structural changes such as crack generation. The anelastic attenuation can be an indicator of the 3D elastic heterogeneity and intrinsic attenuation. Based on the a priori information on the Aji granite, we set V P varying from 3700 m/s to 5800 m/s and Q P from 60 to 200 (Q -1 P from 5.0×10 -3 to 1.7×10 -2 ) for the P-wave. As for the S-wave, V S varies from 2660 m/s to 3550 m/s and Q S from 20 to 90 (Q -1 S from 1.1×10 -2 to 5.0×10 -2 ). Note that all the compar-ison is based on the waveform analysis for ∆σ = 0. We are interested in the relative evolution of the velocity and attenuation parameters from the initial status ∆σ , which can be expressed as δ m, instead of the absolute values of m. δ m can thus give us insight into the relative changes of the rock properties. To simulate realistic wave propagation, for the sake of seeking boundary conditions including silicone rubber over the rock, one needs to carefully mesh the medium in SEM computation. The standard GLL quadrature of SEM requires hexahedral meshes. However, creating such nonstructured hexahedral meshes with acceptable computational cost and high accuracy for complex geometry is challenging. This study used Trellis software (www.csimsoft.com) to mesh the model. In the case of nonstructured meshes, one has to avoid distorted elements, which may cause numerical instabilities.

Traxial compression experiments

Numerical models

Figure 3.3 shows the meshing used in our forward modelling. To avoid the numerical error caused by the contact between more than one medium, we must include a layer as refinement at the boundaries of different materials. Since we use polynomials of degree five as trial functions inside each element, the wavelength should be larger than the element size. Considering the minimum S-wave velocity of 2500 m/s in the granite rock and the maximum target frequency of 2 MHz, the mesh size is set up to be approximately 0.3 mm in all models.

Boundary conditions

The rock samples used in laboratory experiments are roughly cylindrical rocks. However, the sample was cut to form two sharp parallel surfaces, forming an incomplete cylinder. This imperfection of the cylindrical nature can cause a non-negligible effect on the waveforms. Although the rock sample itself is the main target of our simulation, a single-material approximation with free surface boundary would fail in waveform modelling due to the surrounding materials. Materials such as silicone rubber, glue, aluminium frame, shrink plastic attached to the granitic rock may affect the waves propagating in the rock. Therefore, we incorporate complex boundary conditions to find a suitable model setting without losing computational efficiency. Furthermore, we investigate the boundary conditions with an additional physical layer of Shin-Etsu KE45W RTV silicone rubber (Table 3.1) to simulate the transducer waveforms and obtain the preferred model (figure 3.3). The additional boundary As shown in the figure 3.4 a and b, the free-surface boundary rock model cannot reproduce the observations due to its complete reflections from boundaries. Hence, we attempt to modify the boundary conditions by considering the silicone rubber, as in the real experiments: a layer of glue and silicone rubber is attached to the rock sample on its cylindrical side to jacket the sample. This layer isolates the oil from the rock sample and transducers. The physical parameters of this layer are calculated from the characteristics of the silicone rubber used in the experiment (Table 3.1). The layer is added to the cylindrical and the side-cut cylindrical rock model. the more features added, the more similar the synthetics to the observed waveform.

We thus concluded to use model (d) as our boundary conditions.

The additional layer in the model (d) with the characteristics of silicone rubber can provide a buffer-like region that slows down the waves and reduce the amplitude of quasi-surface waves propagating between the two different materials. Figure 3.5

shows the coefficients of P-to-P and S-to-S transmission (Lay and Wallace, 1995), calculated as amplitude and phase at the solid-solid interface from rock to the silicone rubber. When the incidence angle is 0, the P-to-P transmission coefficient is as high as 1.75, which is coherent to our observation that the direct P waves are amplified when they transmit and reach the interface between the two materials, and travel through a low-Q and low-velocity material. The input wavelet of the rock experiment (top) and the numerical modelling (bottom). In order to avoid the noise of the experimental wavelet, we use the filtered wavelet as the input.

Source and transfer functions

In addition to the boundary setting and model geometry, the source time function and the source-receiver relative position will also cause a significant impact on the waveforms. We set the P-wave source as an isotropic explosion source and the S-wave source to be shear stress along the xz direction to mimic the laboratory experiment. The source time function is formally taken from the laboratory experiment input waveform. However, the absolute amplitude of the input signal transmitted to the rock by the piezoelectric transducers is unknown. Therefore, in this study, we estimated the effective source time function by a trial-and-error approach. We concluded that the recorded input waveform low-pass filtered at 2 MHz is the best source time function to be used. This is the same dominant frequency band used in a laboratory experiment (figure 3.6).

Compressional and shear waves are transmitted independently at different source positions during the experiments, and their strain waveforms are recorded at their corresponding antipodal points. In reality, the transducers are glued to the rock sample. However, strains are not always continuous, and we have less knowledge on the coupling between the transducers, silicone rubber, and the rock sample. We thus consider four possible effective source-receiver distributions (figure 3 , the first phases up to the fifth peak (4.0-6.2 µs) have higher amplitudes than the other configurations. Therefore, we use the setting (b) as the source-receiver configuration in this study. We will discuss in detail the possible reasons for this preference in section 3.5.

Choice of a range of seismic parameters for Monte Carlo inversions

After determining the preferred model for the initial state m 0 (figure 3.3), we performed 6000 simulations of different velocity and attenuation models. Since the volume of the rock changed little in the experiment, we did not consider the variation in the rock density caused by the deformation (the density is fixed to be 2650 kg/m 3 ).

We used constant elastic and anelastic parameters of the silicone rubber jacket, whereas we systematically modified the parameters of the rock model to create a large database as shown in table 3.2. We first simulate the P-wave velocity from 4700 m/s to 5800 m/s with a step of 20 m/s, and from 3700 m/s to 4600 m/s with a step of 100 m/s. We ranged S-wave velocity from 2660 m/s to 3550 m/s with a step of 10 m/s (cf. Watanabe and Higuchi, 2015;Zaima and Katayama, 2018).

Based on the above considerations, in our first step of data processing, for P-wave waveform simulation, we vary the elastic parameters of V P , Q P , and the V P /V S while fixing Q S to be constant. During the P-wave analysis, we vary V P /V S from 1.5 to 1.7 with step equals to 0.05 (cf. Zaima and Katayama, 2018). For the S-wave waveform simulation, we change the elastic parameters of V S , Q S , and the V P /V S while keeping Q P constant. During data processing, we derive the best models using equation 3.2, refining the range of parameters to reach the best values. We use our preferred models during the first numerical simulations to localise our inversions.

Processing of observed and synthetic data

The study aims to directly compare the observed and synthetic waveforms to infer the structural changes during the deformation experiment. By carefully treating 3D geometry, attenuation and boundary conditions as we discussed above, our method is advantageous to foster understanding of rock deformation compared to our previous studies (e.g., Lai et al., 2019). We had to introduce transfer functions to compensate unrealistic numerical settings such as 2D assumptions, pure-elastic approximations and free-surface boundary conditions. According to the FFT spectrum (figure 3.8 and 3.9) and time-dependent spectrogram analysis in figure 3.10, the strongest energy of the experimental waveform is concentrated in the frequency band between 1.6 MHz and 3.0 MHz for P-wave records and in the frequency band between 1.0 MHz and 2.8 MHz for S-wave. Compared to S-wave temporal spectrogram that indicates a constant dominant frequency range, the P-wave spectrogram shows lower frequency energy (∼500 kHz) after 7 µs in addition to the most concentrated high-frequency region at the beginning. In this study, however, we focus on the variation in the high-frequency regime to capture the influence caused by inner structure changes. Two Butterworth bandpass filters with a frequency range from 1.5 MHz to 3.5 MHz and from 1.0 MHz to 3.0 MHz are applied on the P-and S-wave data, respectively. After filtering, a Hann window is applied to data to extract the first arrivals. Note that we used the same data processing treatments both on observed and synthetic waveforms. 3.16 shows snapshots of the P-wave strain wavefield in our preferred model for ∆σ = 0 from 2 µs to 7 µs with a time step of 1 µs. These snapshots are crosssections of the wavefield at the same height as the source and receiver positions. At t = 2 µs, the snapshot shows the direct P-wave radiation from the source towards the antipode. No scattering or conversion at the sample surface occurs at this time.

However, the wavefronts are no longer a regular concentric form due to the sidecut surface. The surrounding low-velocity layer outside the sample partially absorbs energy and partially reflects it. This is why the reflection from the sample surface becomes too weak to interfere with the direct P waves. Since our model is not a perfect cylindrical model, the two cut corners of the source side immediately provoke reflections. When the direct wave is reflected at the opposite side (t = 4 µs), the reflected P phase (PP wave) is generated and propagates back to the source side.

Since the receiver side has two cut corners as well, two additional reflected waves are generated and interfere with each other between 5 µs and 6 µs (figure 3.16d and e). However, the external low-velocity silicone layer traps some energy from the direct waves, and the reflections get smaller and smaller. The P-wave propagation snapshot at a later time (t = 7 µs) shows a complex feature due to the contamination of reflected and converted waves developed at the side-cut surface.

Figure 3.17 presents the snapshots of the S-wave wavefield. Most side-cut geometry effects mentioned above can also be seen in S-wave snapshots. The S wave arrives at the receiver at about 6 µs. The strong influence of guided waves along the sample surface differs in the S-wave propagation from P waves. The direct S wave (body wave) at the antipode receiver is contaminated by the waves coming subsequently along with the interface (quasi-surface waves) (figure 3.17e). This indicates the difficulty of extracting information on the sample's interior using S waves alone. Among these three cost functions, zero-lag cross-correlation can localise best models with less deviation since kinematic traveltimes are important in the first place (cf. Fuji et al., 2021). L 1 -norm and L 2 -norm both yield a poorer waveform matching at higher compression stress conditions than at lower stress conditions, which indicates the need for further parameter search, such as 3D heterogeneity and/or attenuation. P-and S-wave waveform fitting results on the dry sample are shown in figure 3.19.

As we discussed in subsection 3.4.1, the S-wave waveform is more difficult to fit than the P-wave waveform. Figure 3.20a shows the comparison of the best-fit V P and V S models of the dry sample as well as the hand-picking results from Zaima and Katayama (2018), indicating the evolution of the P-and S-wave velocity. For the experiment under dry conditions, the evolution of V P and V S in general are similar and follows three stages: (1) increasing stage, where V (= V P ,V S ) increases as the ∆σ increases ; (2) transition stage, ∆V remains small; and (3) decreasing stage, V drops rapidly. The highest P-wave velocity structure of the dry experiment occurs when the difference between the compression stress and the confining pressure is equal to 81-128 MPa. When P-wave velocity reaches the highest, the P wave remains stable.

After ∆σ = 250 MPa, V P drops abruptly, and the sample fails. The peak of V S value occurs at ∆σ = 105 MPa slightly after that of V P (81 MPa). V S slowly increases until the differential stress reaches its peak value and decreases until ∆σ =250 MPa, where the V S starts dropping suddenly. In the dry experiment, the difference between the fastest and slowest V P is 940 m/s (∼18%). As for V S , the difference is about 560 m/s (∼16%). The increase in velocity is larger for compressional waves (∼100 m/s) than for shear waves (∼40 m/s) (figure 3.20a).

Similarly, the velocity evolution of the wet sample (figure 3.20b) follows the three stages of the dry sample. In the wet experiment, the difference between the maximum and the minimum V P is 320 m/s (∼6.1%) and the V S difference is about 425 m/s (∼13.7%). Figure 3.21 shows that increasing compression stress has a more significant impact on V P changes in dry data than in wet data. The dV S of wet data decreases most when having the same stress stage among all the data sets. Furthermore, the turning point of velocity changes of wet data shows earlier (∼60 MPa) than dry data (∼100-150 MPa). The increasing and transition stages of the wet data are much shorter than those of dry data.

Attenuation parameters evolution

The variations of Q -1 P and Q -1 S related to the increasing compression stress are compared in figure 3.22a. For both dry and wet experiments, attenuation shows a different behaviour than velocity during effective stress increases. For dry experiment, Q P varies initially from 250 to 10 (Q -1 P varies from 4 × 10 -3 to 1 × 10 -1 ) and Q S varies from 100 to 20 (Q -1 S varies from 1 × 10 -2 to 5 × 10 -2 ) during the deformation. Unlike velocity change, Q -1 P has a minor variation when ∆σ varies from 0 to 298.2 MPa. Only after ∆σ reaches 254 MPa, Q -1 P starts to rise. When ∆σ exceeds 324 MPa, Q -1 P grows rapidly. Q -1 S has a similar behaviour as Q -1 P during rock deformation. Yet, Q -1 S rises greatly when ∆σ exceeds 298.2 MPa, which is later than the Q -1 P . It shows that the internal structural changes of the rock have a greater impact on the attenuation of P waves in dry cases. For the wet experiment, Q P varies initially from 130 to 170, and from 170 to 20 (Q -1 P varies from 7.7 × 10 -3 to 5.9 × 10 -3 , and from 5.9 × 10 -3 to 5 × 10 -2 ). Q S varies from 80 to 110, and from 110 to 20 (Q -1 S varies from 1.3 × 10 -2 to 9.1 × 10 -3 , and from 9.1 × 10 -3 to 5 × 10 -2 ) along the increase of differential stress. The Q evolution obtained from wet data in particular: the attenuation decreases first and then increases, which could be coherent with the velocity evolution. Since the water-saturated sample breaks earlier than the dry sample, the wet sample could have experienced less cracks and therefore caused the difference between the Q evolution in dry and wet experiments. Since we do not know the absolute amplitude of the wavelet input and the conversion formula of the transducer, we cannot discuss the absolute value of Q and hence we discuss Q -1 variation. Figure 3.22b compares the Q -1 variation between different experiments and shows the relative change of

Q -1 /Q -1 0 . The relative change of Q -1 /Q -1
0 shows that the attenuation of saturated rock is stronger than the dry rock. In dry rock,

Q P < Q S (Q -1 P > Q -1 S
) during the compression. On the contrary,

Q P > Q S (Q -1 P < Q -1 S )
during the compression in water-saturated rock. Under the same compressional pressure, the attenuation of the wet sample is higher than a dry sample. In general, the wet sample is more sensitive to the attenuation variation.

Discussion

In this study, we developed a methodology to invert time-lapse active seismic data in the laboratory in order to infer the changes of the elastic and anelastic parameters.

Here, we first discuss the numerical configuration to yield synthetics sufficiently comparable against the observed data. Thence, we resume the observation and make a further interpretation on our experiments and show the need for 3D waveform inversions in the near future. A model with these cut sides implies more extra reflections than a simple cylindrical model. To conclude, it is essential to consider a three-dimensional geometry.

The surrounding silicone rubber layer can provide a buffer-like region, slowing down and reducing the amplitude of quasi-surface-wave propagation between the two media. Note that the resulting wavefield is amplified when the direct waves reach the interface between the two materials and travel through a low-Q and low-velocity material. The transmission coefficient calculated to confirm the observed phenomenon is shown in figure 3.5.

Observed changes in wave velocity and attenuation

The evolution of homogenised velocity and attenuation in the laboratory under a physical state should provide information on microstructural changes. The simultaneous Monte Carlo inversion for four homogenised elastic and viscoelastic parameters, V P , V S , Q P and Q S , was able to eliminate the cross-talks of these parameters.

While Aji granite comprises an average grain size of ∼0.3 mm (Kudo et al., 1992), the effective wavelength can be estimated as around 1.5-3.0 mm using the apparent dominant frequency (∼1-2 MHz) and shear-wave velocity (∼3000 m/s). The local heterogeneity of the grain size can therefore be initially neglected. This simple estimation justifies the homogeneous modelling in this preliminary study, although we would like to refine the structure using our waveform data.

The variation of wave velocity (figure 3.20) during the deformation estimated from our waveform matching methodology has a similar trend as seen by Zaima and Katayama (2018). Hereafter we call the V P and V S measured by hand-picking the "apparent" V P and V S . Our results show slightly lower P-wave and higher S-wave velocities than the apparent velocities. However, the apparent velocity models and the models obtained by our waveform-matching scheme do not correlate at the large stress stage, especially for S-wave acceleration. This indicates that the hand-picking fails even for first arrivals due to the ambiguity of travel times for highly dispersed waveforms. In contrast, waveform matching can still detect objectively the phase arrivals. Figure 3.21 shows the differential V P and V S obtained using waveform matching concerning the initial status for dry and wet experiments. Since the dry sample resists the external stress more than the wet sample, we can observe that the dry sample's velocities drop down to ∼15-18% while the wet sample shows -5% for V P and -12.5% for V S .

Another advantage of this Monte Carlo waveform inversion for the elastic and anelastic parameters is the direct measurements of Q P and Q S . While Toksöz et al. (1979) studied the spectral ratios of waveform records of rock experiment relative to the waveform of a reference sample, Lockner et al. (1977) and Zaima and Katayama (2018) measured the amplitude ratios of the first and second peak to indirectly infer the attenuation, and did not obtain quality factors. Having included Q P and Q S in the simultaneous simulations, objective functions show local maxima, especially for L 1norm and L 2 -norm. In particular, synthetic waveforms ignoring Q P and Q S are not comparable with the highly compressed experimental data. These differences may be due to seismic dispersion induced by attenuation in viscoelastic media, which certainly has more impact on shear waves than on compressional waves. Kjartansson (1979) modelled the relationship between the velocity and the quality factor and obtained the increase of the velocity with frequency and Q -1 . When looking at the discrepancy between S-wave velocities at the last stage of both dry and wet experiments obtained with hand-picking and waveform matching, including attenuation, "apparent" values are constantly higher than our results. This is because of the waveform broadening due to low Q at the last stage.

Since Zaima and Katayama (2018) measured only the amplitude ratios of the first and second peaks to infer the attenuation, we cannot directly compare our results to theirs. Nevertheless, in their experiments, the amplitude of the S-wave is reduced more than that of the P-wave during deformation under saturated wet conditions, which is consistent with our results (figure 3.22). Both P-and S-wave quality factor evolution in wet data show roughly the same trend as the velocity evolution: the velocity increases while Q increases.

Rock physical interpretation

Many studies present a simple interpretation of the general trend of velocity changes.

At low confining pressure, the micro-cracks in the experimental sample can be regarded as randomly distributed. The second stage is the closing process of horizontal micro-cracks, perpendicular to the maximum stress, causing an increase of velocity.

In the final stage, the vertical cracks form before the failure. However, as the closure of the horizontal micro-cracks has little effect on the energy dissipation of P and S waves, it is difficult to see the difference during the attenuation process.

Experimental studies found that seismic wave attenuation in dry rocks was sensitive to intrinsic material anelasticity, thermoelastic effects, friction at grain boundaries and micro-cracks development (e.g., Simmons and Brace, 1965;Walsh, 1965;Gordon and Davis, 1968;Johnston et al., 1979;Winkler and Nur, 1982;Guéguen and Schubnel, 2003;Tisato and Quintal, 2014). On the other hand, water-saturated rocks could accommodate fluid-related mechanisms that would provoke energy dissipation, which should be responsible for seismic attenuation. Winkler and Nur (1982) declared that the effect of frictional sliding on the attenuation is negligible for porefluid saturated rocks. Hence, they emphasised the importance of the role of intercrack fluid flow, which may be considered as the primary cause of the attenuation.

In their dry, partially saturated, and saturated rock experiments, they showed that the S-wave attenuation increases with the degree of saturation and reaches its maximum at total saturation. However, as saturation increases to more than 90% water saturation, Q P is higher than during the condition of partial saturation, and

Q P > Q S .
These observations can be explained with local fluid flow mechanisms. Guéguen and Schubnel (2003) with the observation of Winkler and Nur (1982), who speculated that these effects are due to intercracks flow. At the beginning, increasing pressure causes a negligible impact on both the P-and S-wave attenuation of dry and wet samples, since there is not enough density of cracks that affect the wave propagation. As the effective stress increases, our results show that Q P > Q S in the wet case, whereas Q P < Q S for the dry case. This may be produced when a compressional pressure is high enough to form micro-cracks with higher density, especially with micro-cracks aligned vertically along the long axis. These highly dispersed waves lead to a large amount of attenuation. The characteristics of the local pore fluid flow can explain our attenuation results obtained from a fully-saturated sample (Winkler and Nur, 1982;Guéguen and Schubnel, 2003).

The contribution to attenuation in the dry rock is assumed to be only due to the friction and the intrinsic aggregate attenuation. Intrinsic attenuation is thought to be caused by energy dissipation due to friction at cracks where the faces are barely in contact, making amplitude changes sensitive to the crack geometry in the specimen (Lockner et al., 1977;Walsh, 1966). Besides, Bonner (1974) presented an increase of shear wave anisotropy in the dry granite caused by cracks oriented parallel or obliquely to the compressional stress. The cracks parallel and subparallel to the axis of maximum compression stress influence attenuation. However, we can see that many laboratory experiments confirmed that Q P < Q S is observed in the dried samples (Johnston et al., 1979). No specific mechanism has been proposed to explain the difference between P-and S-wave attenuation in dry rock.

Figure 3.25 compares the exact and relative V P /V S ratio obtained from waveform matching with the V P /V S ratio of Zaima and Katayama (2018). Several studies have mentioned that apparent V P /V S ratio rise can be found in the seismic data of the subduction zone or laboratory experimental data, which may be linked to the high pore fluid pressure and crack anisotropy when approaching failure (Christensen, 1984;Audet et al., 2009;Peacock et al., 2011;Wang et al., 2012). It is worth noting that there is a difference between the velocity variation curve we corrected through attenuation and the apparent V P /V S ratio (figure 3.25b). Our results present higher V P /V S ratios for the wet experiment when rock approaches rock failure. Even in the dry rock experiment, V P /V S ratios are slightly higher than the apparent V P /V S ratios before rupture. These phenomena may be caused by the sudden rise of S-wave anisotropy (Bonner, 1974), thus influencing the attenuation and wave velocity. This characteristic of the V P /V S ratio rising marks the particularity of the rock right before failure and may provide a robust information on fracture prediction.

Toward 3D waveform inversion

Our method can quickly process large amounts of experimental data and provide the best matching values of V P , V S , Q P , and Q S . However, the synthetic seismograms cannot match the later phases regardless of parameter sets. For example, the fifth peak (∼6 µs) of the synthetic waveforms in figure 3.23a cannot match the data, especially when it coincides with the arrival time of S waves. The difficulty may come from model homogenisation that is not valid for scattered S waves. The second possibility is that our boundary condition is still deficient, causing the S-wave waveforms to match poorly. Since S waves and quasi-surface waves are sensitive to smaller heterogeneity than P waves, it is natural to have this poorer fit in S-wave waveforms.

Therefore, the full waveform inverse problem must inevitably move towards devel-oping heterogeneous models. Our method is an isotropic approximation of the deformed rocks. So far, we can only match the first few peaks in the waveform to get the elastic and viscoelastic parameters. Kovalyshen et al. (2018Kovalyshen et al. ( , 2020) ) estimated anisotropy by inverting elastic wave arrival times for characterising the damage and failure processes in rocks. They also discussed the effective interpretation and correction of waveforms recorded by the finite-size transducers in laboratory experiments. When we extract more information on the weakening of rock materials from the waveform, we must use a three-dimensional heterogeneous model. Indeed, it is only through a heterogeneous model that we could better simulate the structural changes and mechanisms inside the rock during deformation.

Conclusion

By matching the simulated and observed data, we obtained the velocity and attenuation changes as the pressure rises up to rupture in the laboratory. Unlike previous studies, the attenuation estimation no longer requires amplitude analysis but can be directly conducted through our simultaneous inversion schemes, without trade-offs between the unknowns. The "apparent" velocity change could be biased by ignorance of attenuation, and therefore our estimation of V P and V S is more reliable.

Combining seismology and rock physics, the data measured can explain the fracture nucleation during varying stress conditions. Seismology can bring useful information on small-scale rock fractures. Variations of elastic and anelastic parameters revealed in this study show different trends from the usual rock measurements; this observation and analysis may thus allow us to understand the rock failure mechanism not only in laboratory experiments, but also in the shallower part of the Earth's crust.

Chapter4

Calibration of active-source seismic laboratory experiments

Introduction

Measurement of rock properties in the laboratory provides us with a macroscopic and statistical view on the dynamic and/or statistic characteristics of rock (e.g., Lockner et al., 1977;Zaima and Katayama, 2018). Laboratory experiments of wave propagation in rock samples are usually carried out to infer the relationship between changes in external conditions and seismic characteristics such as velocity, attenuation, and anisotropy, to further examine the weakening and fracturing process of the rock. The time-lapse relationship between seismic and rock physical properties helps reveal the process of earthquake nucleation. It is also essential in the industry to monitor hydrocarbon and CO 2 storage using the observations of changes in a stress state and seismic features (e.g., [START_REF] Arts | Monitoring of CO2 injected at Sleipner using time-lapse seismic data[END_REF][START_REF] Stork | Assessing the potential to use repeated ambient noise seismic tomography to detect CO 2 leaks: Application to the Aquistore storage site[END_REF][START_REF] Solazzi | Seismic Signatures of Fractured Porous Rocks: The Partially Saturated Case[END_REF][START_REF] Dupuy | Combined geophysical and rock physics workflow for quantitative CO2 monitoring[END_REF]. Therefore, it is necessary to develop numerical models and conduct laboratory experiments to detail the characteristics of rocks under stressed circumstances.

Measurement of seismic parameters of laboratory-scale rock samples is often challenging. A number of authors have estimated the elastic parameters of laboratory-scale samples using the traveltime and phase velocity of active seismic data using transducer sources (Lockner et al., 1977;[START_REF] Lee | High-frequency normal mode propagation in aluminum cylinders[END_REF][START_REF] Zhang | Finite difference modeling of ultrasonic propagation (coda waves) in digital porous cores with un-split convolutional PML and rotated staggered grid[END_REF][START_REF] Shragge | Full-wavefield modeling and reverse time migration of laser ultrasound data: A feasibility study[END_REF][START_REF] Mineo | Surface waves on cylindrical solids: Numerical and experimental study[END_REF][START_REF] Yoshimitsu | Geometric effect on a laboratory-scale wavefield inferred from a three-dimensional numerical simulation[END_REF]Zaima and Katayama, 2018). Based on the relation between the velocity, frequency, and amplitude of the waveform features, [START_REF] Lee | High-frequency normal mode propagation in aluminum cylinders[END_REF] identified the high-amplitude features using a normal mode approach, which showed that the first arrivals depend on the sample geometry and cannot be estimated from the waveforms easily. Shragge et al. ( 2015) simulated the elastic velocity of cylindrical samples observed in ultrasonic laboratory data and developed a 2D inversion to image fracture in a polymethyl methacrylate sample. [START_REF] Yoshimitsu | Geometric effect on a laboratory-scale wavefield inferred from a three-dimensional numerical simulation[END_REF] found the strong conversion of the P-wave to Rayleigh waves at the corners of the cylinder, which implied that only 3D numerical modelling could identify all the wavelets that reproduce the nature of the observed waveform. There have been few studies that quantitatively analysed active seismic waveforms during laboratory experiments using the 3D numerical modelling method (e.g., [START_REF] Lai | Rock deformation monitoring using Monte Carlo waveform inversion[END_REF]. The lack of published studies on conducting laboratory and numerical experiments on relatively small scales under well-controlled conditions can cause discrepancies in interpreting the experimental data. Therefore, there are studies specifically concerning the impact of the size of ultrasonic transducers (Kovalyshen et al., 2020), the effects of the jacketing around the sample [START_REF] Lai | Rock deformation monitoring using Monte Carlo waveform inversion[END_REF], geometry effects [START_REF] Yoshimitsu | Geometric effect on a laboratory-scale wavefield inferred from a three-dimensional numerical simulation[END_REF], which are necessary for obtaining adequate information on experimental waveforms. Laboratory experiments should be conducted first on metal samples to calibrate the instruments' external effects before treating rock samples that are heterogeneous by nature.

We present here the laboratory and numerical experiments in metal and rock samples. We systematically implemented 1) calibrating amplitudes of recorded signals,

2) deriving source time function of the laboratory experiments, 3) wave propagation analysis using numerical modelling on wave propagation and kernels, 4) waveform fitting of single-material metal samples by applying elastic/anelastic waveform inversion, 5) scattering waveform fitting of rock samples by applying Born approx-imation on random heterogeneity models. Numerical modelling was conducted to understand the phases in the waveforms and the effect of the instruments. In this chapter, we start by introducing the calibration of laboratory experiments.

Methodology

The

Error-bar tests

The error-bar tests were introduced to determine the quality of the experimental data and the consistency of the waveforms recorded at different points on the sample. We proposed four different settings for error-bar tests:

• Correlation tests: A acts as the source, B acts as the receiver at the antipodal point, launching ten times.

• Reciprocity tests 1: A acts as the source, B acts as the receiver at the antipodal point, B acts as the source, and A acts as the receiver.

• Reciprocity tests 2: A acts as the source and B acts as the receiver at the antipodal point, A and B switch their positions and launch again.

• Sensor consistency tests: A acts as the source, and B, C, D, E, F sensors act as the receiver at the antipodal point in turn. Figure 4.10 compares the results of correlation tests using sensor T10 as the source and laser as the receiver (green) with sensor T10 as the source and sensor T15 as the receiver (blue). The waveforms were recorded ten times using one kind of parameter setting. Each trace was used as a reference waveform once for implementing crosscorrelation to find the average waveform stability. In general, the correlation among the waveforms using the same parameter setting is good except for 200 mV as FSV and damping factor 1. The correlation of waveforms using sensor T10 and laser and the correlation of waveforms using sensor T10 and sensor T15 are both high.

Source-time function deconvolution

The source wavelet used in the experiments was the same for all the experiments.

The unit and scale of the waveform recorded by the laser are known to be velocity amplitude. However, the unit and scale of the waveform recorded by the piezoelectric transducer were unknown. Therefore, we later used the transducer-laser test data and combined numerical modelling and frequency-domain deconvolution to determine the exact source wavelet.

u syn (ω) = S ricker (ω) • G syn (ω) (4.1)

d obs (ω) = S real (ω) • G real (ω) (4.2) Since G syn ∼ = G real , then d obs (ω) u syn (ω) = S real (ω) S ricker (ω) . (4.3) F -1 F[S ricker (t) * d obs (t)] u syn (ω) = S real (t), (4.4)
where G is green's function, S ricker is Ricker wavelet and S real (t) is the source wavelet of the experiment we want to obtain. u syn and d obs are the amplitude of the waveforms of the numerical synthetics and observed data, respectively.

However, S real obtained from Eq. 4.4 contains both the source time function and the effect of the source transducer. In order to understand the exact source time function without transducer effects, we implemented the following equations:

u t-l (t) = S(t) * G (4.5) u t-t (t) = S(t) * G * R 2 (t) (4.6) S(t) = S (t) * R 1 (t), (4.7)
where u t-l (t) is the waveform of transducer-laser test, u t-t (t) is the waveform of transducer-transducer test and G is the Green's function. S(t) is the effect at the source side, which combining source time function and the transducer effect. R 1 (t)

and R 2 (t) are the effects between receiver transducer and sample at the receiver side.

Assuming R 1 ≈ R 2 , we can obtain S (t) from Eq 4.7, which is the real source wavelet excited by the oscillator without transducer effect.

Waveform amplitude determination

Conclusions

Using numerical modelling of wave propagation, we were able to calibrate the observed waveform amplitude using the source time function of the laboratory experiments by applying frequency-domain deconvolution to the data and synthetics, and comparing the data with the synthetics in different terms in kinematics. It is the very first study that verifies the physical quantity of the measurements recorded by piezo-electric transducers using numerical methods. We investigated the physics of the waveform amplitude, source input, and the effects of the transducers; thus, we can apply it to all the experiments with similar settings. In the previous studies, due to the uncertainty of the source signal, there were many challenges to match the waveforms accurately. We investigated the stability of the waveform measurements by both the transducer-laser and transducer-transducer tests. These results demonstrate that the sensors are the main controlling factor of the waveforms. Removing the specific effects of the sensor itself may be a necessary step when doing sophisticated analysis or using waveforms for imaging. Since using 200 mV-damping 3 and 4 V-damping 1 have the relatively highest values in sensor consistency tests in figure 4.15 and parameter setting of 200 mV-damping 3 had the better results in the previous tests, we chose 200 mV as FSV and the damping 3 (50 Ω) as the setting in the following experiments. Considering the source wavelet and numerical modelling, the measurements we obtain are no longer limited to the first-arrival picking or relative amplitude studies (e.g., Zaima and Katayama, 2018;[START_REF] Lai | Rock deformation monitoring using Monte Carlo waveform inversion[END_REF], but the full waveforms with reliable data. Hereafter, we can obtain the dependable quality factors of the waves propagating in different rock samples. Kovalyshen et al. (2020) studied the finite size of transducers, showing that the transducer size and geometrical configuration have to be accounted for to ensure the accuracy of the modelling. Elastic parameters and the attenuation estimations could be affected by the transducer contacts between another surface or by the nature of complex radiation patterns due to finite size transducers in ultrasonic frequencies. However, in this study, we retain that these effects are included in the source time function obtained from our calibration process. The physical information obtained from the waveforms can offer new insights into the wave propagation and can further be used in all the studies related to laboratory experiments.

Chapter5

Wave propagation in metals

Introduction

Metal samples were used in this study to investigate the differences in the material characteristics and wave propagation (figure 5.1) since their single-material properties are more simple and, therefore, suitable for numerical modelling. All the samples were made and polished into cylinders. An entire cylindrical copper sample was used as the reference, and two aluminium samples with different shapes were used for understanding the geometry effects. The experimental waveforms were analysed using seismic methods. The impact of cylindrical geometry on wave propagation was studied by the snapshots, seismograms and Fréchet kernels. Since the copper sample has a cylindrical shape with a diameter of 4 cm and a height of 10 cm, it was used as the reference for studying the source wavelets for transducer-laser and transducer-transducer tests. Long aluminium, with a diameter of 4 cm and a height of 21 cm, was used in the study to obtain the waveforms with fewer side reflections. Hollow aluminium, which has a hole with a diameter of 3 cm, was used to detect the effects of waves reflecting within a thinner metal layer. The short aluminium with a small hole of 2.5 mm diameter in the middle was used to limit the resolution using our source time functions. These different metal samples with different characteristics can help us distinguish some features in our observed waveforms and test the resolution of our experimental setting and the ability of our numerical simulation methods by matching the data with synthetics.

Samples

Methodology

We implemented active-one-source experiments on different metals. The experimental data were analysed using seismic methods and were compared with numerical synthetics. The elastic parameters of the metal samples were determined using a Monte-Carlo waveform inversion method.

One-source multi-station experiments on different metals

Experimental waveform analysis

Experimental waveforms were analysed to verify the resolution of our two source time functions and characteristics in frequency spectrums. We introduced the 3D numerical models, using the same geometry and elastic modulus as the real metal samples and compared the synthetic waveforms and experimental data. Figure 5.5

shows the raw data and synthetics, which were used to obtain preliminary waveform matching before applying any signal processing.

The structures we want to observe from the waveforms in metal experiments are the thin-layered structure and the small hole of aluminium samples in the second and third samples in figure 5.1. Therefore, we built two models using the same geometry as the two aluminium with holes and obtained the synthetic waveforms propagating inside the models, using the two source time functions of transducer-laser and transducer-transducer tests. The synthetics waveforms propagating in aluminium models with holes were compared with the synthetic waveforms of entire aluminium models without holes using source wavelet of transducer-laser (figure 4 Even though the differences are not noticeable, through numerical modelling, we confirmed that even the tiny holes with a diameter of 2.5 mm could be identified using the two source wavelets. The thin-layered aluminium cylinder with a thickness of 5 mm causes more reflections than the waves to propagate in a full cylinder, which means that we may identify the thin-layer structures simply through the waveform records. of the two aluminium samples with holes. Compared to the aluminium sample with a small hole, the waveforms features of hollow aluminium sample can be better identified by the direct observation of experimental data. The experimental waveform of hollow aluminium contains more reflected waves after the principal wave peak, where the one-centimetre thin layer may cause additional reflections. However, it can also be seen that the higher frequency source time function of the transducer-laser test can better demonstrate the characteristics in waveforms of the small structures. Figure 5.9 shows the limits of the resolution: the synthetics match the experimental waveforms at the first arrivals but failed to detail the structures. To understand the differences between the experimental waveforms and synthetic waveforms, we turned to look at the frequency spectrum of the two types of waveforms. Fast Fourier transform (FFT) was implemented for the synthetic and experimental waveforms of the two aluminium samples collected using one-source one-station transducerlaser (figure 5.10) and transducer-transducer tests (figure 5.11). No matter in which test, the FFT spectrums of experimental data have two sinks around 3×10 5 Hz and 6×10 5 Hz, which are different from the relatively smooth spectrum of the synthetics.

We then further examined the FFT spectrum of the waveforms propagating in other samples (figure 5.12 and 5.13), the two sinks at 3×10 5 Hz and 6×10 5 Hz appear in the waveforms of each sample. Therefore, we think that these particular features may come from the experimental settings, effects of oscillators, effects of the recording machine, or the cylindrical geometry of the samples. In electronics and signal processing, the full scale represents the maximum amplitude a system can represent. Full-scale Volts (FSV) may be defined as the maximum voltage available for the waveform records. However, we had less knowledge of its specific effects on the recorded waveforms. Figure 5.14(a) shows that the main differences for using the different FSV are shown in waveform amplitude. In figure 5.14(b), the FFT spectrums of the experimental waveforms using 100, 200, 300 and 400 mV as FSV have very similar patterns. We can also find that beyond 1×10 5 Hz, all the data waveforms have higher FFT values and lots of noise. However, from figure 5.15, we can see that the information of the signal should lay between 1-2.5 MHz. Therefore, the waveforms collected from these experiments were filtered by a low-pass filter with a cutoff frequency around 2 MHz.

Waveform inversion for elastic parameters

After implementing calibration with numerical tests and recording the waveforms from laboratory experiments, we launched a waveform inversion to investigate the best average elastic parameters model. Since the metal samples have the characteristics of a single material and low attenuation (high Q), we looked for the velocity model for each experimental sample. Also, since the metals have specific elastic properties, their velocity models can be easily determined using the confirmed information in textbooks. We conducted Monte-Carlo waveform inversion by matching the data with synthetic waveforms with simultaneous changes of parameters: P-wave velocity (V P ), S-wave velocity (V S ) [START_REF] Lai | Rock deformation monitoring using Monte Carlo waveform inversion[END_REF]. We used the source time function we obtained from calibration (chapter 4) and set a range of parameters for launching the numerical modelling tests in SEM. By minimising the objective function, the best-fit waveform with average elastic values of the samples was determined and shown in the table 5.1.

Fréchet kernels

In seismic tomography studies, we collect observations of anomalies between recorded and model-predicted seismograms. We then compute the waveform Fréchet derivatives of seismic observables concerning model parameters to invert the changes in observations due to the model perturbation. Fréchet kernels are also known as banana-doughnut traveltime kernels, which are constructed by using the time-reversed wavefield at one particular receiver as the adjoint source. Kernels calculated in this way may be compared with the finite-frequency kernels using ray-based methods (e.g., [START_REF] Hung | Fréchet kernels for finite-frequency traveltimes-ii. examples[END_REF]. The construction of each kernel is based on the interaction between the time-reversed regular field and the adjoint field. We used the P-wave field for constructing the P-wave kernel as in equation :

∂ u i (t,x x x; 0,s s s) ∂ α(q q q) [s m -3 ] = 2ρ 0 (q q q)α 0 (q q q) t 0 dτ τ 0 dτ S k (τ )G pk,p (t -τ -τ ,q q q; 0,s s s)G li,l (τ,q q q; 0,x x x)

(5.1)

where q q q is one random spatial-temporal location in the model, α is P-wave velocity, ∆V is the mesh volume, ρ is density and G l,l is the diagonals of strain Green's function (cf. chapter 2). x x x and s s s represent the location of receiver and source, respectively. There is a doughnut hole if we take the cross-section of a cigar-shaped kernel.

More details for the theory of Fréchet kernels can be found in Zhao and Chevrot (2011). during the experiments. Since our samples were placed on a plastic platform during the experiments, the coupling between the metal samples and the platform may also cause some effects from one side that the numerical simulation using free-surface setting could not reconstruct.

Results

Waveforms and multiphases of waves

Wave propagation modelling in SEM3D offers us a comprehensive understanding of the wavefield variation in the cylindrical model. In this study, we performed two kinds of analysis and succeeded on discovering multiphase of waves: 1) observing the P wavefield snapshots (figure 5.17) and comparing the synthetic (figure 5.18b)

with observed waveforms provide us additional information on phases of waves and 2) by calculating and analysing the kernels, we further verified our observation of waves phases in snapshots and waveforms.

Figure 5.17 shows the cross-section snapshots of P waves propagating in the cylindrical copper model. From t = 5.6 µs and before 8.2 µs, the direct P wave travels from the source point toward the antipodal receiver point within the sample. After t = 8.2 µs, the waves following the P waves, which are later considered as P P and P S waves, appear when travelling along the surface to 45-90 • . Direct P wave and the P P and P S waves merge around t = 9.6 µs. However, at this time, the direct P wave reaches the edge of the long side of the sample, separates into two parts, and travels toward the two sharp ends of the cylindrical model. After getting the antipodes of the source around 10 µs to 18 µs, P S wave travels back to the source point following the same travelling route along the surface. However, at the same time, P P wave travels back to the source point through the interior of the model. After 18 µs, P P and P S waves reach the edge of the long side, which is at the same side of the source point, P P wave separates into two parts and travels toward the two sharp ends of the cylinder. Still, the P S wave travels back toward the antipodes of the source point through the interior of the model. Direct S, S P and S S waves have the same propagation mode as direct P, P P and P S waves, respectively. S waves propagate slower than P waves and closer to each other, making them more challenging to observe and distinguish using only the snapshots. Therefore, we further performed a velocity waveform seismogram by numerically modelling the waveforms at the surface of the cylinder along the cross-section in every 1 • . By tracing the waves using seismogram and using the observation on snapshots above, we identified the phases of the waves in seismograms (figure 5 The propagation modes of P P and S P waves are similar, therefore, the kernels of the P P and S P waves (figure 5.19B(d) and B(g)) are similar as well, so as P S and S S waves (figure 5.19B(c) and B(f)). Since the P-wave kernel provides less sensitivity to S waves, the kernels for the windows of S waves contain less information than the windows of P waves. The kernel of Rayleigh waves shows the characteristics of surface waves travelling along the cylindrical sample's surface. Even in the homogeneous case, Rayleigh waves arise through interactions at the free surface. The strong kernels surrounding the source and receiver, regularly shown in every kernel, represent the scattering waves reflected at the surface.

Conclusion

Waves propagate in metal samples were modeled using 3D SEM modelling method (figure 5.16). Metal samples have high Q values. Therefore, there were fewer constraints for the waveform matching. Since we can trace the phases of waves arriving Figure 5.18b shows the P-wave strain synthetic waveforms, which was recorded densely at every one degree( • ) around the surface containing the source and receiver points. The seismograms exhibited an overlap of different signals, such as direct P and S waves, reflected and converted body waves, and surface waves. Therefore, it is difficult to disentangle one wave phase to another. First of all, combined with snapshots of wavefields (figure 5.17), we identified the phases of waves and further investigated all the laboratory experiments records. Since our 3D numerical simulations have a comprehensive perspective to observe wavefield propagation, we had some additional observations in the wave travelling scheme (figure 5.18a), with respect to previous studies (e.g., [START_REF] Mineo | Surface waves on cylindrical solids: Numerical and experimental study[END_REF][START_REF] Yoshimitsu | Geometric effect on a laboratory-scale wavefield inferred from a three-dimensional numerical simulation[END_REF].

Second, we proceeded to build Fréchet kernels of the waves in the cylindrical sample to confirm our identification of phases (figure 5.19). This is the first study in- Rayleigh wave arrival are evident in both observed and synthetic waveforms. On the contrary, P P , P S , S, S P and S S waves are hardly to be identified. There are two possible reasons: 1) the piezo-electric transducers we used were designed only for P wave or 2) a continuous vibration occurred between the transducer and the attached surface or in the glued area when the direct P wave arrives at the transducer. If the second reason is true, the place where the continuous vibration breaks may be the arrival. The anomalies of the continuous waveforms after the first arrival can be seen in the window of P P , P S and S. This may also give us some idea about removing unwanted information from the observed waveforms and further emphasising the important features.

Even though synthetic and observed waveforms are not perfectly matched, the differences can be easily explained thanks to our numerical methods. However, when modelling the waves propagating in rocks, the grain size, micro-cracks, heterogeneous, anisotropy, etc., all contribute to the attenuation and scattering of the waveforms. The synthetic waveforms can only match the first arrival and the main peaks of data waveforms by applying the Monte-Carlo waveform matching using P-wave velocity (V P ), S-wave velocity (V S ), P-wave quality factor (Q P ) and S-wave quality factor (Q S ) [START_REF] Lai | Rock deformation monitoring using Monte Carlo waveform inversion[END_REF]. Except for elastic/anelastic settings, some additional constraints or physics are needed to optimise the waveform simulation, which is discussed in chapter 6.

Chapter6

Wave propagation in rocks 6.1 Introduction

Rock is a complex mixture of small grains that align randomly or systematically.

Seismic waveforms are sensitive to the different spatial distribution of grains. Due to rocks' various scale heterogeneities, scatters of seismic waves can lead to significant coda features. Coda waveforms, an ensemble of multi-scattered phases observable after main arrival phases, can be a powerful barometer to monitor minor internal structural changes [START_REF] Singh | Coda Wave Interferometry for Accurate Simultaneous Monitoring of Velocity and Acoustic Source Locations in Experimental Rock Physics[END_REF][START_REF] Obermann | 4-D Imaging of Subsurface Changes with Coda Waves: Numerical Studies of 3-D Combined Sensitivity Kernels and Applications to the $$M_{\mathrm{w}}$$7.9, 2008 Wenchuan Earthquake[END_REF]Zotz-Wilson et al., 2020). Random heterogeneity generates incoherent scattered waves. Scattered waves overlapped with the reflected and converted waves could be regarded as noise in the observed waveforms. However, careful observation of these scattered waves will open an avenue to detail the internal structure. In this study, we model the scattering phenomenon numerically using Born approximation. There were several studies modelling scattering phenomena to understand the attenuation and the coda waves of seismic signals [START_REF] Lo | Ultrasonic laboratory tests of geophysical tomographic reconstruction[END_REF][START_REF] Schultz | Enhanced backscattering of seismic waves from a highly irregular, random interface: P-SV case[END_REF][START_REF] Hung | Fréchet kernels for finite-frequency traveltimes-ii. examples[END_REF][START_REF] Sivaji | A physicalmodel study of the statistics of seismic waveform fluctuations in random heterogeneous media[END_REF][START_REF] Spetzler | A test of ray theory and scattering theory based on a laboratory experiment using ultrasonic waves and numerical simulation by finite-difference method[END_REF][START_REF] Angus | A one-way wave equation for modelling seismic waveform variations due to elastic heterogeneity[END_REF][START_REF] Planès | Imaging multiple local changes in heterogeneous media with diffuse waves[END_REF]. [START_REF] Hung | Fréchet kernels for finite-frequency traveltimes-ii. examples[END_REF] used the cross-correlation function to introduce the frequencydepending time-shift 3D body wave tomography. Elastic waves propagating through different laboratory-scaled rock samples with different scales of heterogeneity have been studied (e.g., [START_REF] Sivaji | A physicalmodel study of the statistics of seismic waveform fluctuations in random heterogeneous media[END_REF]. Seismic scattering problems caused by smallscale and random heterogeneity could be partially expressed by Born approximation [START_REF] Hudson | The use of the Born approximation in seismic scattering problems[END_REF][START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition[END_REF]. [START_REF] Angus | A one-way wave equation for modelling seismic waveform variations due to elastic heterogeneity[END_REF] used a stochastic elastic heterogeneous model to simulate frequency-dependent forward-scattering, which showed the appearance of random caustics and the generation of the seismic coda. [START_REF] Lo | Ultrasonic laboratory tests of geophysical tomographic reconstruction[END_REF] used the Born and Rytov approximations for obtaining geophysical diffraction tomography in the vertical seismic profiling (VSP), and cross-borehole experiments. By applying their technique, they found that diffraction tomography is superior to ray tomography because it can image small scatters comparable to the wavelength. [START_REF] Spetzler | A test of ray theory and scattering theory based on a laboratory experiment using ultrasonic waves and numerical simulation by finite-difference method[END_REF] applied ray theory and scattering theory to laboratory and 2D numerical experiments of random media with different small-scale heterogeneity in two types of granite samples. From some of the previous studies above, the evidence was proved that the scattering theory using kernels is more accurate than the ray-path based studies for predicting time-shifts of the small-scale heterogeneous structures. No assumption on wavelengths needs to be made when using the scattering kernels. However, most studies for monitoring weak changes within the earth's subsurface using multiple scattered/coda waveforms were limited to identify relative weak changes rather than understand the exact physics of these changes. Compared to the standard method of measuring the changes in velocity by picking the firstarrival times or coda wave interferometry, full waveforms and scattering theory can further identify the detailed features and extract the physical information from the entire experimental waveforms. However, at the moment, no research systematically considered and analysed at the same time the laboratory experimental setting, source wavelet input, waves propagation, waveform phases, random heterogeneity.

We combined laboratory measurements and numerical modelling to propose a new method for obtaining information on heterogeneity and anisotropy in rocks. In this study, waves propagating in single-mineral samples were employed for implementing data calibration. They were also combined with Fréchet kernel for understanding the nature of waves propagation. Source time function used in the laboratory exper-iments and the scale of the amplitude of the data waveforms were determined using simulated waveforms. The success of calibrating laboratory measurements and obtaining wave-propagating information from experimental waveforms leads us to remove external effects and further serve the fracture monitoring process and our waveform-matching scheme in rocks. We analysed the waves propagating in marble and granite samples, which have different heterogeneities and no preferred-oriented microcracks. Scatters were added to our synthetic waveforms by applying Born approximation on random heterogeneous models with varying sizes of grain and elastic/anelatic properties to approximate the rocks' properties. The predicted grain sizes and anisotropy in rock samples using active-source seismic waveforms are in good accordance with well-known values.

Methodology

Active-source one-station experiment on different rocks

Three rock samples were used to calibrate and investigate the differences in the scattering features (figure 6.1). All the samples were made and polished into cylinders.

The rock samples were with an approximate diameter of 4 cm and a height of 10 cm.

Three types of rocks were experimentally examined: Carrara marble (CM), Indian gabbro (IG) and Westerly granite (WG). These rocks were selected for their characteristics of uniform grain size and approximately homogeneous properties. Apuan Alps (Carrara, Italy). Carrara marble is a pure, crystalline, monometallic and isotropic aggregate of low porosity (<0.5%). This rock, which is a common standard in rock mechanics, is composed of calcite grains (>98%, 150 µm average grain size) and accessory grains of quartz, pyrite, and epidote of few microns in diameter [START_REF] Fredrich | Micromechanics of the brittle to plastic transition in Carrara marble[END_REF][START_REF] Aubry | Fault Stability Across the Seismogenic Zone[END_REF]. Westerly granite (WG) is a finegrained, calc-alkaline granite from Westerly, Rhode Island. This kind of granite, which is composed of potassic feldspar (33%), plagioclase feldspar (33%), quartz (28%), biotite (4%) and muscovite (2%), and a few accessory minerals (1%), is a reference for rock mechanics experiments under crustal conditions. Its grain size ranges from 0.2 to 0.8 mm [START_REF] Wang | Thermal stress cracking in granite[END_REF][START_REF] Aubry | Frictional Heating Processes and Energy Budget During Laboratory Earthquakes[END_REF]. Indian gabbro (IG) consists mainly of clinopyroxene (37%), plagioclase (33%), hornblende (11%), hematite (5%), ilmenite (4.5%), biotite (4%), and a few amounts of orthopyroxene, orthoclase and quartz [START_REF] Hirose | Growth of molten zone as a mechanism of slip weakening of simulated faults in gabbro during frictional melting[END_REF][START_REF] Togo | Stick-slip behavior of Indian gabbro as studied using a NIED large-scale biaxial friction apparatus[END_REF]. The grains in Indian gabbro have no crystallographic or shape-preferred orientation, and their average size is about 0.51 mm.

Elastic/anelastic parameters estimation

After implementing calibration with numerical tests and recording the waveforms from laboratory experiments, we launched a two-steps waveform inversion for investigating: i) the best model of average elastic/anelastic parameters and ii) the most approximate scattering model. First, we looked for each experimental sample's average velocity and attenuation model. We conducted Monte-Carlo waveform inversion

Experimental waveform analysis

Compared with the Fourier spectrum of waveforms of metal experiments in chapter 5.3.2, the frequency spectrum of waveforms propagating in rocks (figure 6.4 and 6.5) has also the two sinks at the same frequency: 3×10 5 Hz and 6×10 5 Hz. This assured that our previous assumption (chapter 5.

Born approximation and Fréchet derivatives

In order to simulate synthetic waveforms for random 3D heterogeneity added to a homogeneous model, we propose to use Born approximation (the first-order approx-imation of Born series). As we derived in chapter 2, we can generate quasi-3D

waveforms by: u i (t,x x x; 0,s s s) = u 0 i (t,x x x; 0,s s s)

+ N p ∑ I=1
∂ u i (t,x x x; 0,s s s) ∂ α(q q q I ) ∆α(q q q I )∆V I (6.1)

where u u u 0 denotes the synthetic waveforms for a homogeneous model and ∂u u u/∂ α(q q q I ) are the partial derivatives for the P-wave velocity at point q q q I . ∆α(q q q I ) m.s -1 is the P-wave velocity perturbation averaged for a small portion of volume ∆V I m 3 in the vicinity of q q q I and N p the number of perturbed voxels. All the detailed derivations can be found in chapter 2. To verify this equation and its limitations, we performed for the copper model with ∆α = 0.01α for merely a small volume ∆V ∼ In this study, we vary the degree of heterogeneity ∆α to mimic the apparent anisotropy and ∆V to mimic the grain size. Note that all the discussions above can be extended to any elastic and anelastic parameters straightforwardly (e.g. [START_REF] Fuji | A methodology for inversion of broadband seismic waveforms for elastic and anelastic structure and its application to the mantle transition zone beneath the Northwestern Pacific[END_REF][START_REF] Zhao | An efficient and flexible approach to the calculation of three-dimensional full-wave Fréchet kernels for seismic tomography-II. Numerical results[END_REF].

As SEM with six Gauss-Lobatto-Legendre points uses eight different mesh sizes, v m (m 3 ) varying from 9.0×10 -13 m 3 to 1.2×10 -11 m 3 . We define g as the number of consecutive meshes that are grouped to have the same ∆α in order to construct random heterogeneity of a volume ∆V = gv m . Figure 6.3 shows the examples of three random heterogeneous models of v m =5.×10 -4 m, 7.×10 -4 m and 9.×10 -4 m with g varying 1 to 100 and -20%≤ ∆α/α 0 ≤ +20%.

Figure 6.12: Time-lapse P-wave kernels related to the middle perturbation point in model A.

Random Heterogeneous models

Random heterogeneous models were designed using two input parameters: g and δ m. Before applying the random heterogeneity to our models, we built eight numerical models with different mesh sizes from 5.0×10 -4 m to 1.2×10 -4 m (figure 6.13).

As described in chapter 6.2.2, g and δ m were then added to the models using Born approximation. δ m represents the maximum of P-wave velocity perturbation in the random heterogeneity models using a continuous uniform distribution over the stated interval from -δ m to δ m. -4 , (e) 8×10 -4 , (f) 9×10 -4 , (g) 1×10 -3 ,(h) 1.1×10 -3 and (i) 1.2×10 -3 m. After constructing the original models, the parameters for building random heterogeneity are later added to these models.

Figure 6.14: Random heterogeneous model V P distribution using δ m =30% as input parameter.

Results

The results of the rock properties discovered by our waveform inversion are presented in two aspects: i) 0D parameters and ii) heterogeneity degree and grain sizes models. The first part of waveform matching is separated into three subsections by three rock samples. The observed and synthetic data of the transducer-laser test and the transducer-transducer test are also discussed in each subsection. In chapter 6.3.4, we show the best heterogeneous models of the three rock samples according to our waveform matching scheme. Furthermore, through these results we can retrieve the rock properties and the grain sizes of the rocks.

Waveform matching: Westerly granite

Transducer-laser test Figures 6.15 and 6.16 show the effects of δ m and g on the waveforms. These effects on amplitudes will be further added to the synthetic waveform of the homogeneous model to model the waveforms with scatters. From figure 6.15, it appears that the parameters of δ m control the amplitude of the scattering effects exclusively. The higher the values of δ m, the higher the amplitudes. However, in figure 6.16, it results that the parameters of g change the phases of the waveforms. In general, the higher the values of g used for applying the scatters in the models, the higher the amplitudes, and the more significant changes in waveform phases. From these observations, δ m and g were necessary to be considered at the same time when building the random heterogeneous models. Experimental waveforms were matched with all the scattering synthetic waveforms with different grain sizes and scales of random heterogeneity using three different objective function: cross-correlation, L1-norm and L2-norm (figures 6.17, 6.22, 6.24, 6.27, 6.29 and 6.33,). The 2D plots of the waveform matching results were performed by different grain sizes, and considering δ m and g at the same time. In figure 6.17, the 2D δ m versus g map illustrates that the waveform matching is mainly dominated by g. Therefore, we considered g preferentially when searching the best-fit models in the 2D images. After fitting the phases of the waveforms well, we searched the best δ m values to match with the amplitudes. Hence, the cross-correlation waveform matching was considered as first, since the cross-correlation can better match the waveform phases. The waveform matching using L1-norm and L2-norm were considered later, since they are more sensitive to detect the differences in amplitudes. Note that the same process of waveform matching of waveforms propagating in Westerly granite was conducted in the following cases of Carrara marble and Indian Gabbro.

Figure 6.18 shows the several best-fit waveforms determined by observing the highest values in figure 6.17(A) and the lowest values in figure 6.17(B) and (C). Among the waveforms in figure 6.18, we chose one of the best-fit grain sizes and tried to search the proper g and δ m. The synthetic waveform with the best combination of the three parameters which contribute to the scattering effects is shown in figure 6.19, where the cost functions point out these parameters in figure 6.17. We picked the g and δ m for models of each mesh size from the waveform matching calculation in figure 6.22 and plotted in figure 6.23. Unlike the previous results of the transducer-laser test, it is difficult for us to tell which waveform matching between the data and synthetic waveforms shown in figure 6.23 is improved by adding the scatters. Therefore, we cannot conclude the best-fit waveform using the scatter method.

Waveform matching: Carrara marble

Transducer-laser test Transducer-transducer test The scatters added by Born approximation contribute to the match of the coda in synthetic waveforms. The improvement to the waveforms can be seen especially in figure 6.35(d) and (f) for matching the coda after the P-wave arrival. The parameters used for modelling the best-fit waveforms were presented in table 6.2. The conversion between the mesh size of the original models and the mesh volume/grain size volume is shown in table 6.3. As shown in figure 6.36, we can see that the matching between the synthetic waveforms with scatters and the experimental waveforms was not improved using the parameters from transducer-laser tests. Furthermore, the scattering effects were found to act less prominently in transducer-transducer experiments (figures 6.22, 6.27 and 6.33). In the previous chapters, we discussed that the addition of the transducer could be considered as adding a low-pass filter to the recording waveform.

Therefore, the non-significant scattering effects may be caused by the low-pass filtered characteristic of the transducer itself, so that it is difficult to find the best-fit waveforms and their parameters by our methods. In addition, the simulated wave- forms before adding scatters can already match the experimental data. While the addition of a transducer prevents us from using waveform matching to predict the grain properties of rock samples and small structural velocity changes, it can provide a good degree of matching over a large region of velocity structures.

Discussion

One common way to recover the elastic parameters in global-scale seismology is using full-wave inversion. However, this approach is expensive and challenging to apply on the structural inversion at grain size scale and highly heterogeneous properties of the rock material. Therefore, first of all, a forward model and inversion algorithm is built to estimate the average elastic parameters that improve the correlation between the synthetic waveforms and the observed waveforms of rock experiments (figure 6.2). Hereafter, we introduced Born approximation with the random elastic heterogeneous model for simulating rocks with different grain sizes and δV P (figure 6.3). In this study, instead of deriving time-shift variations [START_REF] Sivaji | A physicalmodel study of the statistics of seismic waveform fluctuations in random heterogeneous media[END_REF] Table 6.4: Elastic/anelastic parameters and the parameters used for building random heterogeneous rock models [START_REF] Spetzler | A test of ray theory and scattering theory based on a laboratory experiment using ultrasonic waves and numerical simulation by finite-difference method[END_REF], we extended the idea into 3D SEM numerical modelling and took care of the actual amplitudes of waveforms. We used this technique to propose a waveform inversion by introducing different grain sizes and δV P into the numerical models. We expected to predict the grain size, heterogeneous properties, information on anisotropy and crack voids of the target material using only the waveforms.

Model V P (m/s) V S (m/s) ρ (kg/m 3 ) Q P Q S Mesh
Therefore, here in the discussion, we introduced some equations from the rock physical domain to calculate the theoretical elastic characteristics of our rock samples and further compare the values with our modelled results. Furthermore, we calculated the crack density and anisotropy of the rocks and discussed the discrepancy between the theoretical parameters and the ones obtained in this study.

Predicting heterogeneous and grain sizes

Grain size prediction Table 6.4 presents the elastic/anelastic parameters predicted by our modelling and waveform matching scheme. The grain sizes obtained from the former studies and our results are in table 6.5. In the previous studies, based on micrographs, the grain sizes of Westerly granite, Carrara marble and Indian gabbro are all smaller than our results, in particular the grain size of Carrara marble. By observing the scanning electron microscope micrographs of our samples (figures 6.37 and 6.38) and the EBSD observation (figure 6.39), the grain sizes we predicted using the scattering waveform modelling and matching method seem to be suitable and reasonable for interpreting the order of magnitude of the dominant grain size in rocks. Furthermore, refinement of waveform matching methods, such as the use of waveform matching windows, may improve the over-prediction of overall grain size. In the rock material in the real world, the grain sizes vary so drastically that it is difficult to determine a specific grain size value. However, the results obtained by our method can already provide an efficient numerical way with standards to indicate roughly the average grain sizes. Information on heterogeneity obtained from theoretical and numerical modelled V P and δV P When considering the effective moduli of heterogeneous rocks, Voigt-Reuss-Hill (VRH) average is usually used to estimate the elastic constants of polycrystal material and rocks based on their crystal structure and the elastic constants of the constituent crystallites [START_REF] Jaeger | Fundamental of Rock Mechanics[END_REF].

Consider first a hypothetical rock that is homogeneous, consisting of a single min-eral, and without any crack or pore. If a piece of this rock, with volume V , is subjected to a uniform hydrostatic pressure P over its entire outer boundary, a uniform volumetric strain of magnitude P/K will be induced throughout the sample, and the resulting volume decrease will be given by ∆V = PV /K [START_REF] Jaeger | Fundamental of Rock Mechanics[END_REF]. An "effective" bulk modulus K e is defined by the relation K e = PV /∆V . K e represents the bulk modulus of an ideal rock that undergoes the mean volumetric strain as does the actual rock.

Since the total volume change is the sum of the volume changes of the individual minerals:

∆V = N ∑ i=1 PV i K i = N ∑ i=1 PX i V K i = PV N ∑ i=1 X i K i . (6.2)
The Reuss effective bulk modulus (K R e ) is the weighted harmonic mean of the individual bulk moduli:

K R e = PV ∆V = N ∑ i=1 X i K i -1 . (6.3)
On the other hand, Voigt effective bulk modulus (K V e ), taken the definition of K e = PV /∆V = P/ε ν , is presented as the weighted arithmetic mean of individual bulk moduli:

K V e = P ε v = N ∑ i=1 X i K i . (6.4)
Since Voigt's and Reuss's estimates are the upper and lower bounds of the actual effective bulk modulus, the Voigt-Reuss-Hill estimate of effective modulus K V RH e was introduced, using the average of Reuss and Voigt bulk modulus:

K V RH e = 1 2 K R e + K V e , (6.5)
Equation 6.5 can be widely applied to effective shear modulus and other elastic parameters. Since the bulk modulus and P-wave velocity have a quadratic relation, we calculated the effective V Pe and δV Pe from the single-crystal mineral moduli and velocity in table 6.6, using the Voigt-Reuss-Hill (VRH) average. The δV Pe values in the table 6.6 were calculated from the equation 6.6. Since we have K V RH e and G V RH e , Poisson's ratio (ν e ) was also calculated. Also, by our modelled V P and V S , the modelled K and G were determined using equations 6.7 and 6.8. All the effective moduli and velocities calculated by VRH average are described in table 6.7.

δV P = V max P -V min P V mean P
• 100% (6.6)

V P = K + (4/3)G ρ (6.7) V S = G ρ (6.8) ν = 3K -2G 6K + 2G (6.9)
The V Pe calculated from the VRH average are all higher than the V P obtained from our study. δV Pe are all much smaller than the δ m modelled by our scattering waveforms fitting. Since our rock samples are generally isotropic, the discrepancy between the V Pe and V P particularly gives us a glimpse of the influence of anisotropy due to the grain boundaries and the intrinsic cracks. The differences between the V Pe and V P of the rocks samples from high to low are WG > CM > IG. However, since 1) our numerical modelling included the attenuation and 2) the V Pe and V P are the average mean velocity, the absolute values of the difference cannot be used in a quantitative analysis. Our δ m results, which are ±42%, ±96% and ±90%, can be denoted as δV Pm equal to 84%, 192% and 180%, for WG, CM and IG, respectively. The δV Pm values are much larger than δV Pe , which represent the effects of the anisotropy among the grain boundaries and the cracks. Since the rocks have not only heterogeneity but also voids and cracks, these characteristics imply that the elastic velocity disturbance in our numerical model must be largely increased to meet the 0 value caused by the voids. However, Carrara marble has the highest value for both Table 6.6: Elastic moduli and V P of crystals in WG, CM and IG [START_REF] Ahrens | Mineral physics & crystallography: a handbook of physical constants[END_REF][START_REF] Aleksandrov | The elastic properties of rock forming minerals II: Layered silicates[END_REF][START_REF] Liebermann | Magnetoelastic interactions in hematite: Implications for geophysics[END_REF][START_REF] Selvadurai | A multi-phasic approach for estimating the Biot coefficient for Grimsel granite[END_REF][START_REF] Babuska | Anisotropy of vp and vs in rock-forming minerals[END_REF][START_REF] Liebermann | Elastic Constants of Polycrystalline Hematite as a Function of Pressure to 3 Kilobars[END_REF][START_REF] Silva | Ab initio study of the elastic behavior of MgSiO 3 ilmenite at high pressure[END_REF] Rock/Mineral ρ (Mg/m 3 ) K S (GPa) G (GPa) the δV Pm and δV Pe . Compared with WG, CM and IG impact more on the velocity perturbation structure of the rock because of anisotropy or higher density of voids.

V
However, to conclude the crack density calculation, more equations and analyses are needed, shown in the following section.

Crack density prediction

The effective elastic modulus of rocks is usually accurately estimated if their mineralogical composition is known. However, rocks contain cracks and pores, which cause the elastic modulus to decrease below the values obtained in the non-porous matrix. Referring to the {K m , G m } and {K e , G e } in table 6.8, both the bulk and shear modulus calculated from our study are much smaller than the theoretical ones calculated from the elastic modulus of crystals. The difference in the measured velocities has a concomitant variation on the elastic modulus and the Poisson's ratio ν m . Poisson's ratio indicates the ratio of lateral to axial strain in a flexible environment, while the more resistant material has a lower Poisson's ratio. In our results, WG and CM are revealed to have a lower ν m than the theoretical one, whereas IG has an opposite behaviour. Walsh (1965) indicated that the low Poisson's ratio is commonly found in rocks with many cracks. [START_REF] Shearer | Cracked media, Poisson's ratio and the structure of the upper oceanic crust[END_REF] demonstrated that the aspect ratios of damages also influence the Poisson's ratio, so that the Poisson's ratio of rocks with thin cracks of an aspect ratio less than 1/200 increases, whereas it decreases for thick cracks with aspect ratios between 1/20 and 1/2. However, he also mentioned that the spherical cracks cause a negligible effect on Poisson's ratio. According to the results (table 6.8), WG has more cracks or thick cracks. CM and IG have fewer cracks or only thin cracks with small aspect ratios.

In order to conduct the quantification of cracks from modelled V P , V S , we introduced the equations from [START_REF] Fortin | Effects of pore collapse and grain crushing on ultrasonic velocities and Vp/Vs[END_REF], which were used for determining effective moduli of dry rocks with the assumption of spherical pores and penny-shaped cracks:

K e K m = 1 + ρ 1 -p h 1 -2ν e 1 - ν e 2 + p 1 -p 3(1 -ν e ) 2(1 -2ν e ) (6.10)
G e G m = 1 + ρ 1 -p h 1 + ν e 1 - ν e 5 + p 1 -p 15(1 -ν e )
7 -5ν e (6.11)

where {K e , G e } represent the dry bulk and shear moduli of the crack-and pore-free medium, which are the moduli calculated by VRH average from previous mineral studies; {K m , G m } are the elastic modulus calculated from our numerical modelling.

p is the overall porosity, ν e is the ideal Poisson's ratio (equation 6.9), h is a factor that represents the penny-shaped geometry (equation 6.12) and ρ is the crack density we want to obtain. Our samples are non-porous rocks, thus the overall porosity p equals to 0 for WG, CM and IG.

h = 16(1 -ν 2 e ) 9(1 -ν e
2 ) (6.12)

There are three points in our results that need to be noted: 1) the differences between the two crack densities calculated by equation 6.10 and 6.11; 2) the differences between the crack densities among the three rocks; 3) WG has a lower δV P , as men-tioned in the previous section, but has higher crack densities than the other two kinds of rock.

The differences between the two crack densities, ρ K and ρ G , calculated by equation 6.10 and 6.11 may reflect the orientation and shape of the cracks. Guéguen and Schubnel (2003) displayed that 1) the SH waves have no dispersion when propagating through rock with horizontally aligned cracks, 2) the SV waves have no dispersion when propagating through rock with vertically aligned cracks, and 3) P waves disperse more when propagating through vertically aligned cracks but not as drastically as SH waves. ρ G is mostly dominant by V S , therefore the fact that ρ G has a much higher value than ρ K may indicate that there are more vertically aligned cracks inside the rock.

According to the crack densities, our most intuitive interpretation for the three rock samples is that there are more cracks in WG than in CM and IG, which is in agreement with the ν m values. However, as we previously mentioned, WG has a lower δV P but higher crack densities. One possible reason may come from the anisotropy among the partially aligned grains. Westerly granite may have a smaller anisotropy that causes the more negligible effect to the δV Pm , while the Carrara marble and Indian gabbro have higher anisotropy but smaller crack densities. However, the reason may again be the orientation of the cracks, so that the horizontally aligned cracks in WG cause less influence on P waves and as well for δV Pm .

However, the crack density prediction and the elastic modulus shown here are all calculated from the observed waveforms of a single-source single-receiver experimental configuration. If more observed waveforms from different ray paths could be recorded and dedicated to our modelling and analysing process, we can get average crack density values with higher precision. We can find cracks with various distributions in different parts of the sample through the waveforms obtained from multiple experiments at different angles and the crack density calculations.

Conclusion

This study shows that 1) grain sizes, 2) anisotropy, and 3) quantity of cracks inside the rocks can be obtained by matching the experimental waveforms to the synthetic waveforms simulated from our scattering models. Anisotropy of the rock samples can be calculated effectively through our method. Cracks can even be quantified/predicted using only the waveforms. The average V P obtained from our study are slower than the effective average mineral velocities V Pe . The differences between the velocities are caused by the intrinsic micro-cracks inside the rocks or anisotropy among the grains. The additional scattering and the empty spaces lead to the slow average velocities and also the high δ m and δV Pm , which are much larger than the δV Pe . Our research provides an idea of the gap between the theoretical and the predicted values of δV P and offers possible causes. We can even further use our results

to quantify the micro-cracks and identify the orientation of cracks by equations 6.10 and 6.11. By comparing the crack densities with δV Pm and ν m , it becomes possible to distinguish the effects of cracks or anisotropy. Our study could be influential in that we can predict the physical properties of rocks by only sending active-source seismic signals and matching the experimental waveforms with the synthetic waveforms generated by numerical models. The method is highly likely to be further used to predict the rocks' type effectively and further serve the structural imaging during the deforming and fracturing of the rocks.

Chapter7 Migration imaging on metal samples

This chapter presents seismic migration images applied to observed data for metal samples with small holes. Taking benefit of Fréchet derivatives, we perform migration to detect heterogeneities in a cylindrical sample. Although the migration is merely the zeroth iteration of waveform inversion, it is useful to rapidly localise time-lapse differences in deformation experiments that we plan to perform on rock samples. This step is thus essential as a proof-of-concept study for 3D time-lapse waveform inversion for rock deformation experiments in the future.

Introduction

After i) calibrating the experimental instruments, ii) understanding the source time function, the sensitivity and the resolution of different laboratory experiments, iii)

verifying the heterogeneous properties of the samples, and iv) the elastic/anelastic changes inside the rock samples during the deformation, we have built the basics for the ultimate aim of this study: time-lapse imaging of rock samples with a small number of receivers. The next step in this study is to test models with various structures and image the irregular structure to prepare for the future final goal: monitoring the time-lapse rupture (figure 7.1). In this chapter, we first introduce several transducer arrays to detect overall 3D internal structure changes. We designed various local heterogeneity to implement numerical tests using different transducer arrays after verifying the proper transducer array for imaging inner structures in the samples. We proceed to the numerical and laboratory experiments on metal samples with previously known internal structures. The homogeneous and locally heterogeneous metal samples were used in the laboratory experiments, and the waveforms propagating in these two kinds of samples were analysed simultaneously. Seismic migration was performed on the previously known structure to examine the resolution and accuracy of our imaging method and the experimental setting. The tests presented in this chapter were designed as the last milestone to systematically establish the practical link between laboratory experiments and 3D numerical simulations. Hereafter, we can freely apply this system to all the laboratory experiments on metals or rocks for effectively analysing the waveforms using seismic methods and inverting internal structures using full waveforms.

Methodology

Forward numerical modelling was first implemented to build various local heterogeneous models in SEM3D. The models were established to obtain the influence on waveforms when adding irregular local heterogeneity. The local heterogeneous models are expected to combine with the random heterogeneous models and form a model with global scatters and local anomalies for future studies. Despite the numerical experiments, active-source multi-stations experiments were designed to observe the waveforms recorded all around in a 3D structure, verify our imaging methods, create a sophisticated laboratory signal processing method to optimise inverted images, etc. The numerical tests using the configuration of the laboratory experiments were also implemented in advance to evaluate the feasibility before performing experiments. The seismic migration imaging technique can be considered a way for "migrating" back the reflection events to their proper locations within the medium or in the earth, forming an accurate image for the target structure. In this study, we want to make a detailed use of the amplitude and phase information in the migrated image derived from laboratory experiments in rocks. In this chapter, we describe the mechanics of migration and some of the problems related to small-scale laboratory experiments, such as accuracy and efficiency, velocity estimation and the issues related to numerical modelling.

Local heterogeneity

Defining a region of local heterogeneity in numerical models of SEM3D is technically difficult. We started from designing simple and relatively large structures inside cylindrical models (figure 7.2 and 7.3). The random separated (figure 7.4) and partially combined small structures (figure 7.5) were then investigated. The models with faults in early developmental stages were then built successfully in SEM3D (figure 7.6 and 7.7). The forward modelled synthetic waveforms obtained using the fault models can also give us insight into the characteristics of waves propagating through a preexisting fault. These numerical models were further used for migration imaging to test the capability of the imaging method.

Multi-station laboratory experiments

Two main issues emphasised in this section are 1) the use of previously known structures in laboratory experiments and 2) multi-station configuration. Other instrument settings relevant to laboratory experiments are consistent with the optimal settings discussed in previous chapters.

Active-source multi-transducer experiments and samples with structures A homogeneous aluminium cylinder was used as the reference. The copper part inside the aluminium-copper concentric sample and the hollow part of the recessed aluminium with small and big holes were considered the target structure needed to be distinguished by the waveforms. The small hole and the copper cylinder inside the aluminium-copper sample have a diameter of 5 mm, while the large hole has a diameter of 10 mm. Up to 24 transducers (figure 4.4) were available to be placed on the cylindrical samples, where the attaching surfaces of the transducers were designed to be curved and fit perfectly on the surface of a cylinder. As shown in figure 4.4 and 7.9, all the transducers were numbered, and the positions where we placed the transducers were also marked. Therefore, we can confirm that the factors that affect a signal are the transducers and the connecting point. The same transducer was placed on all the other samples at the same position as the original reference sample. Each transducer acted as the source while the other transducers served as receivers.

The homogeneous reference aluminium and the aluminium samples with the inner structures were subject to the same experiments using the same configuration as figure 7.10.

Transducer arrays

Monitoring time-lapse 3D structural changes is challenging. This is especially true in a centimetre-scaled small sample as we used in this study, because of their relatively large size. Removing the transducer effects or avoiding some transducers being too characteristic to affect the imaging results is the main task of applying this small-scaled laboratory experiment. The configuration of the transducers also plays a key role for the 3D imaging. [START_REF] Brantut | Time-resolved tomography using acoustic emissions in the laboratory, and application to sandstone compaction[END_REF] applied a transducer array with six- Figure 7.12: Twenty-four stations setting for initial test in our compression experiment in laboratory.

Seismic migration

Seismic migration is the zeroth iteration (gradient direction of) full-waveform inversion process (e.g., [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF]) that uses wave equations to "migrates" distortions in waveform records. It is achieved by shifting events to their real spatial locations and bringing diffracted energy back to its original scattering point.

Reverse-time migration is a standard seismic migration scheme that solves the full acoustic/elastic wave equation by extrapolating time and allowing the waves to propagate in all directions. [START_REF] Baysal | Reverse time migration[END_REF] and [START_REF] Mcmechan | Migration by extrapolation of time-dependent boundary values[END_REF] present the method in detail and describe its ability to image all dips with sufficient accuracy, since the wave equation can accurately simulate wave propagation in all directions, including reflections and transmissions. Therefore, this algorithm is robust for velocity modelling with complex geometry. The principal difference of the reverse-time migration from the other migration schemes is the sending back process of backward wavefield recorded at the surface using the wave equation with time running backwards [START_REF] Gray | Seismic migration problems and solutions[END_REF]. Here in this study, we use equation 2.29 that we repeat here:

M (q q q) ∝ ∑ i ∑ r r r ∑ s s s t 0
dt [d i (t,r r r; 0,s s s)u i (t,r r r; 0,s s s)] ∂ u i (t,x x x; 0,s s s) ∂ α(q q q) (7.1)

where M (q q q) is the P-wave migration image at q q q, which is the convolution between the waveform residual between observed data d d d and synthetic data u u u for a reference model and the Fréchet derivatives ∂u u u/∂ α(q q q). Note that we make summation on all the receivers r r r and sources s s s. Detailed discussion can be found in chapter 2.

The imaging scheme defined in equation 7.1 usually generates artefacts with large amplitude near the sharp interfaces, where these noises can even ultimately contaminate the signals. [START_REF] Liu | An effective imaging condition for reverse-time migration using wavefield decomposition[END_REF] discussed the noises generated using equation 7.1 in detail and proposed a new imaging condition to improve the reverse-time migration imaging scheme. In this work, since we have sharp surfaces at the sides and use a relatively small-scaled sample compared to our sensors, we may need to carefully take care of these numerical artefacts for our imaging scheme.

Results

Waveform analysis

Figure 7.13: Experimental waveforms on aluminium cylinder and aluminiumcopper cylinder, traveled from a source transducer (T10) in the middle layer and recorded at the antipodal point by laser using different FSV and damping setting. Even though we considered the combination of FSV=200 mV and damping equal to 3 (50Ω) more proper to be used in our study, we still tried all the pairs of FSV and damping again to see if there is any difference when using homogeneous and nonhomogeneous samples. Surprisingly, the variation of FSV affects the relative waveform amplitudes of aluminium and aluminium-copper samples drastically. When normalizing the waveforms in figure 7.13 as figure 7.14, the experimental waveforms of the two samples appear to be less different in general. Only when plotting the waveforms separately and zooming in to a certain smaller time window as fig-

ure 7.15, some small differences in phases between the waveforms recorded from two samples can be observed.

Since we were taking FSV = 200 mV and damping = 3 (50Ω) as the preferred setting as we discussed in the previous chapter 4.3, only the waveforms with the preferred configurations were used in the following content of waveform analysis. In fig-

ure 7.16, the waveforms of transducer-transducer tests on aluminium and aluminiumcopper samples were plotted, while the source transducer was T10 at position P09 (figure 7.12) and the other 23 transducers were implemented as receivers. The labels in figure 7.16 have two parts, so that "T" represents the transducer number and "P" means the position of the transducer as the points marked in figure 7.12. The waveform amplitudes of the aluminium cylinder recorded at the same middle layer (P09-P16) as transducer T10 are significantly higher than the waveforms recorded at the top (P01-P08) and bottom layer (P17-P24). However, when a thin copper cylinder is inside the aluminium cylinder, blocking the wave propagation direction, the waveforms recorded at the three layers have no apparent differences. It may be related to the relatively low-velocity zone of the copper cylinder (table 5.1) that traps the waves. Since the first arrivals of the waves have no necessity to pass through the thin copper cylinder in the middle, these phenomena are not evident at the other layers. However, the reduction of waveform amplitude caused by the copper may be seen at the later phases of the waveforms.

We then implemented time-frequency analysis on the transducer-laser waveforms, which used the preferred experimental settings, on reference aluminium cylinder and aluminium-copper cylinder in figure 7.13. It is difficult to identify the difference between the spectrograms of waves propagating in reference aluminium (figure 7.17)

or aluminium-copper cylinder (figure 7.18). Therefore, time-frequency misfit functions based on [START_REF] Kristeková | Time-frequency misfit and goodness-of-fit criteria for quantitative comparison of time signals[END_REF] Time-frequency misfit between the the transducer(T10)transducer(T23) waveforms of aluminium and aluminium-copper cylinder.

Local heterogeneity migration imaging using numerical synthetics

Migration imaging was tested in all the local heterogeneous models presented in the previous section 7.2.1. We present an example of the results of migration imaging, which used the numerical models of figure 7.2. Figure 7.22 and 7.23 are the cross-sections of the 3D migration imaging results at 48 mm and 50 mm of the cylindrical model, respectively. The dark blues area perfectly restores the velocity structure anomaly of the target model. However, our small-scale experiments are unavoidable flower-like numerical artefacts caused by the transducer array distribution. This can be a part of the prospects that can be improved.

Discussion and conclusion

This chapter performs an ongoing project that sums up all previous numerical and experimental experiences and advances to a new attempt in this field. Micro-tomography is not a new idea to be applied in the rock samples under pressure. However, it is the first research focused to derive every step in the laboratory experiments in detail using numerical and seismic methods. The previous unknown observation, some instrument or sample effects that cause uninterpretable and unexplainable impact will also have the opportunity to be excluded from now on. Therefore, the microtomography built by migration imaging is highly likely to provide high value in rock physics and seismic nucleation mechanisms. In this chapter, we identified ways to make weird velocity anomalies in 3D models, create new optimal 3D experimental configurations, analyse the effects of intra-sample structures on the signal, illustrate the impact of the detector array on the signal and imaging and implement the simulation and a preliminary migration imaging results. Even though the results we got in our numerical simulation experiments are positive, there are still many improvements and many possibilities in this research. The combination of numerical modelling and seismic analytical methods adds many new information and challenges to these laboratory data, contributing to future studies of imaging the structures of rock samples.

Chapter8

Discussion and conclusions This is the first study combining the small-scale laboratory experiments, numerical modelling methods and seismic waveform inversion altogether, and forming them into a new systematical process for dealing with the data collected in the lab. This research brings novel aspects to traditional rock physical studies. We demonstrated that the experimental data could be further applied in many seismological ways, so that there is more information that can be extracted from the waveforms. Since we have discussed theory and results in each chapter, here, we focus on the most significant achievements and discuss some future perspectives.

Significant achievements

Although several previous studies modelled the rupture process in the seismology field and applied their laboratory experiences back to the field data, the recorded experimental waveforms have never been studied and appropriately treated. Therefore, we started from the very beginning to implement calibration and seismic analysis to understand our experimental data. We simulated synthetic waveforms and Fréchet kernels to understand the wave propagation scheme and the effects of the geometry. We successfully matched the experimental data and the numerically modelled waveforms from a database to invert the time-lapse average elastic/anelastic struc-tural changes during the rock deformation. This is the first practical study that brings the attenuation into the view, combining the experimental and numerical way. After building a system for effectively making velocity and attenuation models for conducting inversion, we advanced our models and inversion by adding scatters using Born approximation to improve the fit of synthetics with the waveform coda in rock experiments. The grain sizes, anisotropy and other properties of rock samples can be predicted by matching the observed and synthetic waveforms. The results shed light on a considerable probability to predict the rock types by sending waves and analysing the seismic signal. The numerical and experimental configuration for applying migration imaging on small-scale samples were investigated and successfully conducted in numerical cases. The initial use of actual experimental data on migration imaging has also been started. This study builds a bridge connecting the two fields: imaging seismology and rock mechanics. We expect to broaden this research into more experimental works and field applications.

Ideas for future investigations

This thesis provides an up-to-date methodology combining rock physics and seismology on a small scale. It may also dedicate to a more considerable scaled seismological work for earthquakes mechanics in the future. In this study, we see several future directions that have not yet been achieved. More efforts are needed to image the time-lapse structural variation in rocks. With the knowledge we have acquired through this research, we discovered many possibilities for significant achievements and discoveries in the future. We, therefore, provide four future outlooks that are likely to be implemented and quite important after this thesis.

Topic A: Optimizing migration imaging process

In this study, one of the unsolved questions is to remove the numerical effects in migration imaging and optimise the imaging process. Applying some filters and other signal processing methods to the experimental waveforms reflects on our synthetic numerical results. cross-section. Then, the transducer array with twenty-four stations (figure 7.12) can be used to image the rupture process during the compression experiments.

Topic C: Full waveform inversion on rock samples

Migration imaging is a kind of inversion without iterations. Iteration may deteriorate the inverted image in the small-scale study. However, implementing only one or two iterations may be a possible way to balance and find the best solution for the imaging. These models containing a discontinuous surface may give us some ideas of waveform variation from one medium to another different medium. With the help of our scattering background model in this study and the local heterogeneous model, we may better understand the subtle changes in the propagation of waves from one rock to another, and thus be able to predict in which type of rock the waves propagate through.

Introduction

Full waveform analysis of active data during deformation experiments could provide more localised 2-D details than those based only on travel time measurements. The detailed time-lapse 2-D images will lead to new insights into microscopic dynamics during weakening and rupture processes. In order to constrain the physics of the weakening and rupture processes of initially unfractured rock in a laboratory scale, and to better understand the earthquake nucleation, we perform waveform modelling to infer to the time-lapse elastic structure changes during laboratory experiments focusing on full waveform details.

Many studies intend to start from finding general physical characteristics from laboratory experimental data. Lockner et al. [1977] obtained the changes of seismic velocity models and attenuation during the deformation by calculating the traveltime. Zaima and Katayama [2018] applied triaxial deformation experiment on a centimetric-scale granite from Aji region, Japan, under dry and wet condition. They found that as an axial pressure applied on the top of the sample increases and the rock sample deforms gradually, the time-lapse data obtained from the receiver show the trend of the velocity change. However, from Lockner et al. [1977] to Zaima and Katayama [2018], they only used P-and S-wave traveltimes picked only for one station, providing averaged wavespeeds of the rock with little knowledge on the temporal stress distribution changes. These waveforms can provide us more information on structural changes than only velocity changing. Thus, in this study, we used the experimental data of Zaima and Katayama [2018] and full-waveform modelling to match the observed waveforms to obtain the velocity models and the waveforms changes during the deformation of the rock. We build a series of velocity models with response estimation and implement an automatic Monte Carlo waveform fitting method to obtain the best-fit velocity models. In addition, we also performed many tests for lateral heterogeneity.

We have tried random models, fault models, etc., to successfully simulate the different distributions of wavefields and the waveforms. From these results, we can construct a rough process of weakening of the rock sample to help understanding the physical characteristics of the earthquake nucleation.

Method and Theory

In order to calculate elastic wavefield in a medium with a curve discontinuities or free surface boundaries, we can use a tensorial formulation of the wave equation and curvilinear coordinate transformations to make the computational uniform grid coincide with the curved interface in physical domain (e.g. Komatitsch and Tromp [1999]). We used finite-difference on collocated Cartesian introducing zig-zag free-surface discontinuities. However, we modified the boundary conditions for optimally accurate collocated finite-difference method based on predictor-corrector scheme (Takeuchi and Geller [2000], Fuji et al. [2016]) for 2-D P-SV medium in this study.

We tried to perform the modelling in the same geometry and conditions as the experiments. The parameters used are dt=0.002 µs, dx=0.005 cm, and dz=0.005 cm. The source wavelet is sent as the real input source used in the experiments with 2 MHz resonant frequency in compressional and shear axes. Since the source is on the side of the cylinder, we put the point source on the top and put one receiver at the bottom of the 2-D circle (Fig. 1A). Here we can use the reciprocity theorem to model synthetic waveforms for the sources of two kinds: compressional and shear, just by using a single force from the receiver to the antipode source position and calculate the compressional and shear strains. We start with homogeneous model since the Aji granite is reported to be nearly isotropic texture with a small compressional wave anisotropy (< 1.3%) and shear wave velocity has no significant changes in different polarisation directions. Homogeneous velocity models are tested in different P-wave and S-wave velocities in order to find a best-fit model for the experimental waveform results during the deformation.

Numerical Examples

We take the velocity model from Zaima and Katayama [2018], comparing it with the observed waveform. It obviously shows the shortcomings of using hand-picking method with only first arrival to infer to the velocity models. We use full waveform for data fitting instead. However, the waveform misfit is too large to go beyond the first arrival fitting for the moment and thus we further explored the way to better model the waveforms. There are many reasons that will cause the differences between synthetics and observed data. In our case of 2-D circle disk forward modelling, 3-D effect, instrumental response, attenuation, finite-surface source, and other unknowns, may cause a huge discrepancy between observed and synthetic data. Hence, we design the response function estimation, which involves all the impact factors into the function, in order to rebuild the similar synthetics to the observed waveform. Since we have:

ℱ #$ % ℱ(' 𝒎) *+, )(.) ℱ(' 𝒎)
,/0 )(.) • ℱ(𝑢 𝒎$ 345 (𝑡))7 = 𝑢 𝒎$ 9:3 (𝑡)

we can have the response of the experiment by taking the relationship between the observed waveform and the synthetic waveform of the same model (m1). Even though the time-lapse P-wave velocity can be slightly changed due to the first setting of P-wave velocity (m1), it will not be as slow as Vp based on hand-picked traveltimes and we value more about the ∆Vp but not the absolute Vp in this study.

Then, we can get the response function R(ω) as follows:

ℛ(𝜔) = ℱ(𝑢 𝒎$ 9:3 )(𝑡) • $ ℱ(' 𝒎) ,/0 )(.) (2)
We can rebuild the waveform similar to the experimental results from only using the numerical synthetics by doing convolution with the response function: ℛ(𝑡) * 𝑢 𝒎A 345 (𝑡) ≅ 𝑢 𝒎A 9:3 (𝑡) We test several norms to automatically and accurately detect the time-lapse structure changes (Fig. 1B). Therefore, we first comply with the eyeball-based trial-and-error method in order to find the best values of velocity models. We then introduce automated Monte-Carlo methods, which are based on waveform misfit, to construct a system for finding the models fast and accurately. We processed a series of velocity with an interval of ∆v = 10 m/s that we have 191 velocity models for autofitting. The automated Monte-Carlo system will chose the synthetic waveform with the smallest error relative to the observed waveform, and chose this model of synthetic data as the best model.

In order that the synthetic waveforms be comparable to observed data, we use a response function (Eq.

3) that contains all the information other than the structural information that we would like to extract (Fig. 2A). The response function is taken from the relation between the zero-displacement (zero-timelapse) waveform traces of observed data and synthetics. We truncate the waveforms before shear waves around in from 3 to 6 microseconds, and use it for response function. Applying the same response function to all the waveforms we have, we can obtain the synthetic waveform with the proper transfer function (Fig. 2B).

In order to be able to investigate the effect of heterogeneities on the waveforms, we perform random heterogeneous media (Fig. 3). After many tests for random heterogeneous media, we are convinced that it must be able to help a lot with understanding the structure of the sample during the deformation. For the purpose of modelling the wave propagating with a fault inside the media (Fig. 4), we created the average media with a fault in it with different width to see the waveform changes during the fault development. We perform both the eye-ball based trial-and-error and L1 norm automated Monte-Carlo fitting to find the best-fit homogeneous velocity models for the three fault models. We superimpose the synthetic waveform for homogeneous models on that of fault models to see the difference. In Fig. 4B, we can see the different feature of the waveforms of the fault models to the homogeneous model that we are convinced that the full waveform can give a clue to the structure inside the sample.

Conclusions

We developed the full-waveform-based elastic imaging methods for laboratory rock experiments. We used a finite-difference code on collocated Cartesian coordinate with zig-zag free-surface discontinuities for boundary for the 2-D circle forward modelling. After comparison, we know that the response function estimation can successfully rebuild the waveform before the attenuation influences by considering all the effects, such as 3-D effect, instrumental response, attenuation and etc. After testing, we obtained the best automated fitting method-L1 norm. With response function and eyeballbased waveform fitting and automated Monte-Carlo waveform fitting, we obtained the averaged timelapse P-wave velocities based on observed waveforms.

In faults modelling test, we know that the existence and extension of faults will lead to the decreasing of the velocity models. The systematic changes in elastic wave velocities in Aji granite during triaxial deformation experiment under dry and wet conditions are related to the development and closure of cracks. In the early stages of deformation, both Vp and Vs increase slightly due to the closure of preexisting cracks that are mostly oriented normal to the direction of maximum compression stress. Elastic wave velocities then decrease as the sample approaches failure, due to the growth of cracks aligned subparallel to the axis of compression. From the confirmation of real data fitting and faults modelling in this study, we can tell that with the wave propagation information and waveform fitting technique, the physical characteristics of the earthquake generation can be inferred. Using full waveforms, this project will further detail a microscopic weakening process during the nucleation of earthquakes.

The large discrepancy between observed waveform and synthetic before convolution with response function showed that it is necessary to take account for more realistic configuration and to design the laboratory experiments by changing the size of the sample, experimental geometry, and etc. Only by extending the study to the 3-D dimension can we understand the weaknesses and rupturing that may occur in real conditions. We look forward to using this method, the new laboratory experiment design, combined with the finite-element method, to perfectly model the actual situation of the earthquake nucleation and learn about the physical characteristics of the earthquake.

Introduction

In order to understand the earthquake nucleation, it is indispensable to study both macroscopic and microscopic behaviors of the rupture process of the Earth's crust and mantle: by global-scale numerical modeling (Fliss et al., 2005;Gabriel et al., 2012); natural earthquakes observation (Di Carli et al., 2010;Lee et al., 2011); and laboratory experiments of rock samples (Benson et al., 2006;Guéguen & Schubnel, 2003;Lockner et al., 1977;Toksöz et al., 1979). Since averaged elastic wave velocities of a loaded rock sample could reflect the localized nucleation and/or growth of cracks, it is important to understand time-lapse changes of an experimental rock during its deformation to get insights on the macroscopic earthquake generation process (Scholz et al., 1973).

There have been several efforts to seismologically image the interior of rock samples. Möllhoff et al. (2010) simulated the preexisting artificial fracture in the rock using 2D discrete elastic lattice method, whereas Lai et al. ( 2019) studied the velocity model changes during rock deformation using the 2D finite differences method. However, 2D numerical modeling is insufficient to describe wave propagation in a 3D rock, even though the transfer function of Lai et al. (2019) or a 2D line source approximation (Igel et al., 2002) seems to be functioning to some extent. In order to surpass the limits, there have been some promising recent efforts of 3D numerical modeling of seismic waveforms in a laboratory scale, albeit a number of difficulties. Yoshimitsu et al. (2016), for instance, modeled the geometry effect on cylindrical aluminum using the 3D finite difference method, for a cylinder with 50 mm in diameter and 100 mm in height with a dominant frequency between 200 and 800 kHz. Their comparisons between synthetics and observed data, being limited to homogeneous metal and plastic experiments, have shown the feasibility of the numerical reproducibility of active seismic data in laboratory. In addition to wave propagation simulation, Brantut (2018) monitored acoustic emission on small-scaled sandstone using a 3D fast-marching method with finite differences and

Abstract

We estimate elastic and anelastic parameters and their evolution during laboratory rock deformation experiments, while developing a Monte Carlo waveform inversion. The transducer-totransducer one-source one-station active seismic data of dry and water-saturated samples are obtained from Zaima and Katayama (2018), https://doi.org/10.1029/2018JB016377. We first performed a trial-anderror estimate of the boundary conditions in order to suppress its influence on waveforms. The synthetic seismic data were generated using equivalent homogeneous models with different combinations of elastic and anelastic parameters with the aid of spectral element method. We compared them with the laboratory experimental data. Based on the comparisons, we obtained the time-lapse variations of seismic velocities and attenuation of rock samples during deformation experiments, which we interpreted as crack developments. Our simultaneous estimation of elastic and anelastic parameters allowed us to detail the dynamics prior to the rock failure.

Plain Language Summary

We developed a novel methodology to monitor the rock deformation in the laboratory using active seismic waveforms, comparing them against 3D synthetic seismograms. Our methodology can invert for the P-and S-wave velocity changes as well as their attenuation quality factors. Changes in velocity and attenuation show the different trends especially just before rock failure, indicating the complexity of microcracks' evolution.

LAI ET AL.

LAI ET AL. 2020) imaged the development of fractures using migration of active seismic data, which they called coda wave decorrelation inversion. The presented studies are encouraging but some of the numerical treatments of the present complex problem should be fostered. The boundary condition employed in finite differences method is unrealistic (Komatitsch & Vilotte, 1998), and the kernels used in the migration are too simplified using generalized ray theory. The aim of the present work is to extend the idea, combine with 3D numerical modeling on wave propagation, and carefully treat the entire waveforms in order to extract information on elastic and anelastic parameters changes.

Here, we would like to develop a methodology which inverts the seismic waveforms for structural changes during laboratory experiments performed by Zaima and Katayama (2018). They performed a triaxial compression experiment on the Aji granite and performed active seismic measurements during the rock deformation. The idea is to robustly and quantitatively infer the subtle structural changes of anelasticity as well as elasticity, and three-dimensional structure in the near future. There are little number of previous studies that attempted "waveform inversion" of laboratory experiments. We use the spectral element method (SEM: Komatitsch & Vilotte, 1998;Chaljub et al., 2007) to model wave propagation in a rock sample during deformation. The spectral element method can accurately handle the boundary condition in 3D complex models (Capdeville et al., 2003;Chaljub et al., 2015;De Basabe & Sen, 2014;Peter et al., 2011). However, due to its computational costs, there are only few studies devoted to the comparison of SEM synthetics with laboratory experimental data (e.g., Pageot et al., 2017;Solymosi et al., 2018). Here, we perform SEM modeling on one-source and one-receiver experimental data, in order to show the adaptability of SEM to the seismic data in laboratory. Throughout 3D simulation, we systematically construct a database of waveforms for different sets of elastic and anelastic parameters: P-wave velocity ( P E V ), S-wave velocity ( S E V ), P-wave quality factor ( P E Q ) and S-wave quality factor ( S E Q ). We then perform Monte Carlo waveform inversion using an objective function which combines three kinds of cost functions: 1 E L -norm and 2 E L -norm of waveform misfit and zero-lag cross correlation. We would like to explore information on rock elastic properties and attenuation changes using numerical simulation and Monte Carlo inversion. Throughout our approach, we will show how numerical models help to understand the observations and to illustrate seismic characteristics of rock deformation.

Data and Methods

In this study, we perform inversions to image elastic and anelastic structural changes during the rock deformation. Hence, we first acquire the active seismic data during the deformation experiment. Meanwhile, we prepare a set of synthetic seismograms with the aid of spectral element method. We then select the preferred models for observed seismograms (compressional and shear strains) for each differential stress step, in order to obtain the time-lapse trend of elastic and anelastic parameters. First, we revisit briefly the experiment configurations and provide known and unknown properties in the experiments. Second, we describe the numerical modeling scheme. Finally, we present our data processing strategy based on the data and modeling conditions and detail the inverse problem implemented in this study.

Experimental Setup

Triaxial compression experiments apply a confining pressure to a cylindrical rock sample wrapped in an impervious membrane, and then loads it axially to failure with dynamic compression. Zaima and Katayama (2018) conducted triaxial compression experiment on the Aji granite, which mainly contains quartz, plagioclase, K-feldspar and biotite, with an average grain size of about 0.3 mm (Kudo et al., 1992). Both the dry and wet conditions were applied during the experiments. The wet experiment was under fluid-saturated condition, where distilled water was used as a pore fluid. The pore pressure of wet experiment was kept constant at 10 MPa during the deformation. The experiment was performed at room temperature with a confining pressure of 20 MPa and a constant strain rate of 10.1029/2021JB021873

3 of 23 source and a receiver, respectively. P waves and S waves were generated and recorded separately using piezoelectric transducers located at the symmetrical sides of the center point at the edge (Figure 1a). For the purpose of gluing the transducers, the two opposite sides of the sample were cut 0.6 mm off and polished, as shown in Figure 1b. Each P-and S-wave waveform with increasing differential stress was obtained by stacking signals for each compression condition during the rock deformation until the failure of the rock sample. When the compression stress approaches the critical state, which varies from one rock sample to another, the rock breaks. All the experiments data and details can be found in Zaima and Katayama (2018).

Forward Modeling of 3D Elastic Waves Using Spectral Element Method

In our work, we use the spectral element code SEM3D [START_REF] Cupillard | RegSEM: A versatile code based on the Spectral Element Method to compute seismic wave propagation at the regional scale[END_REF][START_REF] Delavaud | Simulation numérique de la propagation d'ondes en milieu géologique complexe: Application à l'évaluation de la réponse sismique du bassin de caracas (Venezuela)[END_REF] to compute synthetic seismograms traveling through the cylindrical rock samples. Here, we briefly review the spectral element method but the readers can refer to Komatitsch and Vilotte (1998) and Chaljub et al. (2007). The SEM is based upon a high-order polynomial approximation of the weak formulation of the wave equation. It combines the accuracy of the pseudospectral method with the flexibility of the finite element method. In this method, the wavefield is represented in terms of high-degree Lagrange interpolants per elements, and integrals are computed based upon the Gauss-Lobatto-Legendre (GLL) quadrature. Combined with a tensorial formulation, it leads to a perfectly diagonal mass matrix, which in turn leads to a fully explicit time scheme that is suitable to numerical simulations on parallel computers. The method allows a low dispersion error together with an accurate and implicit description of the boundary conditions. Nevertheless, the tensorial formulation imposes to use hexahedral meshes. We model cut and glued rock samples mimicking laboratory conditions as close as possible.

The Galerkin weak form of equation of motion (Geller & Ohminato, 1994;Komatitsch & Vilotte, 1998) reads:

            2 2 ( ) ( ), t t t T H u g (1)
where E T is mass matrix, the stiffness matrix E H relates the elasticity and anelasticity, E g represents the source term, and E u is the discretized displacement for each element/point in the model. We use the free surface boundary (natural boundary condition) for the whole medium. However, in the reality, the rock sample was surrounded by a silicone jacket with high attenuation and slow velocity that could play a role as the partial-absorbing boundary with respect to the rock sample.

In order to model anelastic attenuation in our time-marching spectral element methods, SEM3D uses standard Zener linear solids (SLS) to approximate a nearly constant quality factor E Q for P and S waves [START_REF] Carcione | Wave propagation simulation in a linear viscoacoustic medium[END_REF]Emmerich & Korn, 1987;Liu et al., 1976;[START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF]. A Zener body consists of a spring in series of the medium properties in connection with their relaxed stage at time   E t . SLS assumes that Q does not depend on frequency. In this study, we use three SLS on a band of relaxation frequency ranging from 33.3 kHz to 3.3 MHz.

Meshing

In order to simulate realistic wave propagation, for the sake of seeking boundary conditions including silicone rubber over the rock, one needs to carefully mesh the medium in SEM computation. The standard GLL quadrature of SEM requires hexahedral meshes. However, it is known to be challenging to create such nonstructured hexahedral meshes with acceptable computational cost and high accuracy for a complex geometry. In this study, we used Trelis software (https://urldefense.com/v3/__https://coreform. com__;!!N11eV2iwtfs!6Ub8yKpVonc2RehbF-2k5BPhPC0B8mjBgvP9GS2FRSoHsV0lCzzuiSmjvL04a2M$) to mesh the model. In the case of nonstructured meshes, one has to avoid distorted elements, which may cause numerical instabilities.

Figure 1c shows the meshing used in our forward modeling. To avoid the numerical error caused by the contact between more than one medium, we must include a layer as refinement at the boundaries of different materials. Since we use polynomials of degree five as trial functions inside each element, the wavelength should be larger than the element size. Considering the minimum S-wave velocity of 2,500 m/s in the granite rock and the maximum target frequency of 2 MHz, the mesh size is set up to be approximately 0.3 mm in all models.

Boundary Conditions

The rock samples used in laboratory experiments are roughly cylindrical rocks. However, the sample was cut to form two parallel sharp surfaces, hence forming an incomplete cylinder. This imperfection of the cylindrical nature can cause a non-negligible effect on the waveforms. Although the rock sample itself is the main target of our simulation, a single-material approximation with free surface boundary would fail LAI ET AL.

10.1029/2021JB021873

5 of 23 in waveform modeling due to the surrounding materials. Materials such as silicone rubber, glue, aluminum frame, shrink plastic attached to the granitic rock may affect the waves propagating in the rock. Therefore, we incorporate complex boundary conditions to find the suitable model setting without losing computational efficiency. Furthermore, we investigate the boundary conditions with additional physical layer of Shin-Etsu KE45 W RTV silicone rubber (Table 1) to simulate the transducer waveforms and obtain the preferred model (Figure 1c). The comparison of the four models tested for boundary condition is described in Appendix A. The additional boundary conditions successfully mimicked quasi-surface waves generated at the rock-rubber discontinuity.

Source and Transfer Functions

In addition to the boundary setting and model geometry, the source time function and the source-receiver relative position will also cause a significant impact on the waveforms. To mimic the laboratory experiment, we set P-wave source as an isotropic explosion source and S-wave source to be a shear stress along x  E z direction. The source time function is formally taken from the experiment input waveform and then filtered as described in Appendix B.

During the experiments, compressional and shear waves are transmitted independently at different source positions and their strain waveforms are recorded at their corresponding antipodal points. In the reality, the transducers are glued to the rock sample. However, strains are not always continuous and we have less knowledge on coupling between the transducers, silicone rubber, and the rock sample. We thus consider four possible effective source-receiver distributions: (a) both the source and the receiver are slightly shifted inside the region of the rock model; (b) both the source and the receiver are slightly shifted to the region of the surrounding model; (c) the source is slightly shifted to the region of the surrounding model and the receiver is inside the region of the rock model; and (d) the receiver is slightly shifted to the region of the surrounding model and the source is inside the region of the rock model. The four source-receiver position scenarios are tested against the observed data for    0 E (Figure 2). Compared to setting (b), the simulated waveforms of settings (a), (c) and (d) do not match the observed waveform. When the source and receiver are both located slightly outside the rock region (setting b), the first phases up to the fifth peak (4.0-6.2 μ E s) have higher amplitudes than the other configurations. Therefore, we use the setting (b) as the source-receiver configuration in this study. We will discuss in detail the possible reasons for this preference in Section 4.

Processing of Observed and Synthetic Data

In the study, we aim to directly compare the observed and synthetic waveforms in order to infer the structural changes during the deformation experiment. By carefully treating 3D geometry, attenuation and boundary conditions as we discussed above, our method is advantageous to foster understanding of rock deformation compared to our previous studies (e.g., Lai et al., 2019). We had to introduce transfer functions to compensate unrealistic numerical settings such as 2D assumptions, pure-elastic approximations and free-surface boundary conditions.

According to the time-dependent spectrogram analysis in Figure 3, the strongest energy of the experimental waveform is concentrated in the frequency band between 1.6 and 3.0 MHz for P-wave records and in the frequency band between 1.0 and 2.8 MHz for S-wave. Compared to S-wave temporal spectrogram that indicates a constant dominant frequency range, P-wave spectrogram shows lower frequency energy (  E 500 kHz) after 7 μ E s in addition to the most concentrated high-frequency region at the beginning. In this study, however, we focus on the variation in the high-frequency regime in order to capture the influence caused by inner structure changes. Two Butterworth bandpass filters with frequency ranges from 1.5 to 3.5 MHz and from 1.0 to 3.0 MHz are applied on the P-and S-wave data, respectively. A Hann window is then applied to data after filtering to extract the first arrivals. Note that we applied the same data processing treatments both to observed and synthetic waveforms.

Inverse Problem

We use waveforms for differential stress 

 E :            ( ) ( ) ( ) T T T PP SS d d d (2)
where PP E d and SS E d are time series vectors of compressional and shear strain waveforms recorded at the antipodal receiver (Figure 1a). The superscript E T denotes transpose. We then generate synthetic seismograms ( ) E u m with the same geometrical and source configurations for a set of seismological structural parameters E m . We define a objective function to maximize:

S T ( , ) ( ) ( ) | ( )|| ( )| ( ) ( ) | ( ) | | m d u m d u m d u m d u                   1 1 ( ( )| m 2 (3)
with  E are difficult to fix from one experiment to another. This should be discussed more in the upcoming studies but here in this study we use our ad-hoc values as discussed in results section.

As the Aji granite is fine-grained and nearly isotropic, we can assume a homogeneous medium with a wavelength of 2-3 mm. We use four parameters to represent the model vector  ( , , , )

T P S P S E V V Q Q m
. This parameterization is to homogenize 3D an-/elastic heterogeneities. The time-lapse change of these global parameters will be related to the microscopic short-wavelength structural changes such as crack generation. The anelastic attenuation can be an indicator of the 3D elastic heterogeneity and intrinsic attenuation.

Based on the a priori information on the Aji granite, we set P E V varying from 3,700 to After determining the preferred model for the initial state 0 E m (Figure 1c), we performed 6,000 simulations of different velocity models and different attenuation models. Since the volume of the rock changed little in the experiment, we did not consider the variation in the rock density caused by the deformation (the density is fixed to be 2,650 kg/ 3 m E

). We used constant elastic and anelastic parameters of the silicone rubber jacket whereas we systematically modify the parameters of the rock model for the purpose of creating a large database as shown in Table 2. We first simulate the P-wave velocity from 4,700 to 5,800 m/s with a step of 20 m/s, and from 3,700 to 4,600 m/s with a step of 100 m/s. We ranged S-wave velocity from 2,660 to 3,550 m/s with a step of 10 m/s (cf. Watanabe & Higuchi, 2015;Zaima & Katayama, 2018).

Based on the above considerations, in our first step of data processing, for P-wave waveform simulation, we vary the elastic parameters of P E V , P E Q , and the V V P S

/ while fixing S E Q to be constant. During the P-wave analysis, we vary V V P S / from 1.5 to 1.7 with step equals to 0.05 (cf. Zaima & Katayama, 2018). For the S-wave waveform simulation, we change the elastic parameters of S E V , S E Q , and the V V P S / while keeping P E Q constant. During data processing, we derive the best models using Equation 3, refining the range of parameters to reach to the best values. We use preferred models during the first numerical simulations to localized our inversions.

Results

Snapshot of the Wavefield

Figure 4 shows snapshots of the P-wave strain wavefield in our preferred model for 

  0 E from 2 μ E s to 7 μ E
s with a time step of 1 μ E s. These snapshots are cross sections of the wavefield at the same height as the source and receiver positions. At t = 2 μ E s, the snapshot shows the direct P-wave radiation from the source towards the antipode. No scattering or conversion at the sample surface occurs at this time. However, due to the side-cut surface, the wavefronts have no longer a regular concentric form. The surrounding low-velocity layer outside the sample partially absorbs energy and partially reflects it. This is why the reflection from the sample surface becomes too weak to interference with the direct P waves. Since our model is not a perfect cylindrical model, the two cut corners of the source side immediately provoke reflections. When the direct wave is reflected at the opposite side ( t = 4 μ E s), the reflected P phase (PP wave) is generated and propagates back to the source side. Since the receiver side has two cut corners as well, two additional reflected waves are generated and interfere with each other between 5 μ E s and 6 μ E s (Figures 4d and4e). However, the external low-velocity silicone layer traps some energy of the direct waves, and the reflections get smaller and smaller. The P-wave propagation snapshot at a later time ( t = 7 μ E s) shows a complex feature due to the contamination of reflected and converted waves developed at the side-cut surface.

Figure 5 presents the snapshots of S-wave wavefield. Most of side-cut geometry effects mentioned above can also be seen in S-wave snapshots. The S wave arrives at the receiver at about 6 μ E s. What differs the S-wave propagation from that of P waves is the strong influence from guided waves along the sample surface. The direct S wave (body wave) at the antipode receiver is contaminated by the waves coming subsequently along the interface (quasi-surface waves) (Figure 5e). This indicates the difficulty of extracting information on the interior of the sample using S waves alone. V is 940 m/s (∼ 18% E ). As for S E V , the difference is about 560 m/s (∼ 16% E ). The increase in velocity is larger for compressional waves (  E 100 m/s) than for shear waves (  E 40 m/s) (Figure 8a). Similarly, the velocity evolution of the wet sample (Figure 8b) follows the three stages of dry sample. In wet experiment, the difference between the maximum and the minimum P E V is 320 m/s (∼ 6.1% E ) and the S E V difference is about 425 m/s (∼ 13.7% E ). Figure 9 shows that increasing compression stress has larger impact on P E

V changes in dry data than in the wet data. Among all the data sets, the S E dV of wet data decreases most when having the same stress stage. Furthermore, the turning point of velocity changes of wet data shows earlier (  E 60 MPa) than dry data (  E 100-150 MPa). The increasing and transition stages of the wet data are much shorter than those of dry data.

Attenuation Parameters Evolution

The variations of 1 ) along the increase of differential stress. The E Q evolution obtained from wet data is particular: the attenuation decreases first and then increases, which could be coherent with the evolution of the velocity. Since, the water-saturated sample breaks earlier than

E Q to initial 1 E Q ( Q Q   1 0 1 /
) of wet and dry data.

Journal of Geophysical Research: Solid Earth

LAI ET AL.

10.1029/2021JB021873

14 of 23 in water-saturated rock. Under the same compressional pressure, attenuation of wet sample is higher than dry sample. In general, the wet sample is more sensitive to the attenuation variation.

Further Waveform Analysis

In this subsection, we observe the waveform discrepancies between the synthetics and observed data. The first observation is that the first arrivals of both the P-and S-wave waveforms gradually arrived earlier and later as the compressional pressure increases. In Figures 11a and12a, the first three peaks of the synthetic P-strain waveforms fit well with the experimental data. The synthetics do not fit well with the observed data after the fifth peak (i.e., 6 μ E s of Figures 11a and6.8 μ E s of Figure 12a) due to the incapability of accurate shear-wave modeling. When   E is 349.23 MPa in the case of the dry data set, the previously disappearing fifth peak reappeared in Figures 11a, which may be again due to the poorly constrained shear waves. In laboratory, the S waves seem to have less effects on the recorded P-wave strain waveforms. When   E exceeds 373.75 MPa in the case of the dry data, the interior structure of the rock is discontinuous and, hence, difficult to homogenize, what makes our modeling inefficient. In Figures 11b and12b, synthetic S-strain waveforms match well with the first two peaks experimental waveforms. After the third peak, the experimental waveforms have more peaks and smaller amplitudes, causing the synthetic waveforms to fail to match the experimental data. However, the matching of the first few peaks provided enough information on the velocity model and attenuation parameters change.

Discussion

In this study, we developed a methodology to invert time-lapse active seismic data in laboratory in order to infer the changes of the elastic and anelastic parameters. Here, we first discuss the numerical configuration to yield synthetics sufficiently comparable against the observed data. Thence, we resume the observation and make further interpretation on our experiments and show the need for 3D waveform inversions in near future.

Effect of the 3D Model Settings and Geometry

There have been few studies that simulated elastic wave propagation at centimeter scales and compared the synthetics against observed data. Yoshimitsu et al. ( 2016) is one of the few examples and they performed their experiments on a stationary aluminum sample while we did for deforming granites. In addition, they used lasers for both a source and a receiver while we had to use transducers due to the deformation experiments configuration. They observed the body-and surface-wave propagation being influenced by the cylindrical geometry. In our numerical modeling with free surface conditions, we could observe the same phenomenon of large-amplitude quasi-surface waves emerged along the curved surface of the sample after the body waves. However, in our experimental data, these large-amplitude surface waves were not present due to the surrounding glue layer. In order to reproduce their absence, we had to introduce special boundary conditions. While Yoshimitsu et al. (2016) added absorbing boundary conditions on the top and at the bottom, we introduced an additional layer with physical properties of silicone surrounding the target rock model, making waveforms much more realistic (Figure 1c). Furthermore, we found that the two cut sides parallel to compressional pressure are an indispensable condition to take into account. A model with these cut sides provoke an important amount of extra reflections compared to a simple cylindrical model. To conclude, it is important to consider the three-dimensional geometry.

The surrounding silicone rubber layer can provide a buffer-like region, which can slow down and reduce the amplitude of quasi-surface-wave propagation between the two media. Note that when the direct waves reach the interface between the two materials and travel through a low-E Q and low-velocity material, the resulting wavefield is amplified. The transmission coefficient calculated to confirm the observed phenomenon is shown in Appendix A.

Observed Changes in Wave Velocity and Attenuation

The evolution of homogenized velocity and attenuation in laboratory under a physical state should provide information on micro structural changes. The simultaneous Monte Carlo inversion for four homogenized elastic and viscoelastic parameters, P E V , S E V , P E Q and S E Q , were able to eliminate the cross-talks of these Kudo et al., 1992), the effective wavelength can be estimated as around 1.5-3.0 mm using the apparent dominant frequency (  E 1-2 MHz) and shear-wave velocity (  E 3,000 m/s). The local heterogeneity of grain size can therefore be neglected for the first place. However, we would still like to refine the structure using our waveform data but this simple estimation justifies the homogeneous modeling in this preliminary study.

The variation of wave velocity (Figure 8) during deformation estimated from our waveform matching methodology has a similar trend as seen by Zaima and Katayama (2018). Hereafter, we call the P E V and S E V measured by hand-picking "apparent" P E

V and S E V . Our results show slightly lower P-wave and higher S-wave velocities compared to the apparent velocities. However, apparent velocity models and the models obtained by our waveform-matching scheme do not correlate at large stress stage, especially for S-wave velocity. This indicates that the hand-picking fails even for first arrivals due to the ambiguity of travel times for highly dispersed waveforms, whereas waveform matching can still detect the phase arrivals objectively. Figure 9 shows the differential P E V and S E V obtained using waveform matching with respect to the initial status for dry and wet experiments. Since the dry sample resists to the external stress more than the wet sample, we can observe that the dry sample's velocities drop down to  E 15-18 % E while the wet sample shows 5% for P E V and 12.5% 

for S E V .

Another advantage of this Monte Carlo waveform inversion for elastic and anelastic parameters is the direct measurements of P E Q and S E Q . While Toksöz et al. (1979) studied the spectral ratios of waveform records of rock experiment relative to the waveform of reference sample, Lockner et al. (1977) and Zaima and Katayama (2018) measured the amplitude ratios of the first and second peaks in order to indirectly infer the attenuation, failing to obtain quality factors. Including P E Q and S E Q in the simultaneous simulations, objective functions show local maxima, especially for 1 E L -norm and 2 E L -norm. In particular, synthetic waveforms ignoring P E Q and S E Q are not comparable to the highly compressed experimental data. The reason of these differences may be due to seismic dispersion induced by attenuation in viscoelastic media, which certainly have more impact on shear waves than on compressional waves. Kjartansson (1979) modeled the relationship between velocity and the quality factor and obtained the increases of the velocity with frequency and 1 E Q . When looking at the discrepancy between S-wave velocities at the last stage of both dry and wet experiments obtained with hand-picking and waveform matching including attenuation, "apparent" values are constantly higher than our results. This is because of the waveform broadening due to low E Q at the last stage.

Since Zaima and Katayama (2018) measured only the amplitude ratios of first and second peaks to infer the attenuation, we cannot directly compare our results to theirs. Nevertheless, the amplitude of the S-wave is reduced more than that of the P-wave during deformation under saturated wet conditions, which is consistent with our results (Figure 10). Both P-and S-wave quality factor evolutions of wet data show roughly the same trend as the velocity evolution. The velocity increases while E Q increases, and E Q grows while the velocity grows.

Rock Physical Interpretation

The simplest way of interpreting the general trend of velocity changes is presented in many studies. At low confining pressure, the micro cracks in the experimental sample can be regarded as randomly distributed. The second stage is the closing process of horizontal micro cracks, which are perpendicular to the maximum stress, leading to velocity increase. The final stage, the vertical cracks form before the failure. However, as the closure of horizontal micro-cracks has little effect on the energy dissipation of P and S waves, it is difficult to see the difference during the attenuation process.

Experimental studies found seismic wave attenuation in dry rocks was sensitive to intrinsic material anelasticity, thermoelastic effects, friction at grain boundaries and micro-cracks development (e.g., Gordon & Davis, 1968;Guéguen & Schubnel, 2003;Johnston et al., 1979;Simmons & Brace, 1965;Tisato & Quintal, 2014;Walsh, 1965;Winkler & Nur, 1982). Water-saturated rocks, on the other hand, could accommodate fluid-related mechanism that would provoke energy dissipation, which should be responsible for seismic attenuation. Winkler and Nur (1982) declared that the effect of frictional sliding to the attenuation is negligible for pore-fluid rocks. Hence, they emphasized the importance of the role of intercrack fluid flow, which

Journal of Geophysical Research: Solid Earth

LAI ET AL.

10.1029/2021JB021873

17 of 23 may be considered as the major cause of the attenuation. In their experiments of dry, partially saturated, and saturated rock, they showed that S-wave attenuation increases with the degree of saturation and reaches its maximum at total saturation. However, as saturation continues to increase to more than 90 % E water saturation, P E Q is higher than during partial saturation condition and 

P S E Q Q . These observations can be explained with local fluid flow mechanisms. Guéguen and Schubnel (2003) also referred to the importance of local fluid flow and squirt flow to the attenuation and anisotropy. They numerically obtained dispersion for two different transversely isotropic distribution of cracks. When cracks are aligned vertically, and therefore orthogonally to the source-receiver horizontal plane, SH waves are more dispersive than P waves.

Figure 9 compares the velocity evolution of dry and fully saturated wet experiments. The velocity results show that the presence of saturated water accelerates the rock sample to reach its maximum velocity, which we think is the non-crack mode. The presence of water also causes the fast reduction of shear waves. Figure 10 shows the waveform fitting results of dry and wet experiments, showing that attenuation changes with the increasing compression stress under different water conditions. The wet S wave shows a greater decrease in both velocity and quality factor, compared with the P wave, from low pressure until the rupture. These changes in the elastic and anelastic behavior are consistent with the observation of Winkler and Nur (1982), who speculated these effects as intercracks flow. At the beginning, increasing of pressure causes little effect on both the P-and S-wave attenuation of dry and wet samples, because there is no enough density of cracks that affect the wave propagation. As the effective stress increases, our results show that  P S E Q Q in the wet case, whereas  P S E Q Q for the dry case. This may be produced when there is a compressional pressure high enough to form micro-cracks with higher density, especially with micro-cracks aligned vertically along the long axis. These highly dispersed waves lead to a large amount of attenuation. The characteristics and the local pore fluid flow can explain our attenuation results obtained from fully saturated sample (Guéguen & Schubnel, 2003;Winkler & Nur, 1982).

The contribution to attenuation in the dry rock is assumed to be only due to the friction and the intrinsic aggregate attenuation. Intrinsic attenuation is thought to be caused by energy dissipation due to friction at cracks where those faces are barely in contact, making amplitude changes sensitive to crack geometry in the specimen (Lockner et al., 1977;Walsh, 1966). Besides, Bonner (1974) presented an increase of shear wave anisotropy in the dry granite caused by cracks oriented parallel or oblique to the compressional stress. The cracks oriented parallel and subparallel to the axis of maximum compression stress have influence on attenuation. However, we can see that many laboratory experiments confirmed that  P S E Q Q exists in the dried samples (Johnston et al., 1979). No specific mechanism has been proposed to explain the difference between P-and S-wave attenuation in dry rock.

Figure 13 compares the exact and relative P E V / S E V ratio obtained from waveform matching with the P E V / S E V ratio of Zaima and Katayama (2018). A number of studies have mentioned that apparent P E V / S E V ratio rise can be found in the seismic data of the subduction zone or laboratory experimental data, which may be linked to the high pore fluid pressure and crack anisotropy when approaching failure (Audet et al., 2009;Christensen, 1984;Peacock et al., 2011;Wang et al., 2012). It is worth noting that there is a difference between the velocity variation curve we corrected through attenuation and the apparent V P S /V ratio (Figure 13b). Our results present higher P E V / S E V ratios for wet experiment when rock approaches rock failure. Even in dry rock experiment, P E V / S E V ratios are slightly higher compared to the apparent P E V / S E V ratios before rupture. These phenomena may be caused by the sudden rise of S-wave anisotropy (Bonner, 1974) and thus influencing the attenuation and wave velocity. This characteristic of V V P S

/ ratio rising marks the particularity of the rock right before failure, and may provide robust information on fracture prediction.

Toward 3D Waveform Inversion

Our method can process large amounts of experimental data quickly and provide the best matching values of P E V , S E V , P E Q , and S E Q . However, the synthetic seismograms cannot match the later phases regardless of parameter sets. For example, the fifth peak (  E 6 μ E s) of the synthetic waveforms in Figure 11a cannot match the data, especially when it coincides with the arrival time of S waves. The difficulty may come from model homogenization that should not be valid for scattered S waves. The second possibility is that our boundary condition is still deficient, causing the S-wave waveforms to have poor matching. Since S waves and Therefore, full waveform inverse problem must inevitably move towards the development of heterogeneous models. Our method is an isotropic approximation of the deformed rocks, such that so far we can only match the first few peaks in the waveform to get the elastic and viscoelastic parameters. Kovalyshen et al. (2018); Kovalyshen et al. (2020) estimated anisotropy by inverting elastic wave arrival times for characterizing the damage and failure processes in rocks. They also discussed the effective interpretation and correction of waveforms recorded by the finite-size transducers in laboratory experiments. When we start to extract more information on the weakening of rock materials from the waveform, we must use a three-dimensional heterogeneous model. Indeed, it is only through heterogeneous model that we could better simulate the structural changes and mechanisms inside the rock during deformation.

Conclusions

Through the matching of the simulated and observed data, we obtained the velocity and attenuation changes as the pressure rises up to rupture in laboratory. Unlike previous studies, the estimation of attenuation no longer requires amplitude analysis but can be directly conducted through our simultaneous inversion schemes, without trade-offs between the unknowns. The "apparent" velocity change could be biased by ignorance of attenuation and therefore our estimation of P E V and S E V is more trustable. Combining seismology and rock physics, the data measured can explain the fracture nucleation during varying stress conditions. Seismology can bring useful information on small-scale rock fractures. Variations of elastic and anelastic parameters revealed in this study show different trends from the usual rock measurements, these observation and analysis may thus give us the opportunity to understand the rock failure mechanism of not only laboratory experiments, but also the shallower part of the Earth's crust.

Appendix A: Boundary Condition Tests

As shown in Figure A1a and A1b, the free-surface boundary rock model cannot reproduce the observations, due to its complete reflections from boundaries. Hence, we attempt to modify the boundary conditions by considering the silicone rubber, as in the real experiments: a layer of glue and silicone rubber is attached to the rock sample on its cylindrical side to jacket the sample. This layer is used to isolate the oil from the rock sample and transducers. The physical parameters of this layer are calculated from the characteristics of silicone rubber used in the experiment, Shin-Etsu KE45 W RTV silicone rubber (Table 1). The layer is added both to the cylindrical and to the side-cut cylindrical rock model. The additional layer in model (d) with the characteristics of silicone rubber can provide a buffer-like region that slows down the waves and reduce the amplitude of quasi-surface waves propagating between the two different materials. Figure A2 shows the coefficients of P to P and S to S transmission (Lay & Wallace, 1995), calculated as amplitude and phase at the solid-solid interface from rock to the silicone rubber. When the

  Nous combinons d'abord les données observées pendant des expériences statiques et le données synthétiques pour calibrer les signaux enregistrés. Puis, nous déterminons les propriétés des ondes qui se propagent dans les échantillons. Ensuite, nous obtenons des informations sur l'hétérogénéité et l'anisotropie 3D des roches en comparant les formes d'onde observées et les synthétiques. Enfin, les méthodes d'inversion de formes d'onde et d'imagerie de migration pour les sources actives sont mises en oeuvre sur les données réelles afin de reconstruire les changements de la structure localisée. La méthode d'inversion que nous proposons a un potentiel plus élevé que les méthodes classiques de laboratoire qui utilisent uniquement les temps d'arrivée. Mots-clés modélisation numérique, propagation des ondes, mécanique des roches, expérience en laboratoire, nucléation sismique, supervision, méthode des éléments spectraux, inversion de forme d'onde, déformation des roches

  , et les noyaux utilisés dans la migration sont trop simplifiés en utilisant la théorie des rayons généralisés. Le présent travail vise à étendre l'idée, à la combiner avec la modélisation numérique 3D sur la propagation des ondes, et à traiter soigneusement les formes d'ondes entières pour extraire des informations temporelles sur les changements des paramètres élastiques et anélastiques (figure 1.4). Comme nous souhaitons introduire les techniques d'imagerie sismique de pointe dans le laboratoire de mécanique des roches, nous présentons ici un bref aperçu de la sismologie globale et d'exploration. Vingt ans seulement après le premier enregistrement d'un événement télésismique (Tokyo-Potsdam) en 1889 (von Rebeur-Paschwitz, 1889), Gutenberg et Geiger ont utilisé moins de 100 traveltimes pour

  et al. (2004) a comparé les données sismiques mesurées pour estimer et prédire la migration dans le temps du CO 2 après l'injection de CO 2 dans le réservoir souterrain. Riahi et al. (2013) a obtenu les vitesses de phase et les angles d'azimut arrière du traitement des ondes de Rayleigh et de Love par un algorithme de formation de faisceau à trois composantes et les a utilisés pour effectuer une analyse time-lapse d'une installation de stockage de gaz souterraine. Planès et al. (2016) a appliqué la technique d'interférométrie sismique pour surveiller les changements temporels dans les barrages en terre et les digues. La réduction de la vitesse des ondes de surface a été estimée dans une expérience à l'échelle du laboratoire et observée sur un remblai de terrain. Cotton et al. (2018) a développé une migration temporelle 4D prestack pour détecter les changements de vitesse obtenus après la corrélation croisée des levés de base et de surveillance. Cette méthode est appliquée à des données précises sur le pétrole lourd pour révéler la propagation temporelle de la vapeur dans le réservoir. L'inversion time-lapse est ensuite développée en inversant les données mesurées à différents moments pour obtenir les changements des paramètres du modèle. La soustraction des données sismiques et synthétiques donne l'information des variations structurelles dans le milieu. La technique d'imagerie par migration est utilisée dans cette étude sur des données expérimentales pour obtenir la variation structurelle pendant les expériences de déformation.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Schematic diagram of the role of micro cracks in fracture process: A fault plane will occur due to the connection of micro cracks into one continual plane, whereas the plane will generally form at around 30°t o the direction of principle stress. . . . . . . . . . . . . . . . . . . 1.3 Gradual process of fracture nucleation. (a) Stress-strain plot showing the four stages in a compression experiment. (b) The changes in volume and the axial shortening show the dilatancy phenomenon. The dashed line indicates the change point from volume decrease to increase. (c-f) the behaviour of cylindrical rock during the successive stages of a confined compression experiment and accompanying stress-strain plot. (c) Pre-deformation state, showing preexisting cracks. (d) Compression begins, and volume decreases due to crack closure. (e) Crack propagation and dilatancy happen. (f) Merging of cracks along the future shear fracture. (Van der Pluijm and Marshack, 2004) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1.4 The schematic schema of our ultimate time-lapse imaging strategy during rock deformation experiments. The grey rectangle represents the cross-section of the sample. t = 1∆t, 2∆t and 3∆t indicate the deformation stages that the higher ∆t represents that higher deformed condition. Our idea is to build a system that 1) process a huge database of numerical models and their corresponding synthetic waveforms; 2) match the experimental waveforms and synthetic waveforms and obtain the degree (%) of agreement between the experimental data and the simulated waveform; 3) obtain the most potential structural changes among the first few best-matched models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 The time-lapse Z-component wavefield of displacement snapshots. The orange cross denotes the source and the green square denotes the receiver. a) t= 1.4; b) 3.4; c) 5.4; d) 7.4; e) 9.4; f) 12.4 µs. . . . . 2.2 Numerical modelled P-wave velocities as functions of the differential stresses of dry sample using hand-picking and auto-fitting. . . . 2.3 Observed waveform data plots in red dots and the synthetics with response in black with the new velocities in this study, cutting from 2 µs to 7 µs without S-wave arrivals. . . . . . . . . . . . . . . . . . 2.4 2D mesh used in ABAQUS software. The circular model is designed with 2 cm in diameter. . . . . . . . . . . . . . . . . . . . . . . . . 2.5 2D wave propagation in Abaqus. . . . . . . . . . . . . . . . . . . . 2.6 3D wave propagation in cylindrical model. . . . . . . . . . . . . . . 2.7 Wave propagation in the cross-section view along the source and receiver points of the 3D model. . . . . . . . . . . . . . . . . . . . . 2.8 The internal friction plots Q -1 . Three SLSs are used for approximate the attenuation behavior. . . . . . . . . . . . . . . . . . . . . . . . 3.1 The schematic diagram of the geometry of sample, the position of the transducers and the numerical model used in this study. (a) The general view of the experimental setting of the sample. (b) The crosssection of the sample. . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The schematic diagram of the triaxial compression test in Zaima and Katayama (2018). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The distorted hexahedral mesh of the model for SEM. The green region is the rock model. The yellow part represents the surrounding low-velocity material. The mesh size is approximately 0.3 mm. . . . 3.4 The four models tested in this study and their synthetic waveform (a) cylindrical model with free-boundary, (b) side-cut cylindrical model with free-boundary, (c) cylindrical model with surrounding silicone jacketing model and (d) side-cut cylindrical model with surrounding silicone jacketing model. . . . . . . . . . . . . . . . . . . . . . . . 3.5 The transmission coefficient P-to-P (P-SV) and S-to-S (SH) of solidsolid surface between the granite and the silicone rubber plotted with incidence angle. The frequency range is set to the dominant frequency 2MHz. (a) P-to-P (P-SV) amplitude, (b) P-to-P (P-SV) phase, (c) S-to-S (SH) amplitude, (d) S-to-S (SH) phase. . . . . . . 3.6 The input wavelet of the rock experiment (top) and the numerical modelling (bottom). In order to avoid the noise of the experimental wavelet, we use the filtered wavelet as the input. . . . . . . . . . . . 3.7 The source-receiver positions tested on the preferred model (figure 3.3) and their normalised-synthetic waveform compared to the data. The red circle represents the source position; the green triangle indicates the receiver position. The four source-receiver positions: (a) source and receiver are both slightly inside the inner material, (b) source and receiver are both slightly outside the inner material, (c) source is slightly outside, but the receiver is slightly inside the inner material, (d) source is slightly inside, but the receiver is slightly outside the material. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 (a) Waveforms and (b) FFT spectrum of observed P-wave data of the dry experiments. The data of each stress state are plotted with different colors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 (a) Waveforms and (b) FFT spectrum of observed S-wave data of the dry experiments. The data of each stress state are plotted with different colors. . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.10 The (a) P-and (b) S-wave spectrogram of experimental data for 0 MPa differential stress. . . . . . . . . . . . . . . . . . . . . . . . 3.11 Spectral ratio of observed P-wave waveforms of the dry experiments.

3. 18

 18 The results of waveform fitting using different cost functions for (ac) P-and (d-f) S-wave, respectively. (a,d) Zero-lag cross correlation, (b,e) L1-norm and (c,f) L2-norm. . . . . . . . . . . . . . . . . . . . 3.19 The result of waveform fitting using objective function (equation 3.2). The possible trend of change in (a) P-wave and (b) S-wave velocity with increasing differential stress is shown as the high extreme of the scale (dark blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Evolution of V P and V S of (a) dry and (b) wet data. The black and red circles represent the numerical simulated best-fit results of V P and V S , respectively. The diamonds are the results from Zaima and Katayama (2018), which are calculated from the hand-picked first arrivals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Ratios of differential velocities to initial velocity (dV /V 0 ), showing the velocity change compared to the original stage at each stress condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.22 (a) Q -1 of wet and dry data, (b) Ratios of differential Q -1 to initial

  wet and dry data . . . . . . . . . . . . . . . . . 3.23 The result of (a) P-and (b) S-wave strain waveform fitting of dry data using equation 3.2 to find the best-match model with proper V P , V S , Q P , Q S . The black lines are the best-match synthetic waveforms. The red lines are the experimental data. The y-axis indicate the deformation progress and correspond to a condition of ∆σ . The x-axis represents the time of wave propagation. . . . . . . . . . . . . . . . 3.24 The result of (a) P-and (b) S-wave strain waveform fitting of wet data using equation 3.2 to find the best-match model with proper V P , V S , Q P , Q S . The black lines are the best-match synthetic waveforms. The red lines are the experimental data. The y-axis indicate the deformation progress and correspond to a condition of ∆σ . The x-axis represents the time of wave propagation. . . . . . . . . . . . . . . . 3.25 (a)V P /V S ratio change, (b)Relative V P /V S ratio change. . . . . . . . . 4.1 Schematic diagram of the experiment setting: (a) Transducer-laser test and (b) Transducer-transducer test. Source transducers are placed at the midpoint of the long side of the sample, while the laser and the receiver transducer are placed at the opposite side of the source. . . 4.2 Using the laser Doppler interferometer as the receiver on an aluminium sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Metal samples used for the tests (from lest to right): copper, hollow aluminium, short aluminium-1, short aluminium-2, brass, and long aluminium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Twenty-four sensors with numbers (T01-T24). . . . . . . . . . . . . 4.5 Twenty-four points (P01-P24) marked on the surface of the aluminium sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 T10 (P9) sensor as the source and the laser (P13) as receiver. . . . . 4.7 T10 (P09) sensor as the source and T15 (P13) sensor as receiver. . . 4.8 T15 (P09) sensor as the source and T10 (P13) sensor as receiver. . . 4.9 T10 (P13) sensor as the source and T15 (P09) sensor as receiver. . . 4.10 Correlation tests for the waveforms (figure 4.6 and 4.7) using sensor T10 as the source and laser as the receiver (green) and sensor T10 as the source and sensor T15 as the receiver (blue) with (a) 200 mV and (b) 4 V as FSV and different damping factors 1, 3 and 5. . . . . 4.11 Correlation tests for the waveforms (figure 4.8 and 4.9) using sensor T15 as the source and sensor T10 as the receiver (magenta) and sensor T10 as the source and sensor T15 as the receiver (cyan) with (a) 200 mV and (b) 4 V as FSV and different damping factors 1, 3 and 5.4.12 Sensor consistency tests of 200 mV FSV, using sensor T10 as the source and different six sensors T15, T16, T17, T18, T19 and T20 as the receiver. The source sensor was always placed at the point P09, whereas the receiver sensor was placed at the point P13. . . . . . . . 4.13 Sensor consistency tests of 4 V FSV, using the same experiment geometry as figure 4.12. . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Reciprocity tests using a set of waveforms using (a) sensor T15 as the source and sensor T10 as the receiver; (b) sensor T10 as the source and sensor T15 as the receiver. . . . . . . . . . . . . . . . . . . . . 4.15 Sensor consistency tests with six sensors: T15, T16, T17, T18, T19 and T20. The reference trace for the cross correlation was the waveform recorded using the T10 as the source and T15 as the receiver. It is used for cross correlate with all the waveforms recorded recorded by the six sensors (T15 to T20). . . . . . . . . . . . . . . . . . . . 4.16 Source time functions of transducer-laser tests obtained by deconvolution. (a)-(d) are the results using different windows for deconvolution. The green line indicates the theoretical first arrival time. . . . . 4.17 Source time functions of transducer-transducer tests obtained by deconvolution. (a)-(d) are the results using different windows for deconvolution. The green line indicates the theoretical first arrival time. 4.18 Synthetic waveforms of transducer-laser tests using the source time functions in figure 4.16(a)-(d). . . . . . . . . . . . . . . . . . . . . 4.19 Synthetic waveforms of transducer-transducer tests using the source time functions in figure 4.17(a)-(d). . . . . . . . . . . . . . . . . . 4.20 Raw recorded sync signal of the source input signal (red) and its integral (blue) and differential (green) waveform. . . . . . . . . . . 4.21 Four direct laser measurements at the surface of sensors using different settings of full scale Volts(±) and damping: 200mV-damping 3, 4V-damping 3, 4V-damping 1 and 4V-damping 5. . . . . . . . . . . 4.22 Source time function of (a) transducer-laser test and (b) transducertransducer test obtained from implementing deconvolution in frequencydomain. The green line indicates the theoretical first arrival time. . . 4.23 Source wavelet of transducer-laser test and its FFT spectrum. . . . . 4.24 Source wavelet of transducer-transducer test and its FFT spectrum. .

4. 25

 25 figure 6.7. (a) The entire cylindrical model with a small perturbation (0.1%) in the middle (green). (b) Zoom in of the patch of small perturbation. The volume ∆V is 1.64×10 -9 m 3 . . . . . . . . . . . . 6.10 Comparison between the real displacement difference and the difference obtained from partial derivative in our perturbation model B (figure 6.9). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.11 Wavefield difference related to the middle perturbation point of model A (figure 6.7) by time steps of every 4 µs. . . . . . . . . . . . . . . . . 6.12 Time-lapse P-wave kernels related to the middle perturbation point in model A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Fault zone structures[START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF]: (a) a model with a single high-strain core surrounded by a fractured damage zone[START_REF] Chester | Implications for mechanical properties of brittle faults from observations of the Punchbowl fault zone, California[END_REF] and (b) a multiple-core model, where many strands of high-strain material enclose lenses of fractured protolith[START_REF] Faulkner | On the internal structure and mechanics of large strike-slip fault zones: Field observations of the Carboneras fault in southeastern Spain[END_REF] .

Figure 1 . 3 :

 13 Figure 1.3: Gradual process of fracture nucleation. (a) Stress-strain plot showing the four stages in a compression experiment. (b) The changes in volume and the axial shortening show the dilatancy phenomenon. The dashed line indicates the change point from volume decrease to increase. (c-f) the behaviour of cylindrical rock during the successive stages of a confined compression experiment and accompanying stress-strain plot. (c) Pre-deformation state, showing preexisting cracks. (d) Compression begins, and volume decreases due to crack closure. (e) Crack propagation and dilatancy happen. (f) Merging of cracks along the future shear fracture. (Van der Pluijm and Marshack, 2004)

Figure 1 . 4 :

 14 Figure 1.4: The schematic schema of our ultimate time-lapse imaging strategy during rock deformation experiments. The grey rectangle represents the cross-section of the sample. t = 1∆t, 2∆t and 3∆t indicate the deformation stages that the higher ∆t represents that higher deformed condition. Our idea is to build a system that 1) process a huge database of numerical models and their corresponding synthetic waveforms; 2) match the experimental waveforms and synthetic waveforms and obtain the degree (%) of agreement between the experimental data and the simulated waveform; 3) obtain the most potential structural changes among the first few bestmatched models.

  ) Time-lapse elastic properties of cracked granite during deformation inferred from laboratory experiments using 2D waveform modelling, Proceedings of the 81st EAGE conference and exhibition, London. https://doi.org/10.3997/2214-4609.201900662 • Lai, S.-T., Fuji, N., Katayama, I., Bonilla, L. F., & Capdeville, Y. (2021) Rock deformation monitoring using Monte Carlo waveform inversion. Journal of Geophysical Research: Solid Earth, 126, e2021JB021873. doi:10.1029/2021JB021873 • Lai, S.-T., Fuji, N., Schubnel, A., Aubry, J., Capdeville, Y., Bonilla, L. F., & Meschede, M., Do Born-approximated waveforms predict grain size and anisotropy of rock samples? (in prep.) • Fuji, N., Lai, S., Yuan, S., Katayama, I. (2018). One-station time-lapse seismic imaging: concept and preliminary applications, AGU Fall Meeting. • Lai, S.-T., Fuji, N., Bonilla, L. F., Schubnel, A., Capdeville, Y., & Katayama, I. (2021) Time-lapse rock deformation monitoring using waveform inversion, AGU Fall Meeting.

  to reconstruct synthetics. The parameters used are ∆t=0.002 µs, ∆x = ∆z =0.005 cm. The source wavelet is set as the real input source used in the experiments with 2 MHz resonant frequency in compressional and shear axes. Since the source is on the long side of the cylinder, we put the point source on the top and put one receiver at the bottom of the 2D circle (figure 2.1). Here we can use the reciprocity theorem to model synthetic waveforms for the sources of two kinds: compressional and shear, just by using a single force from the receiver to the antipode source position and calculating the compressional and shear strains (cf. equation 2.9). We start with a homogeneous model since the Aji granite is nearly an isotropic texture with a small compressional wave anisotropy (<1.3%), and shear wave velocity has no significant changes in different polarisation directions. Homogeneous velocity models are tested in different P-wave and S-wave velocities to find a best-fit model for the experimental waveform results during the deformation.

Figure 2 . 1 :

 21 Figure 2.1: The time-lapse Z-component wavefield of displacement snapshots. The orange cross denotes the source and the green square denotes the receiver. a) t= 1.4; b) 3.4; c) 5.4; d) 7.4; e) 9.4; f) 12.4 µs.

automated

  Figure 2.2: Numerical modelled P-wave velocities as functions of the differential stresses of dry sample using hand-picking and auto-fitting.

Figure 2 . 3 :

 23 Figure 2.3: Observed waveform data plots in red dots and the synthetics with response in black with the new velocities in this study, cutting from 2 µs to 7 µs without S-wave arrivals.

Figure 2 . 4 :

 24 Figure 2.4: 2D mesh used in ABAQUS software. The circular model is designed with 2 cm in diameter.

Figure 2 . 5 :

 25 Figure 2.5: 2D wave propagation in Abaqus.

Figure 2 . 6 :

 26 Figure 2.6: 3D wave propagation in cylindrical model.

Figure 2 . 7 :

 27 Figure 2.7: Wave propagation in the cross-section view along the source and receiver points of the 3D model.

Figure 2 . 8 :

 28 Figure 2.8: The internal friction plots Q -1 . Three SLSs are used for approximate the attenuation behavior.

Figure 3 . 1 :

 31 Figure 3.1: The schematic diagram of the geometry of sample, the position of the transducers and the numerical model used in this study. (a) The general view of the experimental setting of the sample. (b) The cross-section of the sample.

Figure 3 . 3 :

 33 Figure 3.3: The distorted hexahedral mesh of the model for SEM. The green region is the rock model. The yellow part represents the surrounding low-velocity material. The mesh size is approximately 0.3 mm.

Figure 3 . 4 shows

 34 the waveforms for (a) the cylindrical rock model with free-surface boundary, without surface cut; (b) a cylindrical model with free-surface boundary, with surface cut; (c) a cylindrical model with silicone layer without surface cut; and (d) a cylindrical model with silicone layer with the surface cut. The additional low-velocity layer surrounding the cylindrical rock sample behaves as ad-hoc modified boundary conditions to absorb wave energy while preserving some reflections. We can observe the subtle differences between models (a) and (b) in figure 3.4 caused by the cut sides. However, models (a) and (b) cannot reasonably reproduce the observed waveforms. Models (c) and (d) have similarity until 5.5 µs (the third peak), after which the reflections from the antipodal point are larger than the real data with the model (c). Due to the additional reflections back to the inner rock model in (d), the later phases show faster amplitude decay, and the synthetics for (d) resemble better to the observed datum than (c). The comparison among the four models demonstrates that

Figure 3 . 4 :

 34 Figure 3.4: The four models tested in this study and their synthetic waveform (a) cylindrical model with free-boundary, (b) side-cut cylindrical model with freeboundary, (c) cylindrical model with surrounding silicone jacketing model and (d) side-cut cylindrical model with surrounding silicone jacketing model.

Figure 3 . 5 :

 35 Figure 3.5: The transmission coefficient P-to-P (P-SV) and S-to-S (SH) of solidsolid surface between the granite and the silicone rubber plotted with incidence angle. The frequency range is set to the dominant frequency 2MHz. (a) P-to-P (P-SV) amplitude, (b) P-to-P (P-SV) phase, (c) S-to-S (SH) amplitude, (d) S-to-S (SH) phase.

  .7): (a) both the source and the receiver are slightly shifted inside the region of the rock model; (b) both the source and the receiver are slightly shifted into the region of the surrounding model; (c) the source is slightly shifted into the region of the surrounding model and the receiver is inside the rock model, and (d) the receiver is slightly shifted into the region of the surrounding model and the source is inside the rock model. The four source-receiver position scenarios are tested against the observed data for ∆σ = 0 (figure 3.7). Compared to setting (b), the simulated waveforms of settings (a), (c), and (d) do not match the observed waveform. When the source and receiver are located slightly outside the rock region (setting b)

Figure 3 . 7 :

 37 Figure 3.7: The source-receiver positions tested on the preferred model (figure 3.3) and their normalised-synthetic waveform compared to the data. The red circle represents the source position; the green triangle indicates the receiver position. The four source-receiver positions: (a) source and receiver are both slightly inside the inner material, (b) source and receiver are both slightly outside the inner material, (c) source is slightly outside, but the receiver is slightly inside the inner material, (d) source is slightly inside, but the receiver is slightly outside the material.

Figure 3 .

 3 Figure 3.8: (a) Waveforms and (b) FFT spectrum of observed P-wave data of the dry experiments. The data of each stress state are plotted with different colors.

Figure 3 .

 3 Figure 3.9: (a) Waveforms and (b) FFT spectrum of observed S-wave data of the dry experiments. The data of each stress state are plotted with different colors.

Figure 3 .

 3 Figure 3.10: The (a) P-and (b) S-wave spectrogram of experimental data for 0 MPa differential stress.

Figure 3 .

 3 Figure 3.11: Spectral ratio of observed P-wave waveforms of the dry experiments. The reference is the waveform at 0 MPa compression stress state.

Figure 3 .

 3 Figure 3.12: Spectral ratio of observed S-wave waveforms of the dry experiments. The reference is the waveform at 0 MPa compression stress state.

Figure 3 .

 3 Figure 3.13: Cross-correlation tests for P-wave waveforms of the dry experiments.

Figure 3 .

 3 Figure 3.14: Cross-correlation tests for S-wave waveforms of the dry experiments.

Figure 3 .

 3 Figure 3.15: Cross-correlation tests time-lag matrix.

3. 4

 4 Figure 3.16: The snapshot of the P-wave wavefield at different time steps from (a) 2 to (f) 7 µs, showing P-wave propagation process inside the medium.

Figure

  Figure 3.16 shows snapshots of the P-wave strain wavefield in our preferred model

Figure 3 .

 3 Figure 3.17: The snapshot of the S-wave wavefield at different time steps from (a) 2 to (f) 7 µs, showing S-wave propagation process inside the medium.

Figure 3 .

 3 Figure 3.18 shows the waveform fitting results on the dry sample using different cost functions: (a) zero-lag cross-correlation, (b) L 1 -norm and (c) L 2 -norm. Note that equation 3.2 uses reciprocals of objective (b) and (c) when we combine three objective functions. Regardless of the objective function type, blue indicates better waveform fits and red indicates disagreements between observed and synthetics.

Figure 3 .

 3 Figure 3.18: The results of waveform fitting using different cost functions for (a-c) P-and (d-f) S-wave, respectively. (a,d) Zero-lag cross correlation, (b,e) L1-norm and (c,f) L2-norm.

Figure 3 .

 3 Figure 3.19: The result of waveform fitting using objective function (equation 3.2). The possible trend of change in (a) P-wave and (b) S-wave velocity with increasing differential stress is shown as the high extreme of the scale (dark blue).

Figure 3 .

 3 Figure 3.20: Evolution of V P and V S of (a) dry and (b) wet data. The black and red circles represent the numerical simulated best-fit results of V P and V S , respectively. The diamonds are the results from Zaima and Katayama (2018), which are calculated from the hand-picked first arrivals.

Figure 3 .

 3 Figure 3.21: Ratios of differential velocities to initial velocity (dV /V 0 ), showing the velocity change compared to the original stage at each stress condition.

Figure 3 .

 3 Figure 3.22: (a) Q -1 of wet and dry data, (b) Ratios ofdifferential Q -1 to initial Q -1 (Q -1 /Q -10 ) of wet and dry data

Figure 3 .

 3 Figure 3.23: The result of (a) P-and (b) S-wave strain waveform fitting of dry data using equation 3.2 to find the best-match model with proper V P , V S , Q P , Q S . The black lines are the best-match synthetic waveforms. The red lines are the experimental data. The y-axis indicate the deformation progress and correspond to a condition of ∆σ . The x-axis represents the time of wave propagation.

Figure 3 .

 3 Figure 3.24: The result of (a) P-and (b) S-wave strain waveform fitting of wet data using equation 3.2 to find the best-match model with proper V P , V S , Q P , Q S . The black lines are the best-match synthetic waveforms. The red lines are the experimental data. The y-axis indicate the deformation progress and correspond to a condition of ∆σ . The x-axis represents the time of wave propagation.

Figure 3 .

 3 Figure 3.25: (a)V P /V S ratio change, (b)Relative V P /V S ratio change.

3. 5 . 1

 51 Effect of the 3D model settings and geometryFew studies simulated elastic wave propagation at centimetre scales and compared the synthetics against observed data. Yoshimitsu et al. (2016a) is one of the few examples, and they performed their experiments on a stationary aluminium sample while we dealt with for deforming granites. In addition, they used lasers for both a source and a receiver while we had to use transducers due to the deformation experiments configuration. They observed that cylindrical geometry influenced the bodyand surface-wave propagation. In our numerical modelling with free surface conditions, we could observe the same phenomenon of large-amplitude quasi-surface waves that emerged along the curved surface of the sample after the body waves.However, these large-amplitude surface waves were not present in our experimental data due to the surrounding glue layer. To reproduce their absence, we had to introduce special boundary conditions. WhileYoshimitsu et al. (2016a) added absorbing boundary conditions on the top and at the bottom, we introduced an additional layer surrounding the target rock model with physical properties of silicone glue, making waveforms much more realistic (figure3.3). Furthermore, we found that the two cut sides parallel to the compressional pressure are indispensable conditions to consider.

  also referred to the importance of local fluid flow and squirt flow on the attenuation and anisotropy. They numerically obtained dispersion for two different transversely isotropic distributions of cracks. When cracks are aligned vertically, and therefore orthogonally to the source-receiver horizontal plane, SH waves are more dispersive than P waves.

Figure 3 .

 3 Figure 3.21 compares the velocity evolution of dry and fully saturated wet experiments. The velocity results show that saturating water accelerates the rock sample to reach its maximum velocity. The presence of water also causes the fast reduction of shear wave velocity. Figure3.22 shows the waveform fitting results of dry and wet experiments, showing that attenuation changes with the increasing compression stress under different water conditions. The wet S wave shows a more significant decrease in both velocity and quality factors than the P wave, from low pressure until the rupture. These changes in the elastic and anelastic behaviour are consistent

  active-source seismic laboratory experiments are performed to 1) calibrate and verify the units of the amplitude of observed waveforms; 2) obtain the proper source time function from the experiments; and 3) understand the waveforms and the early coda propagating in metals and rocks. Two experiments were employed in this study to calibrate and investigate the exact components of the waveforms: transducer-tolaser and transducer-to-transducer tests. The laser can measure the absolute displacement (in one direction) without transfer function, and thus it is useful for calibration.Both of the experiments were implemented on different metal samples. All the experiments were conducted in Laboratoire de Géologie, École Normale Supérieure, with a static state at room temperature without applying external stress and conditions.4.2.1 Active-source one-station experimentThe schematic diagram of the experimental setting is shown in Fig 4.1. Two experiments: a) transducer-laser test and b) transducer-transducer test were implemented on all the samples to calibrate the parameters, obtain source time function through deconvolution and implement error-bar tests. The calibration experiments are mainly implemented on a copper sample with a diameter of 4 cm and a length of 10 cm and an aluminium sample with a diameter of 4 cm and a height of 8 cm. The laser Doppler interferometer used in this study is a Vibrometer OFV-5000 Modular Vibrometer Controller, combined with a Vibrometer Sensor Head OFV-503 from Politec Ltd (figure 4.2). The ultrasonic P-wave piezoelectric transducers were mounted at the midpoint of the long side of the cylindrical samples (figure 4.3).The elastic waves propagating through the samples were excited by the electric signal transmitted to the transducers. The vibration perpendicular to the surface was stacked 100 times and recorded by the laser or another transducer at the antipodal position.

Figure 4 . 1 :

 41 Figure 4.1: Schematic diagram of the experiment setting: (a) Transducer-laser test and (b) Transducer-transducer test. Source transducers are placed at the midpoint of the long side of the sample, while the laser and the receiver transducer are placed at the opposite side of the source.

Figure 4 . 2 :

 42 Figure 4.2: Using the laser Doppler interferometer as the receiver on an aluminium sample.

Figure 4 . 3 :

 43 Figure 4.3: Metal samples used for the tests (from lest to right): copper, hollow aluminium, short aluminium-1, short aluminium-2, brass, and long aluminium.

Figure 4 . 4 :

 44 Figure 4.4: Twenty-four sensors with numbers (T01-T24).

Figure 4 . 5 :

 45 Figure 4.5: Twenty-four points (P01-P24) marked on the surface of the aluminium sample.

Figure 4 . 6 :

 46 Figure 4.6: T10 (P9) sensor as the source and the laser (P13) as receiver.

Figure 4 . 7 :

 47 Figure 4.7: T10 (P09) sensor as the source and T15 (P13) sensor as receiver.

Figure 4

 4 Figure 4.8: T15 (P09) sensor as the source and T10 (P13) sensor as receiver.

Figure 4

 4 Figure 4.9: T10 (P13) sensor as the source and T15 (P09) sensor as receiver.

Figures 4 .

 4 Figures 4.6, 4.7, 4.8 and 4.9 show the waveforms used for correlation tests and reciprocity tests 1 and 2, where (a)-(f) are the records using: (a) 200 mV-damping 3, (b) 200 mV-damping 1, (c) 200 mV-damping 5, (d) 4 V-damping 3, (e) 4 V-damping 1 and (f) 4 V-damping 5 as FSV and source input damping.

Figure 4 .

 4 Figure 4.10: Correlation tests for the waveforms (figure 4.6 and 4.7) using sensor T10 as the source and laser as the receiver (green) and sensor T10 as the source and sensor T15 as the receiver (blue) with (a) 200 mV and (b) 4 V as FSV and different damping factors 1, 3 and 5.

Figure 4 .

 4 Figure 4.11 compares the results of the correlation tests using sensor T15 as the source and sensor T10 as the receiver (magenta) with sensor T10 as the source and sensor T15 as the receiver (cyan). This waveform record stability test is the preparation for the following reciprocity tests. Both tests show a high correlation among the

Figure 4 .

 4 Figure 4.11: Correlation tests for the waveforms (figure 4.8 and 4.9) using sensor T15 as the source and sensor T10 as the receiver (magenta) and sensor T10 as the source and sensor T15 as the receiver (cyan) with (a) 200 mV and (b) 4 V as FSV and different damping factors 1, 3 and 5.

Figure 4 .

 4 Figure 4.12: Sensor consistency tests of 200 mV FSV, using sensor T10 as the source and different six sensors T15, T16, T17, T18, T19 and T20 as the receiver. The source sensor was always placed at the point P09, whereas the receiver sensor was placed at the point P13.

Figure 4 .

 4 Figure 4.13: Sensor consistency tests of 4 V FSV, using the same experiment geometry as figure 4.12.

Figure 4 .

 4 Figure 4.14: Reciprocity tests using a set of waveforms using (a) sensor T15 as the source and sensor T10 as the receiver; (b) sensor T10 as the source and sensor T15 as the receiver.

Figure 4 .

 4 Figure 4.15 presents the sensor consistency tests with six sensors: T15, T16, T17, T18, T19 and T20. The reference trace for the cross-correlation was the waveform recorded using the T10 as the source and T15 as the receiver. It is used to crosscorrelate with all the waveforms recorded by the six sensors. Therefore, in every plot in figure 4.15, the column with high consistency represents the correlation value of reference trace with waveforms from the tests using T10 as the source sensor and T15 as the receiver. In the plots with the settings 200 mV-damping 5 and 4 V, we can see the results of six separate columns representing the tests of six different sensors. These development characteristics demonstrate that the sensors are the main controlling factor of the waveforms. Removing the specific effects of the sensor itself may be a necessary step when doing sophisticated analysis or using waveforms for imaging. Since using 200 mV-damping 3 and 4 V-damping 1 have the relatively highest values in sensor consistency tests in figure 4.15 and parameter setting of 200 mV-damping 3 had the better results in the previous tests, we chose 200 mV as FSV and the damping 3 (50 Ω) as the setting in the following experiments.

Figure 4 .

 4 Figure 4.15: Sensor consistency tests with six sensors: T15, T16, T17, T18, T19 and T20. The reference trace for the cross correlation was the waveform recorded using the T10 as the source and T15 as the receiver. It is used for cross correlate with all the waveforms recorded recorded by the six sensors (T15 to T20).

Figures 4 .

 4 Figures 4.16 and 4.17 are the source time functions obtained from frequency-domain deconvolution (chapter 4.3.1) using four different time windows. Since we know that the waveform recorded in transducer-laser experiments is velocity waveform, we numerical modelled the synthetic velocity waveform propagating in the copper sample shown in figure 4.18. After comparing the waveforms of displacement, velocity and acceleration, we concluded that the waveforms recorded in transducer-transducer experiments are also velocity waveforms (figure 4.19).

Figure 4 .

 4 Figure 4.16: Source time functions of transducer-laser tests obtained by deconvolution. (a)-(d) are the results using different windows for deconvolution. The green line indicates the theoretical first arrival time.

Figure 4 .

 4 Figure 4.17: Source time functions of transducer-transducer tests obtained by deconvolution. (a)-(d) are the results using different windows for deconvolution. The green line indicates the theoretical first arrival time.

Figure 4 .

 4 Figure 4.18: Synthetic waveforms of transducer-laser tests using the source time functions in figure 4.16(a)-(d).

Figure 4 .

 4 Figure 4.19: Synthetic waveforms of transducer-transducer tests using the source time functions in figure 4.17(a)-(d).

Figure 4 .

 4 Figure 4.20: Raw recorded sync signal of the source input signal (red) and its integral (blue) and differential (green) waveform.

Figure 4 .

 4 Figure 4.21: Four direct laser measurements at the surface of sensors using different settings of full scale Volts(±) and damping: 200mV-damping 3, 4V-damping 3, 4Vdamping 1 and 4V-damping 5.

Figure 4 .

 4 Figure 4.21 shows the direct laser measurement of the impulse sent by the vibrometer, passing through the piezoelectric transducer and sensor, using 200 mV and 4 V as full-scale Volts for recording leads to output waveforms with different features. The ones using 4 V but different damping are similar in waveform phases but differ in amplitudes. According to figures 4.20 and 4.21, our real source time function is more similar to setting (d) in figures 4.16 and 4.17. Therefore, the best source time functions of transducer-laser and transducer-transducer tests used for our numerical modelling were determined as velocity waveforms and are shown in figure 4.22, whereas the Fourier transform spectrum is shown in figures 4.23 and 4.24, respectively. The source time function obtained from the transducer-laser experiment has a higher dominant frequency of 1 MHz, whereas the transducer-transducer experiment has a dominant frequency around 2.0×10 5 Hz.

Figure 4 .

 4 Figure 4.25 indicates the real source time function (S') obtained from the entire deconvolution process in chapter 4.3.1. R and S in figure 4.25 represent the transducer effect, and the transducer effect combined with source time function (S'). The real source time function sent from the oscillator is 1 MHz. It shows that the transducer can be considered as a low-pass filter. We can also intuitively observe the transducer effects from figure 4.25(a) to (c).

Figure 4 .

 4 Figure 4.22: Source time function of (a) transducer-laser test and (b) transducertransducer test obtained from implementing deconvolution in frequency-domain. The green line indicates the theoretical first arrival time.

Figure 4 .

 4 Figure 4.23:Source wavelet of transducer-laser test and its FFT spectrum.

Figure 4

 4 Figure 4.24:Source wavelet of transducer-transducer test and its FFT spectrum.

Figure 4 .

 4 Figure 4.25: Processes to abstract the exact source time function (S ). Analysing (a) the transducer effects (R), (b) the function combining both the source time function and transducer effect (S) in order to extract (c) the exact source time function (S ) using equation 4.7. (d) is the spectrum of the source time function excited by the oscillator without transducer effects. The green line is central frequency ( f c ) = 1.05 MHz.

Figure 5 .

 5 Figure 5.1 shows all the metal samples used in this study. Three metal samples with different shapes and materials were used primarily for analysis, as indicated in figure 5.1. Since the copper sample has a cylindrical shape with a diameter of 4 cm and

Figure 5 . 1 :

 51 Figure 5.1: Metal samples. Copper cylinder, hollow aluminium and long aluminium were mainly used in this study.

Figure 5 . 2 :

 52 Figure 5.2: Experimental settings of source transducer point (blue), receiver transducer points (green) and receiver laser point (red).

Figure 5 . 3 :

 53 Figure 5.3: Vibrometer OFV-5000 Modular Vibrometer Controller, Vibrometer Sensor Head OFV-503 from Politec Ltd and long aluminium sample.

Figure 5 . 4 :

 54 Figure 5.4: Vibrometer and long aluminium sample. The long aluminium sample is applied with three piezoelectric transducers at one side and laser beam at the other side.

  .22a) and transducer-transducer tests (figure 4.22b), shown in figure 5.6 and 5.7, respectively.

Figure 5 . 5 :

 55 Figure 5.5: Normalized raw experimental waveforms and numerical synthetics using physical parameters of the metals: transducer-laser test on (a) copper, (c) long aluminium and (e) hollow aluminium samples; transducer-transducer test on (b) copper, (d) long aluminium and (f) hollow aluminium samples.

Figure 5 . 6 :

 56 Figure 5.6: Testing the resolution of source time function of transducer-laser test on (a) aluminium with small hole and (b) hollow aluminium samples.

Figure 5 . 7 :

 57 Figure 5.7: Testing the resolution of source time function of transducer-transducer test on (a) aluminium with small hole and (b) hollow aluminium samples.

Figure 5 .

 5 Figure 5.8 and 5.9 compared the experimental waveforms and synthetic waveforms

Figure 5 . 8 :

 58 Figure 5.8: Comparing the normalized experimental data and numerical synthetics of transducer-laser test on (a) aluminium with small hole and (b) hollow aluminium samples.

Figure 5 . 9 :

 59 Figure 5.9: Comparing the normalized experimental data and numerical synthetics of transducer-transducer test on (a) aluminium with small hole and (b) hollow aluminium samples.

Figure 5 .

 5 Figure 5.10: FFT spectrum of transducer-laser experimental data and numerical synthetics of (a) aluminium with small hole and (b) hollow aluminium samples.

Figure 5 .

 5 Figure 5.11: FFT spectrum of transducer-laser experimental data and numerical synthetics of (a) aluminium with small hole and (b) hollow aluminium samples.

Figure 5 .

 5 Figure 5.12: FFT spectrum of transducer-laser experimental waveform on different metal samples. The green and the red lines indicate the two sinks at 3×10 5 Hz and 6×10 5 Hz.

Figure 5 .

 5 Figure 5.13: FFT spectrum of transducer-transducer experimental waveform on different metal samples. The green and the red lines indicate the two sinks at 3×10 5 Hz and 6×10 5 Hz.

Figure 5 .

 5 Figure 5.14: (Top) Original experimental data with 100, 200, 300 and 400 mV as full scale Volts (FSV). (Bottom) The FFT of the data.

Figure 5 .

 5 Figure 5.15: (Top) Original experimental data with 200 mV as full scale Volts (FSV). (Bottom) The short time spectrogram of the data. The dark red area shows the zone of the energy of the signal.

5. 4 . 1

 41 Waveform matchingLaboratory experiments and numerical tests of three different metal samples were conducted and analysed in this study(figure 5.1). The waveforms in figure5.16 were applied with a low-pass filter of cutoff frequency 2 MHz. The waveforms are applied in their realistic amplitude, except that long aluminium and hollow aluminium waveforms were used with amplitude normalisation. All synthetic waveforms' first arrivals can match the observed data; the differences between the observed and synthetic waveforms are more evident in the later coda waves. The differences between the data and synthetic waveforms of the copper sample are shown in waveform amplitude, specifically after the arrivals of shear waves and surface waves (> 135 µs) for both the transducer-laser and transducer-transducer tests. The synthetics of the long aluminium sample fit with data for the first 30 µs after the first arrivals, from the arrival of the direct P waves to the surface waves. In our simulation tests, while the waves reach the antipodes of the source point, the waves separate into two parts and propagate toward the two sharp ends of the sample. Therefore, there is almost no specific feature to be seen after 140 µs in the synthetics of long aluminium because of its longer height. However, in real data, the additional waves are shown in the observed waveforms after 140 µs, which are also shown in the tests of hollow aluminium. It may be caused by the vibration of attached faces between the piezoelectric transducers and the metal samples, except for the additional reflections after surface waves-the data and synthetics of the hollow aluminium match well in the first 30 µs. The waveforms of hollow aluminium show special features compared to a full cylindrical sample. Instead of sharp P waves arrival, the small reflections between the thin aluminium walls arrive as the first arrivals. These observations can only be seen in the transducer-laser test because the transducer-transducer test has a lower sensibility, so that the piezoelectric transducer also acted as a low-pass filter

Figure 5 .

 5 Figure 5.16: Models, data and synthetic waveforms of transducer-laser test and transducer-transducer test for metal samples: (a) copper, (b) long aluminium and (c) hollow aluminium. Waveforms were applied with a low-pass filter with a cutoff frequency equal to 2 MHz. Apart from the waveforms of the transducer-laser test of long aluminium and hollow aluminium samples to which amplitude normalisation was applied, the others are in their original amplitudes.

Figure 5 .

 5 Figure5.17: The cross-section snapshots of waves propagate in the cylindrical copper model. The cross-section coincides with the plane of the line connecting the source and the receiver.

Figure 5 .

 5 Figure 5.18: (a) The waveforms recorded around the surface of the cylinder from 0 • to 360 • and (b) the schematic diagram of the propagating direction of waves.

  .18b) and formulated a theory of waves propagation mode in the cylindrical model(figure 5.18a). By analysing the seismogram, we can even distinguish the first appearances of P P and P S waves around 90 • and 45 • after the direct P wave, respectively. The wave propagation mode can be determined by calculating the arrival times and angles of the curves in the seismograms. When the waves travelling within the models reach the long edge of the sample, the waves would divide into two and propagate toward the two sharp edges. If waves propagate along the surface and encounter the long edges of the sample, the waves would pass through the antipodes along the surface and travel back to the original point. However, after the waves merge at the original source point, they change their direction and propagate straight into the model toward the antipodes again.5.4.3 Kernels of wave propagation schemeP-wave kernel (K α ) was developed to understand the regions which contributed the changes of velocity waveforms. According to the phases identified in figure5.18, we defined seven time windows of: P, P P , P S , S, S P , S S and Rayleigh waves(figure 5.19A).Figure 5.19B and C show the kernels of the waves in the seven windows in 2D and 3D (figure5.19B(a) and C(a)). Direct P wave has a cigar-shaped kernel from source to receiver point. P P and P S waves have sensitivity at 90 • and 45 • , respectively, as observed in snapshots and seismogram (figure5.17).

Figure 5 .

 5 Figure 5.19: (A) Seven windows opened for different phases of waves: (1) P, (2) P P , (3) P S , (4) S, (5) S P , (6) S S and (7) Rayleigh waves. (B)(a) The cross-section of the model in Z direction, (B)(b-h) The P-wave kernels of seven time windows: (b) P, (c) P P , (d) P S , (e) S, (f) S P , (g) S S , and (h) Rayleigh waves. (C)(a) The entire model view, (C)(b-h) The P-wave kernels of seven time windows: (b) P, (c) P P , (d) P S , (e) S, (f) S P , (g) S S , and (h) Rayleigh waves.

  troducing 3D Fréchet kernels to laboratory centimetre-scaled structure, which can convincingly verify the previous observation of seismograms and snapshots. Except for the direct arrivals, we further identify and confirm more phases by their arrival times and kernels, showing us their travelling characteristics. The processing of this considerable number of wave paths could drastically increase the resolution of the laboratory waveforms and tomography for future studies. Using the kernels obtained by each time window of the waveforms, we can predict the time for the wave arrivals in observed waveforms (figure 5.19A). In figure5.19A, the first P-wave arrival and

Figure 6 . 1 :

 61 Figure 6.1: Rock samples (left to right): Carrara marble (CM), Indian gabbro (IG) and Westerly granite (WG).

Figure 6 . 2 :

 62 Figure6.2: Workflow describing the Monte-Carlo inversion of full waveforms to obtain the best parameters of (V P , V S , Q P , Q S , δ m, g).

Figure 6 . 3 :

 63 Figure 6.3: Random heterogeneity models of δ m = 20% with meshing space equals to (top) 5×10 -4 m :(a-c) g=1, g=10, g=100, (middle) 7×10 -4 m (d-f) g=1, g=10, g=100, (bottom) 9×10 -4 m (g-h) g=1, g=10, g=100.

Figure 6 . 4 :

 64 Figure 6.4: Fourier spectrum of transducer-laser experimental waveform on different rock samples. The green and the red lines indicate the two sinks at 3×10 5 Hz and 6×10 5 Hz.

  3.2) for such properties is due to the experimental configuration and cylindrical sample geometry. Unlike the waveform spectrum of the metal samples, the two sinks in the spectrum of waves propagating in the rock are relatively insignificant.Time-frequency spectrogram was implemented on the experimental waveforms ofWesterly granite, Carrara marble and Indian gabbro (figure 6.6). The highest values of the spectrograms are concentrated below 1.5 × 10 6 Hz. Therefore, the low-pass filters for the experimental waveforms of rock samples were designed to have a cutoff frequency of ∼ 1.5 × 10 6 Hz for data processing.

Figure 6 . 5 :

 65 Figure 6.5: Fourier spectrum of transducer-transducer experimental waveform on different rock samples. The green and the red lines indicate the two sinks at 3×10 5 Hz and 6×10 5 Hz.

Figure 6 . 6 :

 66 Figure 6.6: Spectrogram of the experimental waveforms of rock samples: (a) Westerly granite, (b) Carrara marble and (c) Indian gabbro.

1 .

 1 15677×10 -9 m 3 as shown in figure 6.8 and 6.10, which used two different perturbations of model A and model B. The comparison of differences between the perturbation calculated from partial derivative and the real displacement difference is shown in figure 6.7 (model A) and 6.9 (model B). The two examples of the displacement waveform from small perturbation in model A and B calculated by partial derivatives show consistency with the real difference in displacement amplitude.

Figure 6 .

 6 Figure6.11 helps us visualise the time-lapse consequence in wavefield caused by a single perturbation at one small point. Another way to visualise the perturbation is using the time-lapse kernels (figure6.12) that can clearly show the perturbation point and its impacting area within the model.

Figure 6 . 7 :

 67 Figure 6.7: Model A for small perturbation and partial derivative test. (a) The entire cylindrical model with a small perturbation (0.1%) in the middle (green). (b) Zoom in of the patch of small perturbation. The volume ∆V is ∼ 1.2×10 -9 m 3 .

Figure 6 . 8 :

 68 Figure 6.8: Comparison between the real displacement difference and the difference obtained from partial derivative in our perturbation model A (figure 6.7).

Figure 6

 6 Figure 6.9: A different model (model B) for small perturbation and partial derivative test, where the small perturbation is different than the one in figure 6.7. (a) The entire cylindrical model with a small perturbation (0.1%) in the middle (green). (b) Zoom in of the patch of small perturbation. The volume ∆V is 1.64×10 -9 m 3 .

Figure 6 . 10 :

 610 Figure 6.10: Comparison between the real displacement difference and the difference obtained from partial derivative in our perturbation model B (figure 6.9).

Figure 6 . 11 :

 611 Figure 6.11: Wavefield difference related to the middle perturbation point of model A (figure 6.7) by time steps of every 4 µs.

  Figure 6.14 shows an example of using δ m =30% as the input parameter of the random heterogeneous model. g values control the number of surrounding meshes given with the same velocity changes.

Figure 6 . 13 :

 613 Figure 6.13: Original models with different grain sizes/element sizes: (a) 4×10 -4 , (b) 5×10 -4 , (c) 6×10 -4 , (d) 7×10 -4, (e) 8×10 -4 , (f) 9×10 -4 , (g) 1×10 -3 ,(h) 1.1×10 -3 and (i) 1.2×10 -3 m. After constructing the original models, the parameters for building random heterogeneity are later added to these models.

Figure 6 . 15 :

 615 Figure 6.15: Effects of δ m in waveforms of granite transducer-laser experiments. From (a)-(h) the models with mesh size from 5×10 -4 m to 12×10 -4 m are shown.

Figure 6 . 16 :

 616 Figure 6.16: Effects of g in waveforms of granite transducer-laser experiments. From (a)-(h) the models with mesh size from 5×10 -4 m to 12×10 -4 m are shown.

Figure 6 .

 6 Figure 6.17: Waveform matching between the transducer-laser data and synthetic waveforms of all δ m and g parameters using: (A) Cross correlation, (B) L1-norm and (C) L2-norm.

Figure 6 . 18 :

 618 Figure6.18: The best-fit synthetic waveforms determined from figure 6.17 and were compared with the experimental data of transducer-laser experiment in order to find the best-fit final result for the Westerly granite sample.

Figure 6 .

 6 Figure 6.19: The final best-fit waveform determined from figure 6.18. The parameters used in the best-fit model are mesh sizes = 11×10 -4 m, δ m = 42 and g = 200.

Figure 6 . 21 :

 621 Figure 6.21: Effects of g in waveforms of granite transducer-transducer experiments. From (a)-(h) the models with mesh size from 5×10 -4 m to 12×10 -4 m are shown.

Figure 6 .

 6 Figure 6.22: Waveform matching between the transducer-transducer data and synthetic waveforms of all δ m and g parameters using: (A) Cross correlation, (B) L1norm and (C) L2-norm.

Figure 6 . 23 :

 623 Figure 6.23: The best-fit synthetic waveforms determined from figure 6.22 and were compared with the experimental data of transducer-transducer experiment in order to find the best-fit final result for the Westerly granite sample.

Figure 6 .

 6 Figure 6.24 shows the waveform matching results using three cost functions. The differences between the regions of best-matched values (dark blue or purple region) using cross-correlation, L1-norm, and L2-norm are more obvious than the transducerlaser test of Westerly granite. L1-norm and L2-norm show very similar results. However, the best-matched values calculated by cross-correlation sometimes appear where L1-norm and L2-norm obtain the less matched values. It may be because the cross-correlation can only focus on matching the waveform phases, whereas the L1norm and L2-norm fit only the amplitudes. After plotting the best results in every models in figure 6.25 and identifying a best-fit waveform from them as shown in figure 6.26, we took as the best-fit values: mesh size = 12×10 -4 m, δ m = 96 and g = 120 from the analysis of the figure 6.24. Unlike the best-fit values of the transducerlaser test of Westerly granite, the best-fit values of the transducer-laser waveform of Carrara marble cannot be directly pointed out by a single cost function in figure 6.24.Cross-correlation and L2-norm have to be taken into account to determine the g value first, and then we have to look for the δ m. Note that the best-fit parameters are situated at a rare point in the 2D image of waveform matching (figure6.24) without drastic opposite values among the three cost functions. By observing the figure 6.25, the coda waves after the principal arrivals introduced by scattering effects are more evident with the increasing of the mesh sizes from 5×10 -4 m to 12×10 -4 m.

Figure 6 .

 6 Figure 6.24: Waveform matching between the transducer-laser data and synthetic waveforms of all δ m and g parameters using: (A) Cross correlation, (B) L1-norm and (C) L2-norm.

Figure 6 . 25 :

 625 Figure6.25: The best-fit synthetic waveforms determined from figure 6.17 and were compared with the experimental data of transducer-laser experiment in order to find the best-fit final result for the Carrara marble sample.

Figure 6 .

 6 Figure 6.26: The final best-fit waveform determined from figure 6.25. The parameters used in the best-fit model are: mesh size = 12×10 -4 m, δ m = 96 and g = 120.

Figure 6 .

 6 Figure 6.27: Waveform matching between the transducer-transducer data and synthetic waveforms of all δ m and g parameters using: (A) Cross correlation, (B) L1norm and (C) L2-norm.

Figure 6 . 28 :

 628 Figure 6.28: The best-fit synthetic waveforms determined from figure 6.27 and were compared with the experimental data of transducer-transducer experiment in order to find the best-fit final result for the Carrara marble sample.

Figure 6 .

 6 Figure 6.29: Waveform matching between the transducer-laser data and synthetic waveforms of all δ m and g parameters using: (A) Cross correlation, (B) L1-norm and (C) L2-norm.

Figure 6 . 30 :

 630 Figure 6.30: The best-fit synthetic waveforms determined from figure 6.29 and were compared with the experimental data of transducer-laser experiment in order to find the best-fit final result for the Indian gabbro sample.

Figure 6 .

 6 Figure 6.31: The final best-fit waveform determined from figure 6.30.

Figure 6 . 32 :

 632 Figure 6.32: The final best-fit waveform determined from figure 6.31. The parameters used in the best-fit model are: mesh size = 12×10 -4 m, δ m = 90 and g = 200.

Figure 6 .

 6 Figure 6.33 presents a relatively random distribution, so that the cost functions cannot provide consistent, nor opposite values for the best-fit waveforms. The scatters in the transducer-transducer experiment provide high-frequency influences (figure 6.34) so that an additional condition may be needed to obtain a proper coda in waveforms of the transducer-transducer experiments.

Figure 6 .

 6 Figure 6.33: Waveform matching between the transducer-transducer data and synthetic waveforms of all δ m and g parameters using: (A) Cross correlation, (B) L1norm and (C) L2-norm.

Figure 6 . 34 :

 634 Figure 6.34: The best-fit synthetic waveforms determined from figure 6.33 and were compared with the experimental data of transducer-transducer experiment in order to find the best-fit final result for the Indian gabbro sample.

Figure 6 . 35 :

 635 Figure 6.35: Comparison of the transducer-laser experimental data for the three rock samples and their synthetic waveform without and after adding scatters by random heterogeneity with different grain sizes. Original homogeneous synthetics are shown in (a) WG, (b) CM and (c) IG. Synthetic waveforms with scatters are shown in (d) WG, (e) CM and (f) IG.

Figure 6 . 36 :

 636 Figure 6.36: Comparison of the transducer-transducer experimental data for the three rock samples and their synthetic waveform without and after adding scatters by random heterogeneity with different grain sizes. Original homogeneous synthetics are shown in (a) WG, (b) CM and (c) IG. Synthetic waveforms with scatters are shown in (d) WG, (e) CM and (f) IG.

Figure 6 . 37 :

 637 Figure 6.37: Scanning electron microscope micrographs of polished and undeformed Westerly granite (J. Aubry, personal communication, 2022).

Figure 6 . 38 :

 638 Figure 6.38: Scanning electron microscope micrographs of polished and undeformed Indian gabbro (J. Aubry, personal communication, 2022).

Figure 6 .

 6 Figure6.39: EBSD observations on postmortem fault sections performed on Carrara marble[START_REF] Aubry | Fault Stability Across the Seismogenic Zone[END_REF].

Figure 7 . 1 :

 71 Figure 7.1: The schematic diagram of the rupture process during triaxial deformation experiment. Applying a stable confining pressure and increasing the axial stress in a sample (a), cracks will grow gradually and finally it forms a fault as shown in (b).

Figure 7 . 2 :

 72 Figure 7.2: Local heterogeneous model with two connected balls.

Figure 7 . 3 :

 73 Figure 7.3: Local heterogeneous model with two separated balls.

Figure 7 . 4 :

 74 Figure 7.4: Local heterogeneous model with multiple balls.

Figure 7 . 5 :

 75 Figure 7.5: Local heterogeneous model with partially separated and partially connected structures using a set of given coordinates as input.

Figure 7 . 6 :

 76 Figure 7.6: Local heterogeneous faulted model.

Figure 7 . 7 :

 77 Figure 7.7: Perspective schematic of local heterogeneous faulted model.

Figure 7 . 8 :

 78 Figure 7.8: Metal samples used for errorbar tests and migration imaging (from left to right): aluminium-copper, hollow aluminium with a big hole and hollow aluminium with a small hole.

Figure 7 . 9 :

 79 Figure 7.9: Aluminium sample as reference sample (left) and aluminium-copper sample (right).

Figure 7 .

 7 Figure 7.10: Aluminium-copper sample and vibrometer and output converter.

Figure 7 . 14 :

 714 Figure 7.14: Normalised waveforms of figure 7.13.

Figure 7 .

 7 Figure 7.15: Separated plot of waveforms in figure 7.13.

Figure 7 . 16 :

 716 Figure 7.16: Using transducer T10, which was situated at position P09, as the source and recorded waveforms by the other 23 stations. The parameter FSV = 200 mV and damping = 3 (50Ω).

  misfit and FPM( f ) is frequency-dependent phase misfit. In figure7.19, the envelope misfits are more significant. The dominant frequency of 2 -5 × 10 5 Hz for the misfit is shown in FEM analysis and the time points containing the most significant differences in TEM analysis.

Figure 7 . 17 :

 717 Figure 7.17: Spectrogram of waveform on reference aluminum cylinder in figure 7.13 (FSV=200 mV, damping=3).

Figure 7 . 18 :

 718 Figure 7.18: Spectrogram of waveform on aluminium-copper cylinder in figure 7.13 (FSV=200 mV, damping=3).There are also differences among the time-frequency misfits between the waves propagating in aluminium and aluminium-copper samples of transducer-laser tests,

Figure 7 .

 7 Figure 7.19: Time-frequency misfit between the the transducer(T10)-laser waveforms of aluminium and aluminium-copper cylinder.

Figure 7

 7 Figure 7.20:Time-frequency misfit between the the transducer(T10)transducer(T17) waveforms of aluminium and aluminium-copper cylinder.

Figure

  Figure 7.21:Time-frequency misfit between the the transducer(T10)transducer(T23) waveforms of aluminium and aluminium-copper cylinder.

Figure 7 .

 7 Figure 7.22: Migration imaging of balls model (figure 7.2) at 48 mm.

Figure 7 .

 7 Figure 7.23: Migration imaging of balls model (figure 7.2) at 50 mm.

7. 3 . 3

 33 Figure7.24 performed the ideal case of migration imaging, which was calculated by synthetic waveforms of a homogeneous aluminium model and an aluminium-copper model with the same geometry as the real sample (figure 7.25). Migration imaging using the real data collected in the laboratory were also implemented, as shown in figure7.26. However, the results shown in figure7.26 are the most preliminary results before signal processing. Such preliminary results suggest that such smallscale waveform data must undergo some methods for numerical artefacts removal and other data processing before they can be used effectively in migration imaging.

Figure 7 .

 7 Figure 7.24: Numerical case of migration imaging on aluminium-copper concentric homogeneous cylindrical model (figure 7.25).

Figure 7 .

 7 Figure 7.25: Numerical model built according to the aluminium-copper sample used in the laboratory experiment (figure 7.8).

Figure 7 .

 7 Figure 7.26: Migration imaging using real data of aluminium-copper concentric homogeneous cylinder (figure 7.8).

Figure 8 . 1 :

 81 Figure 8.1: Eight stations setting for initial test in compression experiment in laboratory.

8. 2

 2 .2 Topic B: Time-lapse migration imaging on the small-scaled rock during triaxial compression experiments After optimising the migration imaging scheme, precise time-lapse migration imaging experiments on rock samples can be applied in the laboratory. It can start step by step from a eight-station experiment (figure 8.1) in order to test the method in a 2D

8. 2 . 4

 24 Topic D: Numerical modelling and full waveforms analysis on interfaces for different types of rocks In chapter 7.2.1, we introduced models with fault structures (figure 7.6 and 7.7).

Figure 1 (

 1 Figure 1 (A)The time-lapse Z-component wavefield of displacement snapshots. The orange cross denotes the source and the green square denotes the receiver. t=(a) 1.4 ,(b) 3.4 ,(c) 5.4 ,(d) 7.4 ,(e) 9.4,(f) 12.4 µs. (B) Numerical modelled P-wave velocities as functions of the differential stresses of dry sample iva1296 using hand-picking and auto-fitting.

Figure 2

 2 Figure2(A) The rebuilt waveform from 3 to 6 microseconds using the relationship between the observed data and the first synthetic waveform. The error between the observed data and the rebuilt synthetics is calculated by Least Square Error (LSE) method. (B) Observed waveform data plots in red dots and the synthetics with response in black with the new velocities in this study cut from 2 microseconds to 7 microseconds without S-wave arrivals.

Figure 3 (

 3 Figure 3 (A) Random media test with the target circle and colourbar for P-wave velocity from 5500 m/s to 4000 m/s, S-wave velocity 3450 m/s to 2450 m/s, and density ρ in the range of 2500 to 2000. The place outside the circle will be set to zeros when modelling. (B) The synthetic waveform of random media applied with response function. (a) P-wave component (b) S-wave component.

Figure 4 (

 4 Figure 4 (A) Models with different width of faults. The surrounding media is 5000 m/s for P-wave velocity and 3124 m/s for S-wave velocity. The fault zone area is 1000 m/s for P-wave velocity and 624 m/s for S-wave velocity. (B) The waveform of the three models with faults in different width and their besi-fit homogeneous waveform. (a) Waveform of the media with thin fault, (b) Waveform of the media with standard fault, (c) Waveform of the media with thick fault.

  evolution of velocity changes of the rock during its deformation. Sarout et al. (2017) estimated the evolution of crack parameters by inverting elastic waves and permeability data. Zotz-Wilson et al. (

  rock sample was roughly in the shape of a cylinder with a diameter of 20 mm and a height of 40 mm. Two piezoelectric transducers were glued at two opposite sides of the rock sample, and were used as aJournal of Geophysical Research: Solid EarthLAI ET AL.

Figure 1 .

 1 Figure 1. The schematic diagram of the geometry of sample, the position of the transducers and the numerical model used in this study. (a) The general view of the experimental setting of the sample. (b) The cross-section of the sample. (c) The distorted hexahedral mesh of the model for spectral element method. The green region is the rock model. The yellow part represents the surrounding low-velocity material. The mesh size is approximately 0.3 mm.

Figure 2 .

 2 Figure 2. The four source-receiver positions tested in the study on the preferred model (Figure1c) and their normalized-synthetic waveform compared to the data. The red circle represents the source position; contrariwise, the green triangle indicates the receiver position. The four source-receiver positions: (a) source and receiver are both sightly inside the inner material, (b) source and receiver are both slightly outside the inner material, (c) source is slightly outside but receiver is slightly inside the inner material, (d) source is slightly inside but receiver is slightly outside the material. The (b) setting of the source-receiver position is determined to be the best.

Figure 3 .

 3 Figure 3. The (a) P-and (b) S-wave spectrogram of experimental data for 0 MPa differential stress.

1 .

 1 Choice of Range of Seismic Parameters for Monte Carlo Inversions

Figure 6

 6 Figure 6 shows the waveform fitting results on dry sample using different cost functions: (a) zero-lag cross correlation, (b) 1 E L -norm and (c) 2 E L -norm. Note that Equation 3 uses reciprocals of objective (b) and (c) when we combine three objective functions together. Regardless of the objective function type, blue color indicates better waveform fits and red color indicates disagreements between observed and synthetics. Among these three cost functions, zero-lag cross correlation can localize best models with less deviation since

Figure 4 .

 4 Figure 4. The snapshot of the P-wave wavefield at different time step from (a) 2 to (f) 7 μ E s, showing P-wave propagation process inside the medium.

Figure 5 .

 5 Figure 5. The snapshot of the S-wave wavefield at different time step from (a) 2 to (f) 7 μ E s, showing S-wave propagation process inside the medium.

Figure 6 .

 6 Figure 6. The results of waveform fitting using different cost functions for (a-c) P-and (d-f) S-wave, respectively. (a), (d) Zero-lag cross correlation, (b), (e) L1-norm and (c), (f) L2-norm.

  difficult to fit than P-wave waveform. Figure8ashows the comparison of the best-fit P E V and S E V models of the dry sample as well as the hand-picking results fromZaima and Katayama (2018), indicating the evolution of the P-and S-wave velocity. For the experiment under dry conditions, the evolution of are similar, and follows three stages: (a) increasing stage, where  b) transition stage,  E V remains small; and (c) decreasing stage, E V drops rapidly. The highest P-wave velocity structure of dry experiment occurs when the difference between compression stress and confining pressure equals to 81-128 MPa. When P-wave velocity reaches the highest, P wave remains stable. until the differential stress reaches its peak value and decreases until   E = 250 MPa, where the S E V starts dropping suddenly. In dry experiment, the difference between the fastest and slowest P E

Figure 7 .

 7 Figure 7. The result of waveform fitting using objective function (Equation 3). The possible trend of change in (a) P-wave and (b) S-wave velocity with increasing differential stress is shown as the high extreme of the scale (dark blue).

  the increasing compression stress are compared in Figure10a. For both dry and wet experiments, attenuation shows a different behavior compared to velocity during the increase of effective stress. For dry experiment, shows that the internal structural changes of the rock have a greater impact on the attenuation of P waves in dry case. For the wet experiment, P E Q varies initially from 130 to 170, and from 170 to 20 ( from 80 to 110, and from 110 to 20 (

Figure 8 .

 8 Figure 8. Evolution of P E V and S E V of (a) dry and (b) wet data. The black and red circle represent the numerical simulated best-fit results of P E V and S E V , respectively. The diamond shapes are the results from Zaima and Katayama (2018), which are calculated from the hand-picked first arrivals.

/Figure 9 .

 9 Figure 9. Ratios of differential velocities to initial velocity ( E dV / 0 E V ), showing the velocity change compared to the original stage at each stress condition.

Figure 10 .

 10 Figure 10. (a) 1 E Q of wet and dry data, (b) Ratios of differential 1E Q to initial 1 E Q ( Q Q

Figure 11 .

 11 Figure 11. The result of (a) P-and (b) S-wave strain waveform fitting of dry data using Equation 3 to find the best-match model with proper P E V , S E V , P E Q , S E Q . The black lines are the best-match synthetic waveforms. The red lines are the experimental data. The y-axis indicate the deformation progress and correspond to a condition of   E . The x-axis represents the time of wave propagation.

Figure 12 .

 12 Figure 12. The result of (a) P-and (b) S-wave strain waveform fitting of wet data using Equation 3 to find the best-match model with proper P E V , S E V , P E Q , S E Q . The black lines are the best-match synthetic waveforms. The red lines are the experimental data. The y-axis indicate the deformation progress and correspond to a condition of   E . The x-axis represents the time of wave propagation.

  waves are sensitive to smaller heterogeneity than P waves, it is natural to have this poorer fit in S-wave waveforms.

Figure

  Figure 13. (a) P E V / S E V ratio change, (b) Relative P E V / S E V ratio change.

  Figure A1 shows the waveforms for (a) the cylindrical rock model with free-surface boundary, without surface cut; (b) a cylindrical model with free-surface boundary, with surface cut; (c) a cylindrical model with silicone layer without surface cut; and (d) a cylindrical model with silicone layer with surface cut. The additional low-velocity layer surrounding the cylindrical rock sample behaves as an ad-hoc modified boundary conditions to absorb some wave energy, while preserving some reflections. Between models (a) and (b) in Figure A1, we can observe the subtle differences caused by the cut sides. However, both of the two models (a) and (b) cannot reproduce the observed waveforms reasonably. Models (c) and (d) have similarity until 5.5 μ E s (the third peak), after which the reflections from the antipodal point is larger than the real data with model (c). Due to the additional reflections back to the inner rock model in (d), the later phases show faster amplitude decay and the synthetics for (d) resemble better to the observed datum than (c). The comparison among the four models demonstrates that the more features added, the more similar the synthetics to the observed waveform. We thus concluded to use model (d) as our boundary conditions.

Figure A1 .

 A1 Figure A1. The four models tested in this study and their synthetic waveform (a) cylindrical model with free-boundary, (b) side-cut cylindrical model with free-boundary, (c) cylindrical model with surrounding silicone jacketing model and (d) side-cut cylindrical model with surrounding silicone jacketing model.

Figure B1 .

 B1 Figure B1. The input wavelet of the rock experiment (top) and the numerical modeling (bottom). In order to avoid the noise of the experimental wavelet, we use the filtered wavelet as the input.

  L'un des moyens sophistiqués de simuler les dispersions dans les roches consiste à introduire une hétérogénéité aléatoire dans les modèles numériques en utilisant les séries de Born. L'hétérogénéité aléatoire génère des ondes diffusées incohérentes lorsque les ondes sismiques se propagent dans un milieu hétérogène. Les ondes diffusées se superposent au signal source, et des formes d'onde complexes sont observées aux stations. Nous pouvons considérer les ondes diffusées comme du bruit dans les formes d'onde observées. Cependant, les ondes diffusées ne sont pas les mêmes que le bruit, apparaissant de manière aléatoire dans les séries tem-

porelles. Les ondes élastiques se propageant à travers différents échantillons de roche à l'échelle du laboratoire avec des niveaux d'hétérogénéité variables ont été

  .4.2 Waveforms and multiphases of waves . . . . . . . . . . . . 5.4.3 Kernels of wave propagation scheme . . . . . . . . . . . . 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.2 Local heterogeneity migration imaging using numerical synthetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.3 Migration imaging . . . . . . . . . . . . . . . . . . . . . . 7.4 Discussion and conclusion . . . . . . . . . . . . . . . . . . . . . . Significant achievements . . . . . . . . . . . . . . . . . . . . . . . 8.2 Ideas for future investigations . . . . . . . . . . . . . . . . . . . . . 8.2.1 Topic A: Optimizing migration imaging process . . . . . . .

	tiques et d'approcher la physique fondamentale de la roche. Cette étude montre que l'hétérogénéité et la taille des grains peuvent être obtenues en faisant correspondre 6 Wave propagation in rocks List of Figures
	les formes d'onde expérimentales aux formes d'onde synthétiques simulées à partir de nos modèles. L'anisotropie des échantillons de roche peut également être calculée efficacement par notre méthode. Notre étude peut être puissante dans la mesure où nous pouvons évaluer les propriétés physiques des roches en envoyant des signaux 8 Discussion and conclusions 6.1 6.2.1 Active-source one-station experiment on different rocks . . 6.2.2 Elastic/anelastic parameters estimation . . . . . . . . . . . 8.1 A Peer-reviewed articles 1.1 Fault zone structures
	Bibliography

Cette étude présente d'abord notre méthode de modélisation numérique tout au long de nos travaux. Nous présentons ensuite l'inversion de forme d'onde de Monte Carlo pour estimer les paramètres élastiques et anélastiques et leur évolution lors d'expériences de déformation de roches en laboratoire. Les données sismiques actives transducteur à transducteur, une source, une station, d'échantillons secs et saturés en eau sont obtenues à partir de Zaima and Katayama (2018). Sur la base des comparaisons, nous avons récupéré les variations temporelles des vitesses sismiques et de l'atténuation des échantillons de roche pendant les expériences de déformation, que nous avons interprétées comme des développements de fissures. L'analyse systématique en comparant les expériences de laboratoire et numériques dans des exemples de métal et de roche sont rapportés plus loin dans cette thèse. Une modélisation numérique a été réalisée pour comprendre les phases des formes d'onde et l'effet des instruments. Nous avons exploré les phases des formes d'onde expérimentales en modélisant le champ d'onde sismique à l'aide de la méthode des éléments spectraux (SEM) et en calculant le noyau de sensibilité. L'hétérogénéité aléatoire a été introduite dans le modèle afin d'inclure les dispersions dans les formes d'onde synthésismiques de source active et en faisant correspondre les formes d'onde expérimentales aux formes d'onde synthétiques générées par nos modèles. Il est fort probable que cette méthode puisse être utilisée pour prédire efficacement le type de roche et servir à l'imagerie structurelle pendant la déformation et la fracturation des roches. Le chapitre 7 présente la méthode d'imagerie de la migration à l'aide de données sismiques de laboratoire et leurs perspectives pour l'exploration des changements dans le temps pendant la déformation des roches. Le chapitre 3 et les chapitres 4 à 6 sont deux articles scientifiques indépendants.

Contents 1 Introduction 2 Numerical modelling of wave propagation in a centimetric cylinder 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Brief review of linearised seismic wave propagation . . . . . . . . . 2.2.1 Green's function, double-couple source, reciprocity . . . . . 2.2.2 Born approximation . . . . . . . . . . . . . . . . . . . . . 2.3 2D numerical modelling: OPT2D and ABAQUS . . . . . . . . . . . 2.4 Spectral element method and viscoelasticity . . . . . . . . . . . . . 2.4.1 Attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . 3 0D monitoring of weakening process from granite deformation 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Inversion problem . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Traxial compression experiments . . . . . . . . . . . . . . 3.3 Numerical models . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Boundary conditions . . . . . . . . . . . . . . . . . . . . . 3.3.2 Source and transfer functions . . . . . . . . . . . . . . . . . 3.3.3 Choice of a range of seismic parameters for Monte Carlo inversions . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Processing of observed and synthetic data . . . . . . . . . . 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56.2.3 Experimental waveform analysis . . . . . . . . . . . . . . . 6.2.4 Born approximation and Fréchet derivatives . . . . . . . . . 6.2.5 Random Heterogeneous models . . . . . . . . . . . . . . . 6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 Waveform matching: Westerly granite . . . . . . . . . . . . 6.3.2 Waveform matching: Carrara marble . . . . . . . . . . . . . 6.3.3 Waveform matching: Indian gabbro . . . . . . . . . . . . . 6.3.4 Heterogeneous and grain sizes models . . . . . . . . . . . . 6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.1 Predicting heterogeneous and grain sizes . . . . . . . . . . 6.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Migration imaging on metal samples 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.1 Local heterogeneity . . . . . . . . . . . . . . . . . . . . . . 7.2.2 Multi-station laboratory experiments . . . . . . . . . . . . . 7.2.3 Seismic migration . . . . . . . . . . . . . . . . . . . . . . 7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.1 Waveform analysis . . . . . . . . . . . . . . . . . . . . . . 78.2.2 Topic B: Time-lapse migration imaging on the small-scaled rock during triaxial compression experiments . . . . . . . . 8.2.3 Topic C: Full waveform inversion on rock samples . . . . . 8.2.4 Topic D: Numerical modelling and full waveforms analysis on interfaces for different types of rocks . . . . . . . . . . .

Table 3 . 1 :

 31 The description of physical parameter of Shin-Etsu KE45W RTV silicone rubber (Properties: Silicone Rubber, 2020)

	Physical parameter	Value
	Density	1050 kg/m 3
	Bulk modulus	1.5-2 GPa
	Poisson's ratio	0.47-0.49
	Shear modulus	0.3-20.0 MPa
	Young's modulus	1.0-50.0 MPa
	conditions successfully mimicked quasi-surface waves generated at the rock-rubber
	discontinuity.

Table 3 . 2 :

 32 Setting range of elasticity parameters and attenuation parameters

	Model	V P (m/s)	V S (m/s)	Density (kg/m 3 )	Q P	Q S
	rock	3700-5800 2660-3550	2650	10-200	10-90
	silicone rubber jacketing	1390	138	1050	6	4
	increment	5-100	5-20	-	5-20	5-10
					Total models: 6000

Table 5 . 1 :

 51 The elasticity parameters used for wave propagation in metals

	Model	V P (m/s) V S (m/s) Density (kg/m 3 )
	aluminium	6206	3073	2700
	copper	4665	2241	8960

Table 6 . 1 :

 61 The elasticity parameters used for wave propagation in rocks

	Model	V P (m/s) V S (m/s) Density (kg/m 3 ) Q P Q S
	Westerly granite (WG)	4250	2850	2650	18 16
	Carrara marble (CM)	5130	2900	2700	14 13
	Indian gabbro (IG)	6400	3350	2950	26 22

Table 6 .

 6 2: Parameters used for building random heterogeneous rock models using data from transducer-laser experiments

	Model	Original mesh sizes	g	δ m(%)
	Westerly granite (WG)	11×10 -4 m	200	±42
	Carrara marble (CM)	12×10 -4 m	120	±96
	Indian gabbro (IG)	12×10 -4 m	200	±90

Table 6 . 3 :

 63 Mesh size of original models and the volume of grain size

	Mesh sizes input (m) Mesh volume of WG (m 3 ) Mesh volume of CM (m 3 ) Mesh volume of IG (m 3 )
	5 × 10 -4	9.1259 × 10 -13	9.0719 × 10 -13	9.0234 × 10 -13
	6 × 10 -4	1.5679 × 10 -12	1.5532 × 10 -12	1.5502 × 10 -12
	7 × 10 -4	2.4873 × 10 -12	2.4902 × 10 -12	2.4649 × 10 -12
	8 × 10 -4	3.6512 × 10 -12	3.6342 × 10 -12	3.6015 × 10 -12
	9 × 10 -4	5.3007 × 10 -12	5.2876 × 10 -12	5.2596 × 10 -12
	10 × 10 -4	7.0449 × 10 -12	7.0449 × 10 -12	6.9658 × 10 -12
	11 × 10 -4	9.5307 × 10 -12	9.4299 × 10 -12	9.3911 × 10 -12
	12 × 10 -4	1.2574 × 10 -11	1.2632 × 10 -11	1.2518 × 10 -11

Table 6 . 5 :

 65 Grain sizes determined from micrographs and from the scattering waveforms matching presented in this study.

	volume (m 3 )	g	Grain size (mm) δ m (%)

Table 6

 6

	.7: Aggregate Properties of the rock samples
	Rock type/ Average type	V	R	H
	Westerly granite (WG)
	Bulk Modulus (GPa)	58.5695 54.193 56.3813
	Shear Modulus (GPa)	35.7615 34.4746 35.1181
	V Pe (km/s)		6.4947	6.4724	6.4835
	δV Pe (%)		32.317 31.2028 31.7599
	Density (kg/m 3 )		2.6497	2.6473	2.6485
	Carrara marble (CM)
	Bulk Modulus (GPa)	73.3	73.3	73.3
	Shear Modulus (GPa)	32	32	32
	V Pe (km/s)		6.49	6.49	6.49
	δV Pe (%)		32.7	32.7	32.7
	Density (kg/m 3 )		2.717	2.717	2.717
	Indian gabbro (IG)
	Bulk Modulus (GPa)	99.1515 88.7765 93.9640
	Shear Modulus (GPa)	48.0925 42.7594 45.426
	V Pe (km/s)		7.3453	7.2641	7.3047
	δV Pe (%)		26.977 24.8993 25.9382
	Density (kg/m 3 )		3.2368	3.1549	3.1958
	Table 6.8: Crack densities calculated from equations 6.10 and 6.11 and the ideal and
	modelled bulk/shear modulus.
	Modulus/Rock type	Westerly granite Carrara marble Indian gabbro
	K m (GPa)	19.1661	40.7796	76.6902
	K e (GPa)	56.3813	73.3	93.9640
	G m (GPa)	21.5246	22.707	33.1064
	G e (GPa)	35.1181	32	45.4260
	ν m	0.09142	0.2652	0.3113
	ν e	0.2421	0.3094	0.2918
	Crack density ρ K (eq 6.10)	0.5984	0.1891	0.05766
	Crack density ρ G (eq 6.11)	0.4329	0.3004	0.2681

Table 1

 1 The Description of Physical Parameter of Shin-Etsu KE45 W RTV Silicone Rubber(Properties: Silicone Rubber, 2020)

Journal of Geophysical Research: Solid Earth

		10.1029/2021JB021873
	LAI ET AL.	6 of 23

Table 2

 2 Setting Range of Elasticity Parameters and Attenuation Parameters

Journal of Geophysical Research: Solid Earth

LAI ET AL.

Acknowledgments

This research is a part of the PhD thesis of Ssu-Ting Lai, which is co-funded by French-National Program MOPGA (Make Our Planet Great Again) and Université de Paris, Institut de physique du globe de Paris. Numerical computations were partly performed on the S-CAPAD platform, IPGP, France. The authors would like to thank the constructive and fruitful reviews of Editor Y. Bernabe, Associate Editor, G. Rossi, J. Sarout, A. Nicolas and anonymous reviewers.

Data Availability Statement

Data is available through Zaima and Katayama (2018), https://doi.org/10.1029/2018JB016377.

was under fluid-saturated conditions, where distilled water was used as a pore fluid.

The pore pressure of the wet experiment was kept constant at 10 MPa during the deformation. The experiment was performed at room temperature with a confining pressure of 20 MPa and a constant strain rate of 1.3 × 10 -6 s -1 .

Figure 3.2:

The schematic diagram of the triaxial compression test in Zaima and Katayama (2018).

The rock sample was roughly in the shape of a cylinder with a diameter of 20 mm and a height of 40 mm. Two piezoelectric transducers were glued at two opposite sides of the rock sample and were used as a source and a receiver, respectively. P waves and S waves were generated and recorded separately using piezoelectric transducers located at the symmetrical sides of the centre point at the edge (figure 3.1a). by matching the data with synthetic waveforms with simultaneous changes of parameters: P-wave velocity (V P ), S-wave velocity (V S ), P-wave attenuation (Q P ) and S-wave attenuation (Q S ) [START_REF] Lai | Rock deformation monitoring using Monte Carlo waveform inversion[END_REF]. We used the source-time function we obtained from calibration (chapter 4) and set a range of parameters for launching the numerical modelling tests in SEM. The best-fit elastic/anelastic parameters of rock samples are shown in table 6.1. The best-fit waveform with average elastic/anelastic values of the samples was determined by minimising the objective function. Hereafter, we introduced a range of values of δ m and gap (g) to build random heterogeneous numerical models for scattering and then added to our best-fit waveforms by introducing Born approximation, where δ m represents the maximum of P-wave velocity perturbation in the random heterogeneity models and g values control the number of surrounding meshes which are given with the same velocity changes. Originally, there were eight numerical models with different mesh sizes from 5.0×10 -4 m to 1.2×10 -4 m. After processing in SEM using GLL number as 6, the mesh volumes of the eight models resulted from 9.0×10 -13 m 3 to 1.2×10 -11 m 3 . Mesh size of the numerical model, δ m and g, simultaneously controlled the random heterogeneous model as shown in figure 6.3. Using the P-wave kernel, we obtained the variations of velocity waveforms caused by the random heterogeneous P-wave velocity changes.

These waveform changes are considered the effects of scattering in the rock samples with different grain sizes and anisotropy that were then applied to our waveforms previously obtained from simulation using best-fit elastic/anelastic parameters. Minimising the objective function between the observed and synthetic waveforms with scattering effects, our inversion method can provide the best-fit mesh size, δ m and g, thus, the best model of the rock sample. The schematic procedure for full-waveform matching using Monte-Carlo inversion and obtaining δ m and g values is shown in figure 6.2.

AppendixA

Peer-reviewed articles

Summary

In order to understand the physics of earthquake generation, we study the microscopic weakening process in a centimetre-scale laboratory experiment. A series of deformation experiment was performed to the granite sample from Aji region, Japan. During the experiments, active seismic wave propagation was measured. This study numerically models the observed waveform in order to infer to the time-lapse elastic structure changes. We first modify the finite difference code on collocated Cartesian introducing zig-zag free-surface discontinuities.

We applied response function estimation to the synthetic waveform to consider all the effect (i.e. 3-D effect, instrumental response, attenuation). We use trial-and-error method to find the best-fit velocity models using full waveform synthetics with response. Then we use the same method to make a series of velocity models for the automated Monte Carlo scheme. Comparing automated velocity models using different methods, we find that L1 norm is the best fitting method. We further tested fault models that we find the trend of the changes of the velocity model after the creation and during the thickening of the fault. This method will detail a microscopic weakening process during the nucleation of earthquakes using full waveforms, which have an impact on fundamental understanding of earthquake physics.

LAI ET AL.

10.1029/2021JB021873 20 of 23 incidence angle equals to 0, P to P transmission coefficient is as high as 1.75, which is coherent to our observation that the direct P waves are amplified when they transmit and reach the interface between the two materials and travel through a low-E Q and low-velocity material.

Appendix B: Source Wavelet Input

The source time function is formally taken from the laboratory experiment input waveform. However, the absolute amplitude of the input signal transmitted to the rock by the piezoelectric transducers is unknown. Therefore, in this study, we estimated the effective source time function by a trial-and-error approach. We concluded that the recorded input waveform low-pass filtered at 2 MHz is the best source time function to be used. This is the same dominant frequency band used in laboratory experiment (Figure B1).