Jorgen Knut

Eirik Grude Maloy

Guillaume Flekkoy

Megan Daniel

Alain Steyer

Valérie Vidal

Takahiro Hatano

Loic Vanel

Semih Turkaya
email: turkaya@unistra.fr

Renaud Toussaint

Fredrik Kvalheim Eriksen

Guillaume Daniel

Eirik G Flekkøy

Knut Jørgen Måløy

Olivier Lengliné

Numerical Studies of the Acoustic Emissions during Pneumatic Fracturing

Keywords: fracturing, lamb waves, acoustic emissions, power spectral evolution, Hele-Shaw cell Pneumatic Fracturing,, Acoustic Emissions,, Signal Localization,, Microseismicity

Merci!

After the three years (and two months) the first word in my mind to start this part is "Thank you!". Because, starting from the first day in Strasbourg I enjoyed (or at least I did my best to enjoy) my life in this small but heartwarming city. I talk a lot about my adventure in Europe but maybe I did not acknowledge enough the people who inspired me to dive into this life. One of the days I took a lecture in "Introduction to Soil Mechanics" my professor K. Önder Çetin in Ankara told me that it is absolutely possible to study Master's Degree in Europe. His words were in my mind, I was excited with the thought of studying in Europe and I saw my advisor Yalin Arici. He told me that it is possible, but it is difficult to find a scholarship. He encouraged me to make an effort on my French so that one day I may use it in the life I want to live. And most importantly I want to thank Güney Özcebe for recommending me to get a scholarship for my first adventure in Europe

ii feel emotional but I am keeping my words for the toast after the defense. However, I should say that thank you very much for being there (not only professionally) everytime I needed.

Everyone in the lab and in FlowTrans (cheers!) that we spent time together particularly Fredrik, JR, Munem, Cécile, Uli, Maxime, Jamie, Gustavo thank you very much for sharing good moments! It is not polite to finish an acknowledgment without mentioning the family. I want to thank my family for supporting me during all this period.

Teşekkürler.

Chapter 1

Resumé Contents

Problématique et Enjeux

Depuis de nombreuses années, la compréhension de l'influence des fluides dans la déformation des roches est un enjeu majeur pour de nombreux scientifiques. Cette question est présente dans de multiples disciplines scientifiques dont la physique et les géosciences avec notamment l'évaluation des risques naturels (e.g. tremblement de terre, glissement de terrain, éruption volcanique) mais également dans des applications industrielles telles que le piégeage du CO 2 .

Prises séparément, ces 2 disciplines sont bien comprises mais les interactions entre elles restent plus complexes à élucider. De grands thèmes présentent un intérêt particulier tel que l'étude des interactions entre les parties solide et fluide, les variations rapides de la porosité dues à la compaction, les transferts en milieu poreux et la fracturation résultant de l'écoulement à grande vitesse d'un fluide.

Dans le cadre de ce projet de recherche, notre attention s'est focalisée sur l'observation et la modélisation des mécanismes d'interaction entre les parties fluide et solide lors de l'aérofracturation d'un milieu poreux.

La Fracturation dans un Milieu Poreux

Afin d'étudier la fracturation d'un solide poreux liée à l'écoulement d'un fluide, notre dispositif expérimental reproduit le phénomène à basse pression. Le système utilisé comprend une cellule d'Hele-Shaw rectangulaire dans laquelle se trouve un milieu granulaire compacté soumis à des écoulements de fluide interstitiel (voir Fig. 1.1)et il permet l'analyse optique. Ce milieu granulaire est constitué de billes de polystyrène non expansées et de granulométrie contrôlée (diamètre de 80 µm à 1% près). La cellule est fermée de manière imperméable sur 3 de ses bords alors que le 4ème est quant à lui semi perméable, il laisse passer l'air mais pas les grains. L'écoulement d'un fluide est modélisé par une injection d'air à haute vitesse et dont la pression d'entrée est imposée, constante, et contrôlée par une électrovalve reliée à une chambre de compression à l'entrée de la cellule.

L'utilisation d'une caméra rapide (à 1000 images par seconde) ainsi que de capteurs acoustiques enregistrant les vibrations au sein du dispositif nous a permis de caractériser, de manière très fine, l'évolution de la cellule initiale soumise à l'injection d'air à haute vitesse.

Combiner le setup optique et les données acoustiques nous permet de corréler les caractéristiques des activités microsismiques, le mode de fracturation/glissement et la forme de la source déterminée en visualisant la déformation durant ces événements dans la cellule de Hele-Shaw. Grace à ce setup il est possible d'améliorer la compréhension de ce phénomène invisible en milieu naturel car ayant lieu dans le sous-sol (e.g. activités microsismiques dans les réservoirs, miroirs de glissement durant un glissement de terrain, éruption volcanique). Les émissions acoustiques enregistrées sur les parois en verre de la cellule ont été analysées temporellement et dans le domaine de Fourier (Figure 2).

La Fracturation dans un Milieu Poreux 3

Ainsi nous avons pu observer que l'écoulement de l'air creuse des canaux dans notre milieu granulaire et peut, dans certaines conditions, mener à la fracturation de ce dernier et changeant ainsi la perméabilité globale du milieu. L'expérience a été renouvelée de nombreuses fois afin d'éliminer les biais expérimentaux mais également afin de tester l'influence de la pression d'injection et la porosité. Au cours des expériences, les émissions acoustiques générées par plusieurs phénomènes ont été enregistrées. L'évolution du milieu est observée et comparée avec les résultats accélérométriques grâce aux photos enregistrées par la caméra rapide. Nous avons analysé les signaux enregistrés dans le domaine temporel et le domaine de Fourier (voir Fig. 1.3). Pendant l'analyse des données, on choisit une fenêtre temporelle de 5 ms dans le signal acoustique enregistré durant l'injection de fluide et on calcule le spectre de puissance. Les spectres calculés pour plusieurs fenêtres temporelles, sont comparés entre eux pour voir l'évolution des pics caractéristiques pendant que le milieu évolue. Lors de l'apparition des canaux vides, leurs tailles diminuent depuis leurs origines vers leurs extrémités. (voir Fig. 1.2A,B), le spectre de puissance est plus fort dans les basses fréquences (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20). Et puis, quand les canaux se propagent et deviennent plus fins (voir Fig. 1.2C,D) la puissance dans les basses fréquences (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) est diminuée et décalée vers les plus hautes fréquences. Par ailleurs, ce phénomène est observé dans les simulations numériques et dans les activités microsismiques à grande échelle. Initialement, le spectre de puissance a une décroissance qui suit la loi de puissance (y = ax b ). Quand le milieu est fracturé, la diffusion de la pression de fluide dans la zone encore compactée devant les canaux déclenche des évènements acoustiques de type stick-slip similaires aux évènements microsismiques. Les évènements microsismiques ont un spectre de puissance caractéristique. Au lieu d'avoir une décroissance en loi de puissance, il y a plusieurs pics dans différentes bandes de fréquences. Les pics dans le spectre de puissance sont influencés par les interactions entre la partie solide et la partie fluide. L'air injecté dans le système vibre dans les zones vides et quand la pression devient suffisante l'air déforme la partie solide ce qui cause la compaction dans le milieu poreux, la création des fractures et la formation des canaux. Par ailleurs, on a vu que le nombre des évènements acoustiques (stick-slip) correspond à la loi d'Omori, qui a été développée a priori pour les aftershocks de tremblement de terre. Les paramètres de la loi d'Omori modifiée (p = 0.5) sont compatibles avec les paramètres définis pour les évènements volcaniques et les tremors. A partir d'un modèle théorique de la diffusion adapté aux conditions de la cellule de Hele-Shaw, on a montré que la diffusion de la pression de fluide dans le milieu prédit effectivement cette équation d'Omori modifiée et ses paramètres. Il est ainsi possible de prévoir le nombre d'évènements acoustiques durant les expériences. A partir des résultats issus du dispositif expérimental, nous avons publié l'article "Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium.", présenté dans le chapitre 3 [Turkaya 2015]. Dans les simulations numériques, on a étudié les vibrations sur les plaques en verre de la cellule qui sont générées par les variations de forces sur celles-ci. Ces forces sont majoritairement de 3 types différents : les vibrations d'air dans les canaux, la diffusion de pression de fluide dans les pores et le transfert de matière (air) au sein du milieu poreux dans la cellule. Nous avons développé, de manière théorique, une équation générale qui permet de calculer l'évolution de la pression de fluide (air) dans la cellule. Pour ce qui concerne la contrainte de la partie solide nous avons utilisé la théorie de Janssen développée pour les silos. La détermination et la quantification des variations des forces au sein de la cellule permettent la modélisation des émissions acoustiques et de l'expression de ces ondes sur une plaque développée par Goyder et White 1982[Goyder 1980].

Localisation des Emissions Acoustiques

La détermination de la position de la source d'un signal est un problème dans de multiples domaines. Par exemple, la détermination d'une position sur les écrans tactiles, la localisation des épicentres des tremblements de terre sont des applications de cette technique. Dans une cellule de Hele-Shaw, la position de la source d'un événement acoustique peut nous indiquer la source mécanique de la déformation fragile dans le milieu poreux. Nous avons travaillé sur la localisation des ondes sur une plaque en verre (les mêmes plaques utilisées dans un setup d'une cellule d'Hele-Shaw). Nous avons publié un article sur une nouvelle méthode de localisation des ondes sur les plaques [Turkaya 2016]. Dans cet article (voir chapitre 5), on explique cette méthode de localisation, développée pour les ondes de Lamb dans les plaques 2D. Cette méthode est basée sur le fait que l'énergie de la source doit avoir la même valeur après la correction de la perte d'énergie due à la distance et l'atténuation. L'énergie initiale est calculée en utilisant les signaux enregistrés par plusieurs capteurs positionnés sur différents points de la plaque. On a réalisé une L'écoulement du fluide a engendré une forte chenalisation et les canaux créés présentent de nombreuses ramifications. La deuxième figure montre l'amplitude du signal enregistré pendant l'expérience. Il nous renseigne sur l'allure du signal tout au long de l'expérience (durée totale 4 secondes). La troisième figure montre un zoom sur la fenêtre temporelle choisie et qui correspond à la photo de la cellule. Le signal enregistré correspond à un événement de type « stick-slip ». En utilisant ce signal, on peut calculer la puissance spectrale associée. La dernière figure montre la puissance spectrale du signal dans la fenêtre temporelle.

Conclusion 7

étude expérimentale et numérique pour bien comprendre les limites et les avantages de cette nouvelle méthode. Pour l'étude expérimentale on a utilisé 4 « Brüel & Kjaer 4374 miniature piezoelectric shock accelerometers » sensible entre 1 -26 kHz (ce qui peut être augmenté jusqu'à 200 kHz avec une correction de réponse de l'accéléromètre) pour enregistrer les signaux générés par une bille en acier qui est lâchée à 2-3 cm de hauteur vers la plaque (en verre ou en plexiglass et de dimensions 80 cm x 40 cm x 1 cm ). De plus, cette méthode a été testée par les simulations numériques dans 4 situations différentes :

1. Plaque infinie (aucun signal n'est reflété) et non-dispersive 2. Plaque finie et non-dispersive 3. Plaque infinie et dispersive 4. Plaque finie et dispersive Les différentes conditions utilisées dans les simulations numériques ont mis en évidence les différentes contributıons induisant des erreurs dans les estimations. Les résultats expérimentaux et numériques sont comparés avec les autres méthodes (Arrival Time Localization), pour avoir une idée de la faisabilité de la méthode. On compare les méthodes sur différentes fréquences d'échantillonnage, géométrie du réseau de capteurs, gamme dynamique, rapport signal sur bruit et temps de calcul. On montre que cette méthode est robuste et facile à adapter dans différentes conditions. L'erreur moyenne d'estimation de la localisation est d'environ 3 cm avec les équipements peu onéreux utilisés. A la fin du chapitre 5, on montre la possibilité d'adapter cette méthode en conditions 3D à des fins industrielles (e.g. microsismicité, robotiques, médecine) ou dans les systèmes naturels (sismologie, volcanologie, glissement de terrain etc.). Cette adaptation peut donner une direction de recherche pour de futurs projets scientifiques.

Conclusion

D'après les expériences et les simulations numériques, nous avons observé que les interactions entre la partie solide et la partie fluide dans une cellule d'Hele-Shaw ressemblent fortement aux explosions volcaniques et à la fracturation hydraulique industrielle. Par ailleurs, nous avons développé une méthode permettant de localiser la source des ondes microsismiques, ce qui peut être utilisé pour estimer la position des sources durant les expériences, et observer conjointement les déplacements dans ce type de source pour améliorer leur compréhension. Cette méthode marche bien en 2D et peut être aussi implémentée en 3D. Dans ce travail le principal objectif est de tester les équipements expérimentaux, l'installation et les procédures pour analyser les données acquises. Par l'utilisation de ce travail préliminaire, nous établissons des règles à suivre pour notre travail. Les principaux buts de cette recherche et de cette comparaison sont clairement indiqués dans la " newsletter proceedings" publiée dans "FlowTrans Review, April 2014 " .

Par plusieurs essais et erreurs, nous avons développé un protocole expérimental avec les étapes à suivre pour reproduire l'expérience et un exemple de la préparation de cette méthode. Quand nous préparons l'expérience dite « test » en suivant les mêmes étapes à chaque fois, nous pouvons reproduire quasiment les mêmes conditions initiales malgré le caractère aléatoire de l'arrangement des particules qui représente toujours une différence. Les propriétés acoustiques de ce dispositif expérimental (bruit initial -Figure 1.4, les émissions acoustiques dues à l'injection -Figure 1.5 etc.) sont contrôlées. Différentes étapes sont identifiées en utilisant les données optiques puis les propriétés caractéristiques des données acoustiques sont trouvées, elles sont différentes pour chaque étape identifiée. Deux différentes méthodes permettant la localisation du signal sont comparées en utilisant les signaux acoustiques et montrant différentes positions pour la source.

Pour les données optiques, certaines « astuces » peuvent aider à obtenir un bon contraste entre les structures vidées et le matériel granulaire. De multiples tests ont été réalisés pour optimiser la qualité de l'acquisition des données. Finalement on utilise un carton noir sous la cellule de Hele Shaw pour obtenir une couleur noire homogène sur l'aire découpée. Finalement, on convertit les images enregistrées de la structure vidée en images segmentées binaires. A partir de ces images segmentées on peut rapidement extraire les quantités observables de la structure en utilisant un script. Une partie de ces particules est colorée avec de l'encre de chine pour avoir un contraste entre les particules afin que tous les déplacements possibles soient traqués pendant la mesure de la vitesse des particules réalisée en analysant de multiples images en utilisant « Particle Image Velocimetry » (PIV).

Pour conclure, les procédures expérimentales développées permettent d'avoir une expérience reproductible d'aérofracturation dans une cellule d'Hele Shaw. En utilisant ces résultats, c'est possible d'obtenir de nombreux jeux de données optiques et acoustiques à étudier, à comparer et éventuellement à utiliser pour comprendre les mécanismes sous-jacents à cette complexe interaction fluide solide. Par ailleurs, dans les travaux suivants, nous souhaiterions localiser la source de ces évènements (optiques et acoustiques) et éventuellement comparer ces 2 différentes méthodes d'enregistrement. Les enregistrements des signaux sont divisés par intervalle de 5 ms, puis l'évolution de la signature du spectre de puissance est analysée dans le domaine de Fourier.

On détermine que le domaine des basses fréquences diminue avec la fracturation du milieu. Comme la zone chenalisée augmente, le pic de la portion des basses fréquences du spectre de puissance devient plus petit et se décale vers les hautes fréquences.

Quand le réseau de canaux atteint sa forme finale (environ 90%), on commence à voir différents types de structures dans le spectre de puissance. Il y a des évènements discrets et impulsifs qui se produisent à cause de la contrainte du réarrangement à cause de la diffusion de la pression de pore dans le système. Ces évènements sont similaires aux répliques de séismes (aftershocks). Le taux d'apparition de ces évènements diminue à mesure que l'injection d'air continue. Dans Figure 1.6, on voit que ces événements (que nous allons appeler de type 2) ont un spectre à large bande dans lequel l'énergie n'est pas concentrée au niveau des basses fréquences mais se propage sur une large gamme de fréquences. nombre de répliques (aftershocks) en utilisant une loi modifiée d'Omori avec les données expérimentales et microsismiques. Nous montrons ainsi que l'exposant de la loi modifiée d'Omori p 0.5 (voir Fig. 1.7), est compatible avec un processus contrôlé par la diffusion de la surpression, suivant un rapide processus de fracturation -le temps de départ de la loi modifiée d'Omori coïncide avec le temps de croissance du canal principal.

L'Article: Numerical Studies of the Acoustic Emissions during Pneumatic Fracturing

Dans la section précédente, nous avons analysé les données expérimentales du système d'acquisition en utilisant les capteurs acoustiques et une caméra ultra rapide. Même si les outils mentionnés nous donnent beaucoup d'informations à propos du système, il reste quelques questions auxquelles ces analyses ne répondent pas (les vibrations de l'air à l'intérieur du système, les contributions des différents pics du spectre etc). C'est possible d'aborder quelques unes de ces questions en utilisant une approche numérique. Dans cette section, nous expliquons comment obtenir davantage d'informations à propos de ces inconnus en utilisant des modèles numériques, avec la physique basée sur le dispositif expérimental et les observations. Pour établir le modèle numérique, d'abord nous divisons la contrainte totale en 2 contributions : (1) la pression de fluide et (2) la contrainte solide. (1) Pour simuler l'évolution de la pression de fluide, nous utilisons l'équation de la conservation de masse pour un solide et pour un fluide dans une cellule d'Hele Shaw. Ensuite, on introduit l'échange des forces entre la partie solide et la partie fluide en utilisant l'équation de Navier Stokes. A cause de la lenteur de l'écoulement du fluide dans les pores, le nombre de Reynold est plus petit que 1. ( 2) Pour la contrainte solide en jeu entre les grains et la plaques du dessus, on utilise un modèle de friction Coulombien, ensemble avec l'équation modifiée de Jannsen pour les "réglages" horizontaux. Après la combinaison de ces 2 zonages de contrainte, on obtient l'évolution de la contrainte totale à l'intérieur de la cellule de Hele Shaw. En utilisant comment les contraintes évoluent à l'intérieur de la cellule, c'est possible de calculer les vibrations. Nous convoluons ces variations de contraintes avec des fonctions de propagation des ondes dans les plaques de verre constituant les côtés de la cellule de Hele Shaw, c'est à dire fonction de Green développée pour les ondes de Lamb. Finalement nous comparons le spectre de puissance de ces vibrations synthétiques avec les enregistrements expérimentaux (voir Fig. 1.9).

Après avoir contrôlé les différentes contributions du spectre de puissance, nous trouvons que la variation de la pression de pores due à la vibration de l'air dans les chenaux ouverts domine essentiellement le spectre de puissance. Au début, (t < 0.05 s) la pression de l'air explose à l'intérieur de la zone vide de la cellule, et éventuellement atteint et pousse le milieu granulaire, lequel crée la plus grande contribution dans la zone des basses fréquences (f < 20 kHz). Nous trouvons que le principal pic du spectre de puissance dans les résultats expérimentaux et numériques a une très bonne correspondance. Par ailleurs, nous avons remarqué une décroissance de la loi de puissance dans le spectre de puissance des vibrations de l'air (comme mentionné dans les résultats expérimentaux). Cela pourrait être l'influence de la zone des chenaux qui présente une configuration fractale, n'ayant pas de taille caractéristique dans l'espace.

Par ailleurs, comme les grains sont compactés autour du canal principal, la contrainte solide commence à générer des signaux puissants. Depuis que la contrainte solide a des sauts dus à la soudaine compaction et à la soudaine relaxation de la Figure 1.9: La comparaison de résultats des simulations numériques avec les données acquises pendant les expériences de même type. En A, on observe que les pics maximaux (lignes rouges) entre l'expérience et le modèle numérique ont la même fréquence. En B, on voit qu'on obtient les mêmes pics si on utilise uniquement la pression de l'air présente dans les canaux. En C, les flêches verticales rouges indiquent que l'on observe les mêmes structures des événements de type 2 entre l'expérience et la simulation numérique.

16

Chapter 1. Resumé contrainte (si l'érosion se produit après la compaction), ce chargement est augmenté dans toutes les fréquences de la plaque de verre de la cellule de Hele Shaw, laquelle génère éventuellement un spectre de puissance similaire aux Type 2 observés expérimentalement (comme expliqué dans la section précédente). En comparant directement les simulations avec les résultats expérimentaux, nous obtenons que le taux d'apparition des évènements générés par cette contrainte solide sont similaires aux évènements de type 2 observés expérimentalement.

L'Article: Pneumatic fractures in confined granular media

Les données optiques acquises durant les expériences d'aérofracturation qui sont analysées en utilisant une méthode de traitement d'images. Ce travail a été principalement réalisé par Fredrik K. Eriksen sous la supervision de Renaud Toussaint. Les expériences ont été réalisé avec Semih Turkaya. Egalement de multiples discussions scientifiques pour faire évoluer la réflexion et le travail scientifiques avant la rédaction et la publication d'articles, ont eu lieu avec tous les co-auteurs. Tous les co-auteurs ont contribué grâce à des suggestions de références scientifiques, de méthodes scientifiques pour l'analyse et la comparaison et également des conseils pour la présentation des résultats obtenus.

A la suite de ces analyses, nous avons observé que la formation des chenaux visibles dans nos expériences entraine des structures très similaires aux structures observées dans la nature ou dans des expériences similaires telles que dans une cellule d'Hele Shaw vide. Egalement, on observe une caractéristique de la croissance des canaux, spécifiquement la présence d'une zone de croissance active en dehors de la structure qui est stable, à cause (screening) du gradient de pression par la partie la plus évoluée de la structure. Ce phénomène est typique pour un système à croissance Laplacienne.

Les chenaux résultants (voir Fig. fig:fre1) dans notre expérience sont différents des chenaux formés dans des systèmes similaires qui ont les bords de la cellule de Hele-Shaw ouverts. Dans ce système ouvert, les chenaux sont généralement plus réguliers et n'ont pas plus de 1 ou 2 branches (ramifications), alors que dans les systèmes fermés, nous observons des structures très ramifiées.

Nous notons que ces structures ont des dimensions fractales bien définies et que les dimensions fractales des chenaux sont trouvées comme étant relativement stables le long de la partie principale des structures. La méthode "local box counting" indique que lorsque l'érosion est faible ou nulle, la structure finale se termine par une dimension fractale typique de l'ordre de D = 1. 53 -1.60. Cette gamme est similaire aux résultats des observations des ramifications visqueuses dans un domaine poreux saturé. Enfin, les analyses montrent que seulement si les structures croissent différemment avec différentes gammes de pressions, les rugosités de l'interface air-solide sont très similaires. Porous media is defined as a binary mixture of superimposed but immiscible constituents [START_REF] Boer | Development of porous media theories-a brief historical review[END_REF], one of them being made of solid constituents (the matrix), the other one being one or several fluids or gas occupying the pores between the solid constituents. This mixture behaves differently from its components. This has attracted the attention of many scientists throughout history in very different branches of science. Material science, geotechnical engineering, medicine, chemical engineering, and particularly soil mechanics participated to this wide interest. Leonhard Euler a Swiss physicist in 1762 defined the porous medium with following paragraph [START_REF] Euler | Opera postuma mathematica et physica: anno mdcccxliv detecta[END_REF]]: "All bodies in the world are composed of coarse and subtle matter, of which the former is called the proper matter, because the other, in view of its practically infinitely low density, contributes nothing to their mass. Since this mixing of the two types of matter extends to the smallest particles, the particles of space in which there is no coarse matter are called the pores of the body, and of these there are several kinds, depending on their size, because even the smallest particles are still filled with pores. The larger of these pores are however not only filled by subtle matter, but also contain air, and therefore also a little of the coarse matter, however it is customary not to count this to the proper matter, and in the present context it is irrelevant whether there is only subtle matter or also air. However the salient difference that must be considered amongst the pores of any body is that some have an open passage to the external aether, whilst others are surrounded by coarse matter to such an extent that the subtle matter can nowhere escape."

Even though the fluids and the rocks are relatively easier to understand individually, the coupled behavior of porous medium makes the system a challenge [Bear 1972[START_REF] Sahimi | Mechanics of disordered solids. II. Percolation on elastic networks with bond-bending forces[END_REF]. Furthermore, in the cases when there is a dynamic fluid flow this system becomes a complete puzzle. The dynamic and complex interaction between the flow and the porous media -rapid changes in the local porosity due to the compaction and migration of the porous material, fracturing due to the momentum exchange in the fast flow -make understanding of such a complex system a difficult task [Gidaspow 1994, Kunii 1991[START_REF] Goren | [END_REF], Goren 2011]. As the number of the phases increase, the behavior of the system differs. For example, all ingredients of a dry porous rock (i.e. grains, bonds between grains, empty pores) will have a similar motion as the body is getting deformed. However, a partially saturated soil column which is getting deformed by the influence of a fast fluid flow through it, will be more complex to understand since the solid fluid and possibly gas phases react differently to the applied forces. Although, this coupled system is widely present in the nature (e.g. volcanoes, earthquakes, rivers, tides in sedimentary structures etc.) monitoring them for further analysis and understanding is not always possible with only keen senses. With the use of high technology monitoring devices (e.g. high speed cameras, Computerized Tomography (CT) scans, accelerometers etc.) it is possible to squeeze more information out of even the simplest experiments about what is happening inside.

Another challenge is to interpret the recorded information during the experiments, the method of translating the obtained information into physical expressions. As René Magritte (see Fig. 2.1) wrote in his famous painting of a pipe La Trahison des images, "Ceci n'est pas une pipe" which means "This is not a pipe" in French. The idea behind this work is that the pipe here is not the actual object. It is only an image of the reality which can be deceptive if interpreted wrong and taken as the actual object. This argument is valid for the experimental results as well, finding a correct way to interpret experimentally obtained data should be based on sound physical explanations [Viggiani 2007].

In this research, we are using acoustic emissions (AE) and optical recordings obtained via a high speed camera to understand how a porous medium is evolving during continuous fast fluid injection.

Brief History of Research on Porous Medium

It is very difficult to give an exact date when the first work on the porous media started. However, it is possible to say how this problem and the treatment 22 Chapter 2. Introduction evolved through time. The work on flow in porous media is so extensive that several journals (e.g. Journal of Porous Media, Transport in Porous Media ), or special topic sections of journals (e.g. Frontiers in Physics: Flow and Transformations in Porous Media) are dedicated to it. This introduction attempts to illustrate some of the topics related to the our research. The question about the mechanics of materials has existed since people first started to build: mastering foundations of large buildings or bridges, building cities on drained swamps, or even very heavy structures like cathedrals [START_REF] Al-Taie | Materials and the Style of Buildings used in Iraq during the Islamic period[END_REF], Bruzelius 2004]. The ancient Egyptian people needed to understand sand/fluid mixtures to build pyramids [START_REF] Petrie | The pyramids and temples of gizeh[END_REF][START_REF] Fall | Sliding Friction on Wet and Dry Sand[END_REF], Reshafim 2016]. In this thesis, based on the comments of De Boer, we assume that the scientific "treatment" of mechanics started around 16 th century when Galileo Galilei founded the mechanics as a branch of science [START_REF] Galilei | Galileo Galilei. Discorsi e dimostrazioni matematiche, intorno à due nuove scienze attenenti alla mecanica & i movimenti locali / del signor galileo galilei ... ; con una appendice del centro di gravità d'alcuni solidi[END_REF][START_REF] Boer | Theory of porous media: highlights in historical development and current state[END_REF]. In his work, he mentioned that Aristotle had already thought about the strength of materials. Galileo did experimental studies on the strength of materials and projectile motion. The famous book of Isaac Newton, Philosophiae naturalis principia mathematica [1687], in which he compiles all the information coming from early scientists and the scientists from his era to interpret them in a new way [START_REF] Newton | Philosophiae naturalis principia mathematica[END_REF]]. Starting from 1730's, Leonhard Euler worked on mechanics (e.g. rigid body motion, fluid flow in pipes etc.) and later specifically on porous medium [START_REF] Boer | Theory of porous media: highlights in historical development and current state[END_REF][START_REF] Euler | Mechanica sive motus scienta analytice exposita[END_REF][START_REF] Euler | Scientia navaiis seu tractatus de construendis ac dirigendis navibus.(cum tabulis)-petropoli[END_REF][START_REF] Euler | Sur le mouvement de l'eau par des tuyaux de conduite[END_REF], Euler 1755[START_REF] Euler | Du mouvement de rotation des corps solides autour d'un axe variable[END_REF][START_REF] Euler | Opera postuma mathematica et physica: anno mdcccxliv detecta[END_REF][START_REF] Euler | Theoria motus corporum solidorum seu rigidorum ex primis nostrae cognitionis principiis stabilita et ad omnes motus qui in huiusmodi corpora cadere possunt accommodata[END_REF][START_REF] Euler | Leonhard euler's briefe an eine deutsche prinzessin über verschiedene gegenstände aus der physik und philosophie: Auf's neue nach dem französischen bearbeitet. mit einem supplemente, die neuesten ergebnisse und bereicherungen der physik in briefform behandelnd[END_REF][START_REF] Euler | Sectio secundo de principii motus fluidorum[END_REF]]. Although there were several scientists working on mechanics in several different aspects, it took until the end of 17th century to introduce the first notions of porous medium. Euler described what is a porous medium in detail and gave eventually an example of sponge to give a clear definition [Euler 1789]. Coulomb in 1774 developed a method to calculate the earth pressure exerting on the retaining walls [START_REF] Coulomb | Théorie des machines simples: en ayant égard au frottement de leurs parties et à la roideur des cordages[END_REF]]. He considered a failing soil block as a free body to calculate limiting horizontal forces. Since the failure plane is not uniquely defined, various number of failure planes need to be checked to find the most critical state (i.e. the failure plane which exerts the maximum force on the retaining wall). Following this, Reinhard Woltman introduced the concept of volume fraction while describing the mud in his four volume opus in 1794 [Woltmann 1791]. Despite having a very first step in this area of science, unfortunately Woltman did not continue his research on porous medium. Following Woltman in 1856, a French physicist named Henry Darcy started working on saturated porous solids [START_REF] Darcy | Les fontaines publiques de la ville de Dijon: exposition et application[END_REF][START_REF] Darcy | Recherches expérimentales relatives au mouvement de l'eau dans les tuyaux[END_REF]]. He did many experimental studies with flow through natural sands, eventually discovering the famous "Darcy's Law", which states the velocity of the flow inside the porous medium is proportional to the pressure deviation from hydrostatic gradient in a continuous volume.

Around the same time, Adolf Fick developed his diffusion laws independently [START_REF] Fick | Ueber diffusion[END_REF]]. Developing Fick's work, in 1871 Josef Stefan from University of Vienna discovered the diffusion of gas [Stefan 1871]. However, he did not get into the mechanical properties of these mixtures which was investigated in the same city but in a different university, Technische Hochschule of Vienna by two professors Paul Fillunger and Karl von Terzaghi.

In 1913, Paul Fillunger worked on porous media theory of liquid-saturated porous solids. He worked on the uplift problem of saturated porous medium which comes essentially from the problem in dam foundations [Fillunger 1913, Fillunger 1914, Fillunger 1915, Fillunger 1929, Fillunger 1930, Fillunger 1934[START_REF] Fillunger | Erdbaumechanik? Selbstverl. d. Verf[END_REF]. His formulation, although famous for a while, was not entirely correct and thus was subjected to criticism from the science community in the following years. Fillunger's formulation was corrected by Karl von Terzaghi who was also working on the same domain in 1934 which eventually created some conflict between two scientists [Terzaghi 1923, Terzaghi 1924, Terzaghi 1925b, Terzaghi 1925a, Terzaghi 1933, Terzaghi 1936a, Terzaghi 1936b, Terzaghi 1934]. Apart from all these minor details, discovery of the fundamentals of the mechanics of porous medium (uplift, friction, capillarity, and effective stress theory) can be credited to these two professors from Vienna [START_REF] Boer | Theory of porous media: highlights in historical development and current state[END_REF].

In 1935, a well-known Belgian-American engineer Maurice Biot continued Karl von Terzaghi's concept by introducing one dimensional consolidation theory into three dimensional bodies for linear elastic porous media with fluids having no viscosity [START_REF] Biot | Le problème de la consolidation des matières argileuses sous une charge[END_REF]]. Biot in his work [Biot 1941b, Biot 1941a, Biot 1941c[START_REF] Biot | Bending settlement of a slab resting on a consolidating foundation[END_REF], assumed the following properties for the soil:

1. isotropy of the material, 2. reversibility of stress-strain relations under final equilibrium conditions, 3. linearity of stress-strain relations, 4. small strains, 5. the water contained in the pores is incompressible, 6. the water may contain air bubbles, 7. the water flows through the porous skeleton according to Darcy's law.

From 1956, Biot implemented his work into more complex cases and introduced anisotropy and viscoelasticity of the medium into his theory [Biot 1955, Biot 1957, Biot 1956b, Biot 1956a]. Later in the same year, he published his work about elastic waves in porous solids in which he summarized all his early work on the porous medium [Biot 1956c].

The work on porous medium continued with some experimental studies in a Hele-Shaw cell (named after Henry Selby Hele-Shaw) [Hele-Shaw 1898], consisting in two transparent plates having an infinitesimal gap (quasi 2D). This system allows visualization of Darcy flow through porous medium and has since been widely used in experimental studies.

It has been mentioned [START_REF] Taylor | The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I[END_REF]] and experimentally confirmed [Lewis 1950] that when two inviscid fluids of different densities are moved in a direction perpendicular to their interface, there is an instability in the interface if the less dense fluid is being injected into the dense one from underneath.

In 1958, Saffmann and Taylor, worked on the instability of the interface between the filling fluid initially saturating the porous medium and the second fluid (immiscible with the first one) being pressure driven into it [Saffman 1958]. If the driving fluid is the least viscous of these two fluids, the fluid-fluid interface is unstable. At the end of XX th century a very large number of people investigated this phenomenon, both from a theoretical and experimental point of view [Bensimon 1986[START_REF] Homsy | [END_REF][START_REF] Tryggvason | Numerical simulations of the Rayleigh-Taylor instability[END_REF][START_REF] Difrancesco | [END_REF], Cross 1993, McCloud 1995[START_REF] Coussot | Philippe Coussot. Saffman-Taylor instability in yield-stress fluids[END_REF][START_REF] Batchelor | [END_REF], Mora 2009, Mora 2010].

Following fluid into fluid instabilities, one of the first works on fluid flow into compressible porous media was conducted by Jönsson et al. in 1992. They explained the problem numerically using steady state flow equation into theoretical permeability and compressibility relationships to model the behavior [Jönsson 1992]. In the same year, a similar study, using gas-gas mixtures, was done by [START_REF] Mulder | [END_REF]].

J.C.Geminard and his colleagues worked on the mechanical response of a wet granular layer against shearing in low stress [Géminard 1999]. Following him, the work on granular materials continued with Losert et al. They worked on hardening of the medium when loaded stepwise in shear [Losert 2000]. Gostiaux et al. worked on air injection inside a porous medium and on the shape of bubbles and their movement through the medium while analyzing the changes in the material [Gostiaux 2002]. Making it more complex and getting more out of a similar experiment, Divoux et al. implemented acoustic sensors to the system [Divoux 2008]. They focused on the acoustic emissions when bubbles are bursting on the interface after traveling through a non-newtonian fluid. Acoustic emissions due to bubbles are not limited to this work. Starting from 2006, Vidal et al. worked on the acoustic emissions of the soap bubbles when they are forming a thin film [Vidal 2006, Vidal 2010]. The main idea is to link these events to the naturally occurring events during volcano eruptions. Varas et al. combined some of these studies, by creating a crater on a porous medium surface due to ascending air injection inside [Varas 2009]. They showed that the size of the crater is invariant of gas dynamics but it depends logarithmically on time of injection. Besides the work focusing on the flow in porous medium, there is a large body of work on mechanical waves in porous media [START_REF] Tuncay | Wave propagation in poroelastic media saturated by two fluids[END_REF], Berryman 2000, Wang 2000, Madeo 2010, Yang 2015]. Seismic attenuation due to flow inside the medium is investigated by Pride et al. [START_REF] Pride | [END_REF]]. This phenomenon is further investigated for the cases with patch saturation by Masson et al. [Masson 2011]. In addition to the works on poroelasticity mentioned here, there are plenty of work including the irreversible deformation due to fluid flow (poroplasticity and poroelastoplasticity).

The flow inside a porous medium may slowly deform it as well. Koehn et al. worked on dissolution of uniaxially stressed crystals with the effect of brine. Under the effect of uniaxial compression, the crystals develop dissolution patterns on its free surface [Koehn 2004]. They worked on numerical simulations of structure development in grain boundaries during dissolution-precipitation creep [Koehn 2006]. Soil pore volume change with applied pressure internal and externally is studied by [START_REF] Berli | Deformation of pores in viscoplastic soil material[END_REF] as well [START_REF] Berli | Deformation of pores in viscoplastic soil material[END_REF]]. They presented a viscoplastic model to quantify the rates of deformation.

Root-soil mechanical interactions are investigated to understand soil stability on hill slopes using experiments in laboratory and simulations using fiber bundle model to determine the influence of different factors (e.g. geometry, water content, soil type) [Schwarz 2010, Cohen 2011, Schwarz 2011[START_REF] Giadrossich | Mechanical interactions between neighbouring roots during pullout tests[END_REF]]. This slope model is developed further by Lehmann et al. [Lehmann 2012]. Ruiz et al. worked on the soil penetration requirements for earthworms and/or plant roots to travel through [Ruiz 2015, Ruiz 2016]. They developed a model to estimate the energy requirements and limits for this phenomenon. [START_REF] Capelli | [END_REF] worked on how the acoustic waves are travelling through snow, experimentally on different snow columns to characterize the propagation of reference signals and numerically using Biot's theory [START_REF] Capelli | [END_REF][START_REF] Capelli | [END_REF]. Recently Faillettaz et al. published research on a natural hazard warning system based on acoustic emissions during failure [START_REF] Faillettaz | Codetection of acoustic emissions during failure of heterogeneous media: New perspectives for natural hazard early warning[END_REF]]. This real-time early warning system is based on the elastic waves generated by failure in a medium modeled using fiber bundle model.

The fact that coupling the mechanical response of the porous medium with fluid flow and chemistry is one of the challenging questions needing to be solved and is mentioned by John Rudnicki in the beginning of 21st century [Rudnicki 2000]. On this type of question on the mechanical aspect, Johnsen et al. worked on a coupled fluid-grain flow with very fine grains [Johnsen 2006]. They studied fluid (gas) -grain interactions to understand systems with several applications, as fluidized beds [START_REF][END_REF], Davidson 1995, Lim 1995, McNamara 2000], flow in tubes where dense solid and gaseous parts interplay in a nontrivial way [Baxter 1989, Måløy 1998, Le Pennec 1998, Flekkøy 1998, Gendron 2001, Flekkøy 2001], and ticking hour glasses [START_REF] Wu | [END_REF], Le Pennec 1996]. Or and Tuller worked on the flow in fractured porous medium to determine the hydraulic conductivity of rough surfaces generated by fracturing [Or 2000]. Vinningland et al., starting from 2007, worked on a similar fluid-grain flow in free fall experimentally and numerically where the aim is to study the granular/gas analogue of the Rayleigh-Taylor instability in the flow of such mixtures [Vinningland 2007b, Vinningland 2007a, Vinningland 2010, Vinningland 2012], i.e. how a body of heavy grains falls into a clear fluid. After, in 2007 and2008, Johnsen et al. studied the granular analogue of the Saffman-Taylor instability: in a horizontally placed rectangular Hele-Shaw cell, they injected air from the side of a shorter edge into the porous medium confined on the long edges of the rectangular cell but free to move towards the other shorter edge [Johnsen 2008a, Johnsen 2008b]. In this setup they studied the global characteristics of pattern formations and the dynamic processes during injection of air inside a non-cohesive porous medium. Following this, Niebling et al. [Niebling 2010a, Niebling 2010b, Niebling 2012a, Niebling 2012b] worked on mixed flow and air injection experiments. They worked further on the mixed flow instabilities by adding a different wetting fluid inside the system. Also, they developed a numerical model to explain the fracturing process due to air injection inside a porous medium, similar to the experiments studied in this thesis. There are already some of the numerical 25 work using the continuum approach on this domain was carried out by Rozhko et al. andValko et al. [Rozhko 2007, Valkó 1994]. However, the advantage of the granular/fluid representation is that it drops some assumptions for the continuity of the displacement, and allows to reproduce highly branched patterns observed during fluid injection in solids. Ghani et al. conducted similar research on numerical modeling of hydrofractures in rocks using discrete element method combined with Darcy based pore-pressure diffusion in continuum [Ghani 2013, Vass 2014, Ghani 2015]. The main objective of this model is to represent large scale formation of joints and faults due to fluid flow.

Brief History of the Research on Acoustic Emissions

Acoustic emissions (AE) are naturally occurring phenomena caused by the radiation of stress waves in a continuous medium. They occur when a material endures a sudden stress redistribution due to some changes in its internal structure (such as bond breaking, granular crushing, fracturation in rocks) under the effect of one or more internal or external agents (e.g. chemicals, mechanical forces, temperature) [START_REF] Huang | [END_REF]]. In metals "tin cry" can be said to be the very first acoustic emission heard, which happened during the Bronze Age after learning smelting to purify tin from its form existing in the nature. The earliest work about acoustic emission in the literature was published in London in 1678 under the title of " The Works of Geber" [Holmyard 1928]. In his work, Geber (Jabir ibn Hayyan) noted that "tin" (called Jupiter at that time) releases a "harsh sound" or "crashing noise". Another example of tin cry was mentioned in the work of [START_REF] Mm | Watts' dictionary of chemistry[END_REF]]. They said that "When a bar of tin is bent, a cracking sound may be heard due to crystals in the inner parts of the bar breaking against one another." It took a while until the scientists commented on acoustic emissions coming from the rock samples. In 1924 in Russia, scientists found that during shear deformation in heated rock specimens create a noise similar to a clock "tictoc" [Classen-Nekludowa 1929, Classen-Nekludowa 1927, Classen-Nekludowa 1928]. Then, in 1925, in United States Robert J. Anderson mentioned that the acoustic emissions generated during aluminum alloy yielding have their pitch varying depending on the thickness of the specimen [Anderson 1926b, Anderson 1926a]. Thin sheets produce high pitched (similar to Japanese glass chimes) sounds while thick sheets produce low pitch sounds like grunting. However, in these mentioned works none of these sounds were experimentally recorded.

Acoustic Instrumentation

Instrumentation is a necessity to have more information out of an experiment and possibly to record it for further analysis. Robert Hooke, a scientist from Royal Society of London worked in the second half of 16th century, thought that using devices to increase the range of the senses -particularly hearing -may lead to new inventions in science. For him acquiring information to understand nature can be obtained by three methods [Hooke 1969]:

• By the Help of the Naked Senses

• By the Senses assisted with Instruments, and arm'd with Engines

• By Induction or comparing the collected Observations, by the two preceding helps, and ratiocinating from them Furthermore, he stated his ideas (which eventually lead to the discovery of the acoustic horn and the stethoscope) as: "There may be also a possibility of discovering the internal motions and actions of bodies by the sound they make. Who knows but that as in the watch we may hear the beating of the balance, and the running of the teeth, and multitudes of other noises: who knows, I say, but that it may be possible to discover the motions of the internal parts of bodies, whether animal, vegetable, or mineral, by the sound they make. That one may discover the works perform'd in the several offices and shops of a man's body, and thereby discover what instrument or engine is out of order, what works are going on at several times, and lies still at others, and the like; that in plants and vegetables one might discover by the noise the pumps for raising the juice, the valves for slopping it, and the rushing of it out of one passage into another, and the like."

To hear the sounds propagating in a different medium than air (such as in water or in the ground), Leonardo da Vinci had some ideas as well. He said [START_REF] Vinci | Leonardo Da Vinci[END_REF]: "If you cause your ship to stop, and place the head of a long tube in the water, and place the other extremity to your ear, you will hear ships at a great distance from you. You can also do the same by placing the head of the tube upon the ground, and you will then hear anyone passing at a distance from you."

Under the direction of these bright ideas, the invention of the stethoscope was made by the French physician René Théophile Hyacinthe Laennec, whose study focused on the diseases of chest and heart in 1816 [Lindsay 1973[START_REF] Laennec | A treatise on the diseases of the chest, and on mediate auscultation[END_REF].

Instrumentation to amplify the acoustic emissions was not only limited with the stethoscope. There are experiments in the early 19th century in Russia using optics and acoustics together to see if the optically monitored plastic deformation coincides with 10.000 times amplified acoustic emissions [Classen-Nekludowa 1929, Classen-Nekludowa 1927, Classen-Nekludowa 1928]. The problem with these studies was the lack of organization in presenting the knowledge. The first scientific report on planned acoustic experiment was presented in Tokyo in 1933 [START_REF] Kishinouye | Fuyuhiko Kishinouye. An experiment on the progression of fracture[END_REF][START_REF] Kishinouye | Fuyuhiko Kishinouye. 49. Frequency-distribution of the Ito Earthquake Swarm of 1930[END_REF]]. The research of Kishinouye was focusing on linking the acoustic emissions during deformation of wood (amplified and recorded) to compare with the crustal deformation which is leading to earthquakes. After some more years of research in the same domain, he recorded the very first acoustic emission waveforms using the oscillograms made by him. Around 1940, in the United States, acoustic emissions is used in research about mines to predict and control rock and mine failures [Obert 1941[START_REF] Obert | [END_REF], Obert 1977, Lockner 1993].

Developments in acoustic emission studies created the necessity to use better instruments to achieve higher quality datasets from the same experiments. Thus, in 1950, Mason worked with transducers made of piezoelectric crystals [START_REF] Mason 1948wp Mason | Ultrasonic observation of twinning in tin[END_REF][START_REF] Warren | Piezoelectric crystals and their applications to ultrasonics[END_REF]. In his work he recorded acoustic emissions while deforming specimens on a quartz crystal transducer which were covering a frequency range from 1-2 kHz to 5 MHz. Playing with the broadband frequency data eventually gave a way to the question, "'which part of this frequency band is the real information and which part is the noise?"'. In the same year with Mason, Josef Kaiser, a graduate student at the Technische Hochschule München in Germany, worked in a more comprehensive way to understand acoustic emissions. In his study, he determined which frequency range corresponds which type of process within the material [Kaiser 1950, Kaiser 1952, Kaiser 1953, Kaiser 1957a, Kaiser 1957b]. In addition to that, he studied different materials to compare these behaviors. His most important finding is the irreversibility of the acoustic emissions. He mentioned that the plastic deformation, which is known to be irreversible, causes acoustic emissions which are irreversible as well.

In 1964, Dunegan et al. tried to eliminate the noise by using the instruments more efficiently. They used instruments with a narrow range higher than the audible range. They concluded that the use of a narrow frequency bandwidth between 30 and 150 kHz helps to get rid of mechanical background noise and helps to interpret experimental results in a more correct way [Dunegan 1963[START_REF] Engle | [END_REF]].

The research on monitoring rock failure and fracture development continued with the work of various researchers [START_REF] Scholz | Microfracturing and the inelastic deformation of rock in compression[END_REF], Lockner 1977, Main 1989, Lockner 1992[START_REF] Liakopoulou-Morris | Microseismic properties of a homogeneous sandstone during fault nucleation and frictional sliding[END_REF], Eberhardt 1998, Zang 1998, Lei 2000, Amitrano 2003[START_REF] Baud | [END_REF], Lei 2004]. The unique potential of acoustic emissions to be a non-destructive testing procedure has been seen by [START_REF] Tatro | Acoustic emission from crystalline substances[END_REF][START_REF] Tatro | Acoustic emission from crystalline substances[END_REF]]. More recently, acoustic emissions are compared with digital image correlation to track the crack propagation inside rock samples in laboratory [Hall 2006].

Acoustic emissions are used to track strain heterogeneities in Argilite rocks [Valès 2007]. In 2009, Cohen et al. published their work on Fiber bundle model for multiscale modeling of landslides. They said that it is possible to study the signatures of the acoustic emissions using this model as an early detection system. In 2012, by Michlmayr et al. this model is developed further by including energy bursts and delayed release of acoustic emission [START_REF] Michlmayr | Fiber bundle models for stress release and energy bursts during granular shearing[END_REF]]. Following this work, they continued investigating different mechanisms generating acoustic emissions during shear deformation [START_REF] Michlmayr | Mechanisms for acoustic emissions generation during granular shearing[END_REF].

There are some studies that look at the sources directly, both optically and acoustically, in various problems to characterize the different source mechanisms [START_REF] Grob | [END_REF], Grob 2011, Lengliné 2011b, Lengliné 2011a, Lengliné 2012].

Farin et al. conducted some experimental studies on rockfalls and avalanches where he monitors those phenomena using acoustic emissions [Farin 2015]. There are several articles with the experiments on paper fractures using acoustic emissions created during crack propagation [Stojanova 2013a, Stojanova 2013c, Stojanova 2014, Stojanova 2013b].

Microseismic events in large scale are directly related to the acoustic emissions generated inside a lab scale porous medium as well. Daniel et al. worked on the reasons of the perturbations on the aftershock activity following the Izmit earthquake in 1999 [START_REF] Daniel | [END_REF]]. They also worked on earthquake aftershock triggering using 3d models [Daniel 2008]. Fluid-triggered earthquake swarms in the East African Rift (similar to hydrofracturing based acoustic emissions in rocks) in real scale is investigated by Lindenfeld et al. in 2012[Lindenfeld 2012]. Michlmayr et al. worked on shear induced deformation in porous media [START_REF] Michlmayr | Shear-induced force fluctuations and acoustic emissions in granular material[END_REF]]. In their experimental setup, they recorded acoustic emissions to study the characteristics of shear zone forming in dense porous medium. Fluid front generating acoustic emissions as it progresses through a porous medium is studied by Moebius et al. [Moebius 2012a, Moebius 2012b, Moebius 2014a]. This macroscopic transport motion leads to multiple abrupt pore invasion events generating acoustic emissions having different characteristic properties. Similarly, during an experiment of drainage it is possible to optically detect these pore invasion events [Moebius 2014b].

The fact that the emitted signals are directly related to the phenomenon occurring inside the porous media, makes recording these acoustic emissions a very important monitoring method for these systems. Particularly for the cases which are not visible from outside (i.e. borehole drilling, rock failures, earthquakes) acoustic emissions contain a lot of information about the evolution of the mechanical properties with time. This phenomenon is even possible for large scale emissions like earthquakes. Duputel et al. worked on centroid time-delays to characterize source properties to identify unique earthquakes [START_REF] Duputel | [END_REF]].

In the work reported in this thesis, the main focus is analyzing the evolution of the acoustic emissions in the systems which can be optically monitored such as a Hele-Shaw cell will help linking physical phenomena to the characteristics of the recorded signals. This could eventually lead to the development of a better signal based monitoring techniques for real scale events (i.e. volcano eruptions, oil extraction, CO 2 sequestration)

Brief History of the Signal Localization

Localization of the source of acoustic signals is a widely applicable technology used in many different academical and industrial areas. From robotics to medicine, from telecommunication to earth science it is possible to find applications of signal localization [Gershman 1995, Valin 2003, Elnahrawy 2004, Malioutov 2005, Zhu 2007, Fink 2015[START_REF] Garnier | [END_REF]].

Chapter 2. Introduction

Even though it is trivial to locate impulsive, non-dispersive signals without any reflections [Aki 2002], as the signal and propagating medium get more complex localization becomes more difficult. Mainly we can subgroup the localization methods into three main types. There are localization methods available based on the time difference of arrival of the signal (TDOA), based on the energy of the received signal or based on the time reversal of the signal.

Time Difference of Arrival

TDOA can work in broadband sources but it requires a very accurate measurement of the signal to find even the tiniest differences in time delay. Carter in 1980 published a special issue about different methods to estimate the travel time of an acoustic signals and thus the delay between two different receivers [Carter 1980]. In 1987 Smith worked on closed-form least-squares solutions of time difference localization with noniterative techniques [Smith 1987]. He also included the maximum likelihood into his localization formulation in one case. He made numerical simulations to compare the normalized error obtained from these different methods and their statistical performances. Brandstein et al. focused on finding different error regions for different array setups. He used Monte Carlo simulations in a context of speech in a video-conference to follow the speaker [Brandstein 1996]. Following this work, he also worked on closed form solutions of estimators. He proposed linear intersection estimator method which is robust and fast enough to be used in real time, where search-based location algorithms are inappropriate [START_REF] Brandstein | [END_REF]]. Yao et al. in 1998 worked mainly on beamforming techniques for localizing signals.

He also used time delay localization in his work to compare performances of these two methods [Yao 1998].

TDOA is extremely sensitive to signal to noise ratio. The delay on the arrival time (or simply the arrival time on each sensor) has to be measured (thresholding, modeling the shape of the first arrival etc.) very precisely. In addition to that, signal velocity has to be computed carefully since the dispersion of the medium can make the travel time complicated to calculate. It has to be considered that the geometrical spreading which leads the signal to noise ratio to decrease with the distance traveled. This causes unavoidable uncertainties in estimations [Grabowski 2014].

Time Reversal Localization

Time reversal localization is based on the idea of sending the signals recorded at different sensors back to a point and stack them to see the maximum amplitude of the sum over all points. For the estimation, stacking on a correct source point should give a larger amplitude compared to a false one. This method is tested to locate seismic source location with simple wave velocity models by McMechan in 1982[McMechan 1982]. This work is continued by McMechan and Chang in the next years [McMechan 1983, Chang 1991]. Rietbrock and Scherbaum in 1994 made first attempts of localization using time reversal using seismic waves in local scale.

Reversing signals back to a source in local scales (or lab scale experiments) are easier than the global scale since there are fewer inhomogeneities (or at least a controllable amount of) in the system [Rietbrock 1994[START_REF] Larmat | Time-reversal imaging of seismic sources and application to the great Sumatra earthquake[END_REF]]. Several works conducted by Fink et al. proved that time reversal is a reliable technique to locate signal sources [Ficek 1997[START_REF] Fink | [END_REF], Fink 2003]. Furthermore, Tourin et al. in his work in 2001 pointed out that small inevitable perturbations during recording and re-emitting phases do not create significant estimation error [START_REF] Tourin | [END_REF]]. In 2006, Larmat et al. worked on imaging of seismic sources using time reversal with different numerical models [START_REF] Larmat | Time-reversal imaging of seismic sources and application to the great Sumatra earthquake[END_REF]]. One of the most important advantages is that time reversal localization can easily be applied to very large datasets without picking particular phases in different recordings since it uses the entire waveform [Larmat 2008]. Moreover, Larmat et al. continued working on localization of real tremor data in global scale models [Larmat 2009, Larmat 2010]. Very recently, in 2016, Bacot et al. published their work on application of time reversal by changing the properties of the medium instantaneously to generate the mirror of the waves [Bacot 2016]. In their experimental work they show with water waves that instantaneous wave mirroring in 2D propagation causes waves to refocus on the source point.

Energy Based Localization

Energy based signal localization methods follows the main principle of received signal strength based localization (RSSL), which relies on the fact that the energy decays with distance [Meesookho 2008]. For a brief history of energy based methods, it is worth to mention their predecessors.

In 2000, Bahl and Padmanabhan from Microsoft Research Team worked on a system called "RADAR", a radio frequency based system which is using signal strength information gathered in multiple sensors to triangulate the position of the user [Bahl 2000]. Girod et al. focused on their research to find a better way to estimate the range of these signals, which is essential to have a high quality estimation [START_REF] Girod | [END_REF]]. Another work on RSSL was conducted by Flathagen and Korsnes in Norway. In 2010, they worked on localization in wireless sensor networks using signal strength measured by different sensors by using different schemes to make the localization more accurate [START_REF] Flathagen | Localization in wireless sensor networks based on Ad hoc routing and evolutionary computation[END_REF]].

Going back to energy based localization, Li and Hu worked on energy based localization using a microsensor array [Li 2003]. They mentioned that the potential advantage of energy based method is the requirement of low intersensor communication and robustness against noise and parameter based perturbations. This work is followed by the work of Sheng and Hu in the same year.

They compared several different methods (Maximum Likelihood, Nonlinear Least Squares etc.) to estimate the source position after computing the signal energy [Sheng 2003b]. In 2005, Sheng published another part of his work where he compares the maximum likelihood implemented energy based localization method extensively with the other available energy based methods [START_REF] Sheng | [END_REF]]. In the following year, 32 Chapter 2. Introduction from Michigan University, Blatt proposed a new method instead of maximum likelihood and nonlinear least squares. In this proposed method he mentioned that it is possible to converge to the optimum result extremely fast [Blatt 2006].

Among these mentioned methods, energy based approaches are simpler and more economical to obtain from the signal dataset [Sheng 2003b[START_REF] Hu | [END_REF], Sheng 2003a]. Apart from the types of localization methods explained, methods based on direction of arrival and methods based on beamforming can be added to the list. However they are not in the scope of this study.

In chapter 4, a new approach of localization based on the signal energy is presented. Here in this work, we used an exhaustive search to find the best estimation over a grid map. It is possible to apply this estimated source energy homogeneity (ESEH) method to localize microseismic events in aero-fracturing of a porous medium in a Hele-Shaw Cell [Niebling 2012b, Turkaya 2015] (see chapter 2), or other type of 2D experiments such as fracture propagation experiments in plexiglass [Lengliné 2011b] or crack propagation in paper [Stojanova 2014]. Furthermore, it is possible to implement this method into a 3D medium at a larger scale which then can be used in the industrial applications such as geothermal activities [Cornet 1998[START_REF] Evans | [END_REF], Šílenỳ 2014], or well stimulation [Belayouni 2015].

Motivations and Objectives

In this work we are trying to simplify the complex structures of porous medium, decompose them into its most basic components and try to achieve the knowledge that is between different partitions, as Leonardo Da Vinci said [Da Vinci 2013]: "In order to arrive at knowledge of the motions of birds in the air, it is first necessary to acquire knowledge of the winds, which we will prove by the motions of water in itself, and this knowledge will be a step enabling us to arrive at the knowledge of beings that fly between the air and the wind".

In this research, an experimental setup is prepared to study the aerofractures due to a fast fluid (air in this setup) flow. This flow causes channeling and fracturing inside the porous medium ultimately creates a completely fractured system with increased permeability (see Fig. 2.

3).

The aim is to understand the mechanics of this phenomena and eventually model the soil-fluid interactions in the ground during this aero/hydro fracturing. To understand the mechanical behavior during the experiment mentioned above, we monitor it via taking pictures with a high speed camera (1000 fps) and recording the acoustic emissions using accelerometers and piezoelectric transducers (see Fig 2.4 and chapter 2 for more details). Furthermore, the information gathered from experiments is used to calibrate poromechanics based simulations.

The research is presented in this thesis over three main chapters and a final fourth chapter about our conclusions: This chapter starts with the preliminary findings of the work on Hele-Shaw cell using acoustic emissions and optical recordings during pneumatic fracturing. Description of experimental setup, signal analysis in Fourier domain, source localization of the signals recorded on the Hele-Shaw cell, pattern formation during fracturing, granular particle tracing are described briefly as a summary of the early work. The part of this first section related to the image analysis and pattern formation are the done with essential contributions of Fredrik K. Eriksen -the author of this PhD manuscript participated in the experimental development and manuscript writing. In the second part of the chapter we investigate the relationship between the acoustic emissions generated during air injection inside a Hele-Shaw cell, and the mechanics of porous medium. First, the power spectral evolution of the acoustic emissions during fracturation process is investigated. Then, according to the characteristics (i.e. peak shape, height, frequency) of different power spectra two different types of events are defined (i.e. Type-1 and Type-2). These different events are detected using different key parameters in detection algorithm (i.e. bandwidth and amplitude of the threshold). Following the event detection, the physical explanation of their frequency of occurrence are investigated using the Omori Law and the Modified Omori Law (for details on these theories see section 3.3). These different acoustic emissions (and their evolution in Fourier domain during continuous injection) are compared using the optical images to understand the mechanics of the source creating acoustic emissions. Using theoretical models, the physical causes of these acoustic emissions are investigated and explained. Moreover, a systematic way to estimate the number of acoustic events is developed. Finally, these physical explanations, Modified Omori Law, power spectral evolution are compared with large scale microseismic activities (e.g. 2004 eruption process of the Arenal volcano, Costa Rica). Later in the same chapter, we discuss the numerical simulations to explain power spectral signature of different source mechanisms during air injection inside the porous medium.

Experimental findings obtained in the chapter 2 are investigated and explained more in detail using numerical simulations. Using poroelasticity formulation (see chapter 3 for detailed formulation and background of this theory) the fluid stress map inside the Hele-Shaw cell is obtained. Then using the Janssen's solid stress formulation for silos [START_REF] Janssen | Versuche über getreidedruck in silozellen[END_REF]] modified for horizontal quasi 2D medium in Hele-Shaw cell [Sandnes 2011], the solid pressure exerted on the glass plates of the Hele-Shaw cell is also obtained. Furthermore, the stress variations over time are used to generate acoustic signals on the glass plates using Goyder and White's far field approximation of plate waves [Goyder 1980].

Having acoustic emissions experimentally and numerically, it became possible to compare the different mechanics generating acoustic emissions separately.

First, different stress contributions to the main recorded signal are investigated. Signals that are due to solid stress, air vibration inside the carved area, and air vibration inside the pores are analyzed separately to identify which mechanical phenomenon is dominating the systemic response in different phases.

Eventually, it is concluded that during the phase where the air bursts inside the porous medium (t < 0.05 s) air vibration inside the empty area and pores are effective. Following this, when the compaction occurs, gradient of solid stress generates more energetic signals. Similar characteristics of the power spectral evolution mentioned in chapter 2 were seen in the simulations as well.

As the system gets more stable after the initial burst, the amplitude of the power spectrum diminishes and eventually gets back to the initial noise level.

In the end of this chapter, we have the summary of the mechanical behavior of solid-fluid interactions during fast fluid flow into a porous media. We will conclude this chapter by showing that our findings are corresponding very well with what we observe in the nature. We explain with experimental and numerical results that how the different earthquakes are occurring and what are their characteristics. Furthermore, we will explain how these characteristics are evolving -how they correspond to the physical parameters of the fluid and solid.

• Analysis of mechanical phenomena during experiments using optical images and image processing methods. In this section, we investigate the different patterns in the channel network after fingering and fracturing using the optical data. The variation of these patterns due to the variations in pressure, porosity, and friction are studied. Different parts of these formed patterns are compared to find self-similar shapes. In the second part of this chapter, we examine the deformation of the porous medium optically. Compaction and transfer of the porous medium is observed and compared with the displacement maps obtained using particle image velocimetry (PIV). Finally, a numerical model accounting for the fluid pressure inside the porous medium is developed. This numerical model is used to find the pressure maps corresponding to the displacement (and therefore stress maps) obtained from PIV analysis. Using the theoretical and experimental results, the effect of injection pressure variation inside the porous medium are investigated. This chapter are based on the work mainly done by Fredrik K. Eriksen. We conducted the experiments together and we regularly discussed during the analysis and reporting phases.

• Developed energy based wave localization method for the waves on thin plates generated via simple sources. In this chapter, we present a new technique to localize signals on thin plates. We first experimentally test this new technique and we compare with the other methods that are currently being used. Then, using a numerical model we compared different possible sce-

38

Chapter 2. Introduction narios to find the strong and weak points of our method in comparison with the others. Based on the experimental and numerical studies we concluded that this method is working better than the other methods when it is possible to define a clear primary signal (signal recorded before the arrival of the reflected waves). Under the conditions when signal is noisy or badly sampled, even though there is a loss in resolution, the results obtained using this new technique are reliable. Then in the second part of the chapter, we discuss the application of the mentioned energy based wave localization method on the stick-slip motion based acoustic emissions in a Hele-Shaw cell. In this particular section we use the localization method explained earlier in the chapter. The result of this localization is compared with the optical localization using digital image correlation between the images right before the start and right after the end of the acoustic signal.

Chapter 3

Investigating Solid-Fluid Interactions inside the Hele-Shaw Cell

Introduction

Fluid induced brittle deformation of a porous medium exists in many different areas.

From an espresso machine to volcanoes, from food industry to construction it is possible to see this phenomenon in different areas of everyday life. In this chapter, we explain solid-fluid interactions using an analogue model developed in a linear geometry with confinement and at low porosity to study the instabilities that occur during fast motion of fluid in dense porous materials: fracturing, fingering, and channeling. Our research focused on these complex fluid/solid mechanical systemsin a rectangular Hele-Shaw cell with three closed boundaries and one semi-permeable boundary -using two imaging techniques: optical imaging using a high speed camera (1000 fps), and passive listening using high frequency resolution accelerometers and piezoelectrical sensors.

In the analysis phase, we compute the power spectrum of the acoustic signal in time windows of 5 ms, recorded by shock accelerometers Brüel & Kjaer 4374 (Frq. Range 1 Hz -26 kHz) with 1 MHz sampling rate. We found that the peaks in the low frequency range (f < 20 kHz) diminish while the medium fractures. The evolution of the power spectrum is compared with the optical recordings. The power spectrum initially follows a power law trend and when the channel network is developed, we see stick-slip motions inside the porous medium that generate acoustic events having a characteristic signature in the power spectrum. These changes in the power spectrum are strongly influenced by the size and branching of the channels, compaction of the medium, vibration of air in the pores and the fundamental frequency of the plate which is also investigated in chapter 3. Furthermore, the number of these stick-slip events, similar to the data obtained in hydraulic fracturing operations, follows a Modified Omori Law decay with an exponent p value around 0.5. An analytical model of overpressure diffusion predicting p = 0.5 and two other free parameters of the Omori Law (prefactor and origin time) is developed. The spatial density of the seismic events, and the time of end of formation of the channels can also be predicted using this developed model.

The air injection via an inlet is done during experiments to create fracturing, channeling in this porous medium. As mentioned, these experiments are observed using optical and acoustic measurement devices. To interpret this acquired data more efficiently, numerical simulations are developed using fluid mechanics and poroelasticity. These numerical simulations have two steps: (1) respecting the boundary between the fluid and solid parts (validated with experimental results), solid and fluid forces exerting on the glass plates are computed. This allowed us to compute a state of stress field at a particular instant in time during air injection. (2) This variation of forces in time are convoluted with a green function of Lamb waves on the glass plates, synthesizing acoustic emissions compatible with the experimentally recorded ones.

This chapter has four subsections: First, the preliminary results of the experiments with Hele-Shaw cell are discussed. These experiments are monitored using optical and acoustic data. In this section, the early results, the method of inter-pretation, the possibilities and the expected results are discussed. The part related to the treatment of the images are done by Fredrik K. Eriksen, which is the topic of the following chapter. Second, the acoustic emissions during air injection inside the porous medium are investigated thoroughly. The characteristics of the signals, the acoustic events, their frequency of occurrence are investigated and compared with real scale events. In the third section, we tried to understand these acoustic emissions further using a numerical simulation. In this model, the stress variations are computed inside the Hele-Shaw cell then combined with a Green function to generate Lamb waves on the plate. Finally in the fourth section, the conclusions of what is obtained in this chapter and the possible future work is discussed.

Furthermore, different causes of these mentioned stress changes (air vibration in the carved area, changes in the effective stress due to fluid-solid interactions) are separately analyzed and compared with the experimental results to understand how the response of the porous medium is evolving in time. In the analysis phase, power spectrum of different timewindows (5 ms) obtained from the recorded signal are computed. We found that (as expected from experimental results) in a synthetic dataset as well, the peaks in the low frequency range (f < 20 kHz) diminishes while the medium fractures.

Experimental Setup

In this experiment we are using a Hele-Shaw cell made of two glass plates (80 cm × 40 cm × 1 cm) placed on top of the other, separated by a 1.5 mm distance. Using aluminum spacers between the glass plates (where steel clamps are located) we ensured that this aperture is equidistant across the cell. To protect the glass plates from stress focus on the clamping points, rubber sheets are placed between the clamp screws and the glass plates. For the experiments particular to this work, we completely sealed three boundaries of the Hele-Shaw cell using steel clamps after placing a double-sided rubber sealing tape. In this stage, we fill the cell after placing it vertically and grains are poured inside. Then, we make the fourth boundary as semi-permeable using 50 µm steel mesh which allows fluid to exit the system but keeps the solid grains inside the cell. Following the placement of semi-permeable boundary , we placed the Hele-Shaw cell vertically, in a way that the semi-permeable boundary stays at the bottom side to de-compact the grains and recompact them at the bottom under gravity: this homogenizes the solid fraction through the cell, and ensure the reproducibility of the experiments by reaching similar initial solid fraction in each experiment. Another important goal of this rotation process is to provide a small rectangular buffer empty of grains around the air injection inlet to avoid having point injection force over the medium: this allows to study the subsequent instability, i.e. the destabilization of an initially flat front perpendicular to the average flow. After a couple of rotations, the Hele-Shaw cell is placed horizontally. In these experiments we filled the cell with non-expanded polystyrene grains (80 µm diameter ± 1%) called Ugelstad spheres. The density of the spheres is 1.05 kg/dm 3 .

Chapter 3. Investigating Solid-Fluid Interactions inside the Hele-Shaw Cell

A mass of 170 g of grains is required to fill the cell. By filling the cell with this amount, we have an initial solid fraction of 52±5% which is close to the 57% random loose pack for monodispersed grains. Depending on the aperture, the electrostatic forces and humidity may be effective on the real experimental beads but of course not for theoretical hard spheres, which can lower the solid fraction. Furthermore, the measurement of 1.5 mm plate separation is subject to an error of 10% which has a direct effect on the solid fraction error bar. To have some tracers and texture inside the porous medium (to be used for digital image processing) we colored some of the grains using Indian Ink. These textures and traces (and their movement) help us to identify displacements during optical data analysis. The upper plate of the Hele-Shaw cell has a pressure inlet which is used for injecting the pressurized air. This inlet is located at the bottom end of the plate, midway from the long edges and 3 cm inside from the steel clamp. The solid-air interface is placed 1-2 cm away from the inlet to avoid pressure localization close to the inlet. The injection of constant air pressure, provided by a pressurized air tank, is started and stopped via an electrovalve placed on the pipe very close to the air. This pressure is proportional to a Heaviside function in time (i.e. a step function, inlet pressure is kept constant after the start). This inlet pressure is monitored using a pressure sensor placed on the air inlet, allowing to check if the real injection pressure is within 5% of the value required during the experiment. The vibrations of the glass plates on the Hele-Shaw cell are recorded using two different sensors. We used piezoelectrical sensors and accelerometers. The data recorded on the Euro Physical Acoustics Micro 80 piezoelectric sensors which are mostly sensitive in the range (200 kHz -900 kHz) are amplified with a Signal Preamplifier. The data recorded on the Brüel & Kjaer 4374 shock accelerometers which are mostly sensitive in the range (1 Hz -26 kHz) are amplified with a Brüel & Kjaer Nexus Charge Amplifier -Type 2692-A. The amplified/conditioned signal is transmitted to the computer via a Ni-DAQ mx PCI-6133 acquisition card with multiple channels at 1 MHz sampling rate. In addition, synchronized with the acoustical data, images of the Hele-Shaw cell are taken via a Photron SA5 high speed camera at up to 1000 frames per second, transferred to a computer and stored numerically. A TTL signal is used as a trigger to initiate the injection and the data acquisition via the camera and the acoustic sensors, thus, enabling time synchronization between the apparatus. The ambient lab noise is also recorded for reference and investigated using camera and accelerometer recordings prior to air injection. After recording, the signals are corrected by using the response function of the accelerometers provided by the manufacturer and cross-checked using an independent calibration test at the lab.

3.3. Article: Aerofracture through a double looking glass, mixing optics and acoustics 43

Article: Aerofracture through a double looking glass, mixing optics and acoustics

This study is conducted 50% -50% by Semih Turkaya (analysis and interpretation of acoustic data) and Fredrik K. Eriksen (analysis and interpretation of optical data) under the supervision of Renaud Toussaint. Then, scientific discussions to develop the work during the redaction are done with co-authors. All co-authors contributed with some suggestions on references, scientific methods on analysis and comparison, and presentation of the obtained results. This work was published in the FlowTrans Review, April 2014 ( http://www.flowtrans.net/newsletters/ ). In this work the main objective is to test the experimental equipment, setup and procedures to analyze the acquired data. Using this preliminary work, we obtained a guideline for our work. The main focus of research and the comparison are clearly indicated in this newsletter proceedings. By trial and error, we have developed stepwise experimental procedures and sample preparation methods. When we prepare the experimental samples by following the same steps each time, we can create fairly reproducible initial conditions, although the randomness of the particle arrangement will always account for some differences. The acoustic properties of the experimental setup (initial noise, noise due to injection etc.) are checked. In different stages which are identified using the optical data, the characteristic properties of the acoustic data are found to be different as well. Two different methods of signal localization are compared using the acoustic signals, showing different source positions. For the optical data, the improvements that can help to obtain a good the contrast between the emptied structure and granular material are tested. Finally we used black cardboard below the Hele-Shaw cell to obtain a homogenous dark color on the carved areas. Then, we converted raw image data into segmented binary images of the emptied structure. From these segmented images we can quickly extract observable quantities of the pattern using a script. Some of the particles are dyed with Indian Ink to have contrast between particles so that any possible displacement can be tracked during Particle Image Velocimetry (PIV) analysis using the subsequent images. As a conclusion, the experimental procedures that enable to have reproducible aerofracturing experiments in a Hele-Shaw cell are developed. Using these results, it is possible to obtain several optical and acoustic datasets to investigate, compare and eventually understand the mechanics behind these complex fluid-solid interactions. Furthermore, in the following work, we would like to locate the source of these acoustic events (optically and acoustically) and eventually compare these two different monitoring methods.

Aerofracture through a double looking glass, mixing optics and acoustics

Semih Turkaya1 , Fredrik Eriksen 1,2 , Renaud Toussaint 1 , Knut Jørgen M åløy 2 , Eirik G. Flekkøy 2

Introduction

The characterization and comprehension of rock deformation processes due to fluid flow is a challenging problem with numerous applications in many fields. This phenomenon has received an ever-increasing attention in Earth Science, Physics, with many applications in natural hazard understanding, mitigation or forecast (e.g. earthquakes, landslides with hydrological control, volcanic eruptions), or in the industry, as in CO2 sequestration, production and reservoir exploitation, or borehole stability problems.

Even though fluid flow, rock deformation [1] or granular dynamics [2] are vast phenomena to understand individually [3], the coupled behaviour of flow in deformable porous media with a dynamic fluid flow makes the system even more intricate to comprehend. The dynamic interaction between flow and the porous media, rapid changes in the local porosity due to the compaction, dilation and migration of the porous material, fracturing due to the momentum exchange in fast flow, make understanding of such a complex system a challenge. [4,5] 2 State of the Art There have been several studies carried out recently related to the solid-fluid interactions. Vinningland in 2006 studied the coupling between air and the granular material falling down under the effect of gravity in a Hele-Shaw cell [6,7,8]. From 2006 to 2008, Johnsen studied coupled behaviour by injecting air into non-saturated and saturated medium and observed the multiphasic flow, decompaction and fluidization in a Hele-Shaw cell [9,10,11]. Following these studies, Niebling [12] studied coupling between glycerol/water solution and granular material similar to Vinningland's work. During these works, numerical models have been developed and the results compared successfully with the experiments [9,12,13,14].

In this study, we present an experimental approach to understand coupled behavior by aerofracturing very fine granular medium inside a Hele-Shaw cell. The experiments are monitored using a combination of techniques, both from geophysics and from experimental fluid mechanics to obtain information to understand solid-fluid interaction.

Description of the Experimental Setup

The experimental setup consists of a rectangular Hele-Shaw cell with three closed boundaries and one semi-permeable boundary which enables the passage of the fluid but not the solid particles. In the current stage of the research, air is used as the intruding fluid.The Hele-Shaw cell used in the 4 experiments consist of two rectangular glass plates with a size of 80cm x 40cm which are placed on top of each other. For the initial experiments the structure used is the three boundaries are impermeably sealed and one of the short edges is semi-permeable (i.e. fluid can go out but grains remain inside the cell). The air injection is done from a point which is close to the sealed short edge (2 cm) and on the central line with respect to the longer edges. The Figure 1 [a] shows a Hele-Shaw cell in the mentioned experimental setup. The density of the grains is 1.005g/cm 3 . Depending on the solid fraction used in the experiment the approximate amount of material used is 150g. For image processing, some of the grains are dyed using the Indian ink. This method increases the contrast between the pixels and help to provide better input data to the image correlation methods.

During the experiments, the fluid is injected into the system from the point opposite to the semipermeable boundary so that the fluid penetrates into the solid and makes a way via creating channels, fractures or directly using the pore network to the semi-permeable boundary. In the Figure 1 [c] channeling in a horizontal experimental setup due to air injection can be seen.

Experiments are monitored via a high speed camera (up to 1000 fps) and various acoustic sensors which are sensitive to the different frequency ranges. Therefore, in the end it is possible to visualize the solid-fluid interaction by processing images and acoustic emissions. These experimental data can be used to have information about the mechanical properties of the solid partition by analyzing separately and checking the compatibility frequently between the results obtained.

Seismic Interpretation

During those interactions between the solid and the fluid phase, microseismic and acoustic signals are emitted which are recorded by various types on sensors -namely Piezoelectric Shock Accelerometer (Freq. range: 1 Hz -26 kHz) and Piezoelectrical Sensors (Freq. range: 100 kHz -1 MHz) -placed on the glass plates. To avoid misunderstanding, in the latter parts of this report Piezoelectric Shock Accelerometer will be referred as Accelerometer.

The acoustic sensors are so far primarily used to understand the characteristic elastic properties of the glass plates and the Lamb waves that are propagating in them. As the project continues the acquired signals will be used to have information about the physical interaction between fluid and solid particles. Several controlled signals (using a signal generator, steel ball drop, hammer hit) are produced to understand the response of the experimental setup and verify the scripts for source localization prior to the air injection experiments. After defining the parameters to characterize the glass plates, the procedure involving the analysis of the acoustic emissions due to fractures started. Recorded signals created by channeling and fracturing are compared and investigated further in both time and frequency domains. Moreover, by using different techniques localization of the acoustic emissions are done, compared and developed. Localizing the source can give vital information about the systems for aero-fracturing for the conditions where the visual data is not available.(e.g. aerofracturing in a borehole).

Signal Localization

Here in this project, several methods have been tried out and compared for localization. First, the classical method using the arrival time delay between different sensors is tried out. Arrival times at each sensor are calculated using the numerical modeling in addition to the experimental data. The numerical model is constructed using the experimentally obtained group velocity of the Lamb waves and the model is used to calculate all possible arrival time delays for all possible sources in the medium. The experimentally obtained delay in the arrival time is then compared with those results to find the best fit between the experimental results. Both data is correlated using the least squares error method to find possible points which correspond to the location of the source of the signal received in the receivers.

Localization via signal energy is based on the fact that attenuation of the signal energy in glass plates is homogenous and very low. Therefore the energy recorded every position in the plate should be equal (or at least very close) to each other. Using the fact that the Equation 1 depends on the distance between the source and the receiver, the location can be found with a trial and error procedure. Below in the two flow charts the methods are also explained stepwise. [15] [16]

Es = 2πRhρ V g |a 2 (ω)| 2 ω 2 exp(γR)dw, (1)
where R is the distance between source and the receiver, h is the half-thickness of the plate, a(ω) is the amplitude in fourier domain, V g is the group velocity, ω is the angular frequency, γ is the attenuation coefficient which is taken as 0 for glass plates used in this experiment.

To validate the localization methods, several set of experiments are conducted by dropping a steel ball to a fixed location in the glass plates. The results of the two methods were similar to each other and were giving results with a 2cm precision. As it can be seen in the Figure 3 methods are successful to find the location of the signal source.

A limitation of the arrival time based method relies on the following approximation: the localization using the arrival time delay makes the assumption that the group velocity of the Lamb waves is the same everywhere in the glass plates, although it can vary slightly. On the other hand, a limitation of the other method arises from the following: calculating the energy without taking into account reflecting waves from sides of the glass plates is another difficulty.

For comparison, both methods are applied to the signal received during an experiment to locate the event.

Above in the Figure 4 the comparison of the two localization methods on a signal received during the air injection can be seen. It should be noted that both methods give results very close to each other as in the validation phase. It is not easy to conclude that one is better than the other. This would require a comparison of both results with a third technique, as e.g. the localization results obtained from image processing. However, for now it can be said that arrival time delay localization is faster to compute than the localization using the signal energy.

Pattern segmentation

Granular fracturing in our experiments is seen as regions or channels where the initial granular packing has been displaced by air and emptied of grains. The bottom of the cell is black, so that an emptied region is seen as an area of dark pixels. On the other hand, the granular media is white and is seen as an area of brighter grey pixels. Therefore, the granular packing and emptied pattern, as shown in the top-middle image of figure 5, are in principle easy to distinguish with a pixel intensity threshold.

To highlight the empty structure formed during an experiment as a change in the granular packing we make a modified image sequence where all the frames in the raw-data are subtracted pixel by pixel by the image of the initial packing. The resulting difference images display the emptied structure as a bright region and non-displaced parts as a dark region, shown in figure 5. After highlighting the emptied pattern, the images are converted into segmented binary images. This means that the pixels we want to measure in the images are given the pixel value "1" and the rest of the pixels are given the pixel value "0". The emptied pattern is the interesting part in this case, and can be selected as the pixels having a value above a certain threshold intensity. The bottom image in figure 5 shows an example of a binary pattern image.

Measurable time-dependent quantities of the emptied structure include basic spatial characteristics as indicated in figure 6. These quantities are the total width w(x, t) of the structure at the position x and time t, the area A(t) of the pattern at time t, the number of fingers n(x, t) crossed by a perpendicular line at the position x and time t, and the length l(t) of the most advanced finger at time t. Other optically obtained quantities of the pattern that can be investigated over time are the fractal dimension, the splitting characteristics of fingers and the resting times between events in e.g. stick-slip motion. Figure 6: Indications of the measurable quantities of the pattern: the total width w, area A, length of the most advanced finger l and the number of fingers n. The "1" pixels are white and the "0" pixels are black. The x-direction is from the closed short-side toward the semi-permeable short side, and the time t is the experimental time of the frame.

Granular particle tracing

As we have tracer particles in the granular packing, we can perform an image correlation technique called PIV (Particle Image Velocimetry) on the image sequences. This technique will give us the average velocity field of the tracer particles on the time between two subsequent images. Currently we have not yet performed PIV on a frame-to-frame basis in order to investigate quantities, but the idea is that we can use the displacement fields to estimate and locate changes in packing density, track the front between the compacted/uncompacted areas in the solid and investigate particle movements in the solid. We have tested the compability of the image sequences with PIV by trying the software on selected frame pairs, and by the look of the results we can conclude that the software is able to track the particles. Figure 7 shows an example of such a PIV test result.

A quantity we can measure of the compacted region, i.e. the region of the solid where particles have been displaced, is the length of the compaction. This is measured in the x-direction (flow direction) from the emptied pattern to the displacement front, and indicates how deep the fluid-solid interaction penetrates the solid. It is of interest to compare the length of the compaction region with the length of the emptied structure along lines in various parts of the compaction. A similar comparison can be done between the area of the compaction region and the area of the emptied structure. We can also look at the progress of the longest compaction length over time. Another quantity we can try to investigate is the particle density distribution around the emptied structure. 6 Results

Signal Analysis in Fourier Domain

The acoustic recordings are also investigated by using several techniques. Below in the figures a horizontal Hele-Shaw cell setup is given. During the experiment 1.5 bar pressurized air is injected into the cell during 6 seconds starting from t = 1sec. The changes in frequencies can be seen by cutting pieces in time domain and visualizing them after Fourier transform.

Figure 8 shows the initial situation in the cell. The uppermost plot shows the signal and the red rectangle shows the signal part to be focused in the lower two images. For the plot showing partial signal in Fourier domain after Fourier transformation Hanning window is also applied. Please note that the low frequency amplitude in the Fourier domain is slightly higher. This can be due to the random noise present in the lab environment.

Figure 9 shows the cell and the signal received when the pumping starts. Bulk movement of grains due to air injection creates 4 peaks in the Fourier domain around 1 kHz,30kHz, 50 kHz and 200 kHz.

When channeling stops the only movement is due to the stress redistribution in the grains which creates the micro events which is similar to tremors. Since there is not a bulk movement the peak in the 1kHz range is getting smaller. However, as long as the injection continues the peak exists. Furthermore, it can be seen that the peaks at 30kHz and 200kHz have disappeared when channeling reaches its final form, and only a steady flow of air remains through the system. So far we have tested the experimental equipment, setup and procedures of aerofracturing in a dry, non-consolidated powder in a rectangular and closed Hele-Shaw cell with a fluid-permeable short edge. By trial and error, we have developed stepwise experimental procedures and sample preparation methods. When we prepare the experimental samples by following the same steps each time, we can create fairly reproducible initial conditions, although the randomness of the particle arrangement will always account for some differences.

For the acoustic emissions, we have investigated the acoustic properties of our experimental equipment and evaluated two methods of signal location of signals emitted by different sources. We have also received and investigated signals emitted by fracturing events during experiments.

When we consider the optical data, we see that the contrast between the granular material and emptied structure is good, and we have successfully converted raw image data into segmented binary images of the emptied structure. From these segmented images we can quickly extract observable quantities of the pattern using a script. The tracer particles dyed with India ink also prove to be sufficiently contrasted to be tracked during PIV, so that we can obtain particle displacement fields between frames of the image sequences.

The conclusion is that we have developed experimental procedures that works and have a setup which enables us to perform reproducible aerofracturing experiments. We are able to obtain a lot of optical and acoustic data to investigate and compare.

We have done the initial work to establish experimental methods, and are now ready to conduct repeated experiments with given initital conditions. This will provide us with several sets of optical and acoustic data, which we can analyze separately and also compare to each other. Specifically, we want to compare acoustic signal localization with optical observations near the located source. We also aim to systematically change initial conditions, such as injection pressure and initial packing density, and investigate the effects on the acoustic signals and observed quantities. We will also use a pressure sensor at the air injection site in order to have a reading on the pressure of the intruding air, which is nearly constant in the air-filled region of the cell.

Article:

Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium 57

Article: Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium

This study is mainly conducted by main author under the supervision of Renaud Toussaint. The experiments are co-operated with Fredrik K. Eriksen. Then, scientific discussions to develop the work were done with all co-authors. All co-authors contributed with some suggestions on references, scientific methods on analysis and comparison, and presentation of the obtained results. In this work, as explained in section 3.2, a purpose-built Hele-Shaw cell was designed to enable both optical and acoustic recordings associated with controlled fracturing of a porous medium via air injection. 4 accelerometers are used to record signals emitted during the experiments. These recorded signals were first divided into 5 ms time windows, then the evolution of power spectral signature was investigated in Fourier domain. We determined that the the low frequency content f < 10 kHz diminishes as the medium is fractured. As the channeled area increases, the bump in the low frequency portion of the power spectrum gets smaller and shifts towards higher frequencies. As the channel network reaches its final form (about 90%) we start to see different types of structures in the power spectrum. These are discrete and impulsive events that are occurring due to stress rearrangements induced by pore pressure diffusion into the system. These events are similar to aftershocks. The rate of occurrence of these events are diminishing as the injection continues. These events (that will be called Type-2 events) have a broadband spectrum in which the energy is not focused on the low frequency but spreads on a wide span of frequencies. Both power spectra evolution and diminishing frequency of occurrence are present in real scale microseismic data. In addition, it is possible to estimate the number of aftershocks by using a Modified Omori Law in experimental and microseismic data. We show that the exponent of this Modified Omori Law, p ≈ 0.5, is compatible with a process controlled by overpressure diffusion, following a fast fracturing process -the time of start of the Modified Omori Law coincides with the end of the main channel growth These signatures, inflection point in the average frequency, appearance of impulsive events of high frequency, and starting point of an Omori Law, can be used straightforwardly in large scale microseismic monitoring of fluid injection and well stimulation. The permeability of the medium can also be directly estimated from the prefactor of the Omori Law.

Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium Semih Turkaya1 *, Renaud Toussaint 1 , Fredrik K. Eriksen The characterization and understanding of rock deformation processes due to fluid flow is a challenging problem with numerous applications. The signature of this problem can be found in Earth Science and Physics, notably with applications in natural hazard understanding, mitigation or forecast (e.g., earthquakes, landslides with hydrological control, volcanic eruptions), or in industrial applications such as hydraulic-fracturing, steam-assisted gravity drainage, CO 2 sequestration operations or soil remediation. Here, we investigate the link between the visual deformation and the mechanical wave signals generated due to fluid injection into porous medium. In a rectangular Hele-Shaw Cell, side air injection causes burst movement and compaction of grains along with channeling (creation of high permeability channels empty of grains). During the initial compaction and emergence of the main channel, the hydraulic fracturing in the medium generates a large non-impulsive low frequency signal in the frequency range 100 Hz-10 kHz. When the channel network is established, the relaxation of the surrounding medium causes impulsive aftershock-like events, with high frequency (above 10 kHz) acoustic emissions, the rate of which follows an Omori Law. These signals and observations are comparable to seismicity induced by fluid injection. Compared to the data obtained during hydraulic fracturing operations, low frequency seismicity with evolving spectral characteristics have also been observed. An Omori-like decay of microearthquake rates is also often observed after injection shut-in, with a similar exponent p ≈ 0.5 as observed here, where the decay rate of aftershock follows a scaling law dN/dt ∝ (tt 0 ) -p . The physical basis for this modified Omori law is explained by pore pressure diffusion affecting the stress relaxation.

Introduction

Fluid flow [1,2], rock deformation [3] and granular dynamics [4] by themselves are very large scientific domains to investigate individually [5]. However, the idea of putting them together via a system of deformable porous medium with a fluid flow makes the phenomena even harder to understand. Rapid changes in the porosity of the medium due to fluid flow, channeling and fracturing via momentum exchange with the flow make understanding the mechanics of the system a challenge [6][7][8][9]. Hydraulic fracturing of the ground is a good example for this coupled behavior of solid and fluid phases. First, the pressure of the flow creates fissures and cracks which changes the permeability of the initial rock. Then, a flowing mixture of fine sand and chemicals helps maintain this cracked state by penetrating the newly opened areas. By jamming and/or cementing the newly-formed channels and cracks, possible relaxation after injection is prevented. Thus, a more permeable state of the rock is preserved after injection for various types of industrial applications. Recently, various well-stimulation projects have attempted to use pressurized gas (N 2 , CO 2 ), instead of water, to trigger fracturing within reservoirs for several reasons (e.g., to avoid wasting water, to sequestrate CO 2 , environmental risks due to chemicals etc.) [10][11][12]. In this study, contrary to conventional fracturing methods, the fractures are induced using air injection.

Monitoring, predicting and controlling fracture evolution during hydraulic-fracturing, steam-assisted gravity drainage, or CO 2 sequestration operations is a key goal [13][14][15][16]. One possibility for monitoring is to use generated acoustic emissions during those operations. In the hydraulic fracturing industry, the typical monitoring devices consist of geophones and seismometers. However, the interpretation of the signals during fast deformations of porous media due to fast fluid flow is not simple. Particularly, the measurements of deformations are usually difficult to achieve in an opaque medium, and the source of the seismic waves and acoustic emissions can be complex. The study of microseismicity during well operations is routinely done in the industry, but its interpretation is often delicate [17][18][19].

In this paper, we present an experimental study using a purpose-built setup allowing channeling and fracturing due to fluid flow, where we can observe the deformations optically using a fast camera and transparent setup, and simultaneously record the mechanical waves emitted by the complex channels and fractures created. Both signals, optical and acoustic/microseismic, are then analyzed. They display a complex evolution of the source geometry, and of the spectral characteristics. The experimental setup designed to achieve this consists of a rectangular Hele-Shaw cell filled with 80 microns diameter grains, mixed with fluid (air). The linear cell has three lateral sealed boundaries and a semi-permeable one enabling fluid (but not solid) flow. During the experiments, air is injected into the system from the side opposite to the semi-permeable boundary so that the air penetrates into the solid and at high injection pressures makes a way to the semi-permeable boundary via the creation of channels and fractures -or at low injection pressures, directly using the pore network.

For a similar system of aerofractures in a Hele-Shaw cell, a numerical model was conducted by Niebling et al. [20]. These models were also compared with experiments for further development and validation [20][21][22][23][24]. Same kind of experiments-but without acoustic monitoring-with a Hele-Shaw cell have also been conducted. Johnsen et al. worked on the coupled behavior by air injection into the porous material both in fluid saturated and non-saturated cases to study multiphase flow numerically and experimentally [23,25,26]. Aero-granular coupling in a free falling porous medium in a vertical Hele-Shaw cell was studied numerically and experimentally by Vinningland et al. [27][28][29]. Varas et al. conducted experiments of air injection into the saturated porous media in a Hele-Shaw cell [30,31] and in a cylinder box [32]. Eriksen et al. and McMinn et al. worked on injecting gas into a saturated deformable porous medium [Eriksen et al.,submitted;33]. Sandnes et al. classified different regimes of fingering of porous media in a Hele-Shaw cell [34]. Rust et al. developed a closed-system degassing model using volcanic eruption data [35]. Holtzman et al. also studied air induced fracturing where they identified different invasion regimes [36]. Furthermore, a recent study was conducted by Eriksen et al. where the air injection causes bubbles in a fluid-grain mixture [37].

The equivalent of microseismicity monitoring in the lab is the tracking of acoustic events. Hall et al. compared recorded acoustic emissions with the digital image correlation to track crack propagation in the rock samples [38]. Valès et al. used acoustic emissions to track strain heterogeneities in Argilite rocks [39]. Some studies have started to look at sources directly, both optically and acoustically, in various problems to characterize the different source mechanisms [40][41][42][43][44]. Farin et al. conducted some experimental studies on rockfalls and avalanches where he monitors those phenomena using acoustic emissions [45]. Stojanova et al. worked on fracture of paper using acoustic emissions created during crack propagation [46][47][48]. During the current experiments, acoustic signals are recorded using different sensors (shock accelerometers and piezoelectric sensors). Those signals are compared and investigated further in both time and frequency domains. Furthermore, during the experiments, photos of the Hele-Shaw cell are taken using a high speed camera. Thus, it is possible to visualize the complex branched patterns arising due to the solid-fluid interaction and to process images to gather information about the strain and strain rates, and investigate the mechanical properties of the solid partition.

Experimental Setup

A Hele-Shaw cell is made of two glass plates (80 × 40 cm) placed on top of the other, separated by 1.5 mm distance. The plates are separated via aluminum spacers placed close to the edges to provide equidistant spacing across the cell. For the experiments particular to this study, we completely sealed three boundaries of the Hele-Shaw cell and made one semi-permeable boundary using 50 µm steel mesh which allows fluid to exit the system but keeps the solid grains inside the cell. One of the plates has an inlet, which is used for injection of pressurized air, located at the bottom end of the plate, midway from the long edges 3 cm inside from one of the clamps close to the shorter edge (Figure 1A). The cell is filled with non-expanded polystyrene grains (80 µm diameter ±1%) called Ugelstad spheres (see details in Toussaint et al. [49]). The density of the spheres is 1.005 g/cm 3 . A mass of 170 g of grains is required to fill the cell, corresponding to an initial solid fraction of 52 ± 5% which is close to the 57% random loose pack for monodispersed grains [50]. Solid fraction of the grains in similar systems was investigated thoroughly by Johnsen et al. and found to be as low as 44% [23]. The difference between our experimental value and the theoretical value for an infinite box could be due to finite size of the cell which causes steric effect between the grains and the flat boundaries of the container, making the solid fraction 10-15% less [51]. Additionally, the measurement of 1.5 mm plate separation is subject to an error of 10% which has a direct effect on the solid fraction error bar. Electrostatic forces and humidity are also effective for real experimental beads and not for theoretical hard spheres, which can lower the solid fraction.

Some of the grains are colored with Indian ink to provide markers and texture, allowing a better resolution of the displacement measurements based on optical data. Images taken via high speed camera are used for digital image correlation to have full field measurements of displacement, velocity, strain rate etc. [21,[52][53][54]. The solid-air interface is placed 1-2 cm away from the inlet to avoid pressure localization close to the inlet. A schematic diagram of the Hele-Shaw cell is shown in Figure 1. The sides of the Hele-Shaw cell are clamped using steel clamps after sealing with double-sided rubber sealing tape. To protect the glass plates from stress focus on the clamping points, rubber sheets are placed between the clamp screws and the glass plates. Before placing the semi-permeable boundary, the cell is placed vertically and grains are poured inside. Following this, a semi-permeable filter is placed on the 4th edge. Then, the Hele-Shaw cell is positioned vertically, in a way that the semipermeable boundary stays at the bottom side to decompact grains and homogenize the solid fraction through the cell. Another important goal of this rotation process is to provide a small rectangular buffer empty of grains around the air injection inlet to avoid having point injection force over the medium. After the filling phase, the Hele-Shaw cell is placed horizontally.

Air injection is started and ceased via an electrovalve placed on the pipe very close to the air inlet, (Figure 2). This air pressure is provided via a pressurized air tank. Injection pressure over time is constant. It starts like a step function and is monitored using a pressure sensor placed on the air inlet allowing to check if the real injection pressure is within 5% of the value required during the experiment.

During the experiments, acoustic signals are recorded using different sensors. The data recorded on the piezoelectric sensors which are mostly sensitive in the range (200-900 kHz) are amplified with a Signal Preamplifier. The data recorded on the shock accelerometers which are mostly sensitive in the range (1 Hz-26 kHz) are amplified with a Brüel and Kjaer Nexus Charge Amplifier-Type 2692-A. The amplified/conditioned signal is transmitted to the computer via a Ni-DAQ mx PCI-6133 acquisition card with multiple channels at 1 MHz sampling rate (Figure 2). In addition, synchronized with the acoustical data, images of the Hele-Shaw cell are taken via a Photron SA5 high speed camera transferred and stored numerically. A TTL signal is used as a trigger to initiate the injection and the data acquisition via the camera and the acoustic sensors, thus, enabling time synchronization between the apparatus. Ambient lab noise is also recorded for reference and investigated using camera and accelerometer recordings prior to air injection. After recording, the signals are corrected by using the response function of the accelerometers provided by the manufacturer and cross-checked at the lab.

Experimental Observations

At large enough injection pressures, the fluid makes its way by creating channels and fractures toward the semi-permeable boundary as seen in Figure 3. In the beginning, the solid-air interface of the porous medium (closest to the injection point) moves more or less homogeneously, with the appearance of only large scale curvature of the interface (Figure 3B). Then after roughly 150 ms, some thin finger-like carved formations of thickness around 2 mm start to appear at several points (marked with yellow circles in Figure 3C). As injection continues, those fingers penetrate further in the medium. They get larger and wider with the help of the air pressure (Figure 3D). In addition to that enlargement, fingers branch out into thinner fingers. In the end, a tree-like branched channel network is created inside the porous medium. As a result of those fractures and channels, the surrounding material is displaced and the porous medium is compacted. Fracturing, channeling and fluid interaction inside the porous medium has its effects on the granular part of the medium. These interactions also result in granular transportation and compaction which involves inter-granular interactions as well as interactions of the solid grains with the confining glass plates. Initially the solid fraction is homogenous inside the plate. However, during the experiment the solid fraction of the grains increases in the close vicinity of the channels and fractures (up to the maximum possible solid fraction up to 63% for loose packed medium [50]. As observed in numerical models of such system in Niebling et al. [24], this solid fraction may decrease with distance from the fingers. After the experiment, depending on the pressure duration, initial solid fraction, and overpressure, some parts of the medium may have remained at the initial solid fraction, in other words may not have been compacted at all.

The acoustic events recorded during the experiments arise presumably due to sources of various types. While some vibrations are happening solely due to the air pressure fluctuation inside the channels, some others are generated by the stress increases on the plates due to granular compaction, intraporous air pressure vibrations, granular shocks and variation of the friction forces. These sources excite the confining plates, transporting mechanical waves to the sensors, i.e., the source types that are convoluted with the response of the Hele Shaw cell structure. Eventually, what is recorded is not just a signal created by a simple source, but a systemic response (i.e., signals interacting with plates and clamps, reflecting from edges, refracting through interfaces and eventually having different characteristic properties) to the many individual processes happening inside the plate during the whole period of the experiment. Signature of the signals recorded during experiments did not depend significantly of the sensor type and location.

Even though many individual acoustic events are superposed in time in the recorded signals, and are influenced by a systemic response, this does not mean that their specific signature is lost. Superposed signals may hide their signatures in the time domain, however their influence in the power spectrum may still be noticeable. In the following section, the evolution of the power spectral signature with varying solid-fluid interactions (e.g., compaction of the solid with fluid pressure, channeling, diffusion of the overpressure of the injected fluid through the pore spaces etc.) is shown. First, the power spectrum of several snapshots in time (i.e., Fast Fourier Transform, FFT), taken from different experimental stages, are presented. Then, they are analyzed and compared with each other. The flowchart in Figure 4 describes the analysis procedure.

Results

Power Spectral Evolution

The first time window analyzed, occurred prior to injection (i.e., at a state of rest). Figure 3 shows an image of the Hele-Shaw cell, start of the injection, channel initiation can be observed at several points, as highlighted by the yellow circles. Due to the interaction between the solid and fluid phases inside the Hele-Shaw cell, some mechanical signals are generated. The power spectrum of the signal recorded during this snapshot is presented in the bottom right panel of Figure 6. As channeling continues, the bulk movement of the grains takes place together with fluid motion. This causes emergent acoustic emissions which are predominantly in the low frequency range (less than 10 kHz). Vibration of the plates, granular friction and stress differences on the plates due to compaction create waves that are also in the power spectrum. All those contributions from different source mechanics gives a shape to the presented power spectra, similar to a power law decay having an exponent b = -1.61.

In addition to this slope, it is possible to find the mean frequency of the power spectra using Equation (1).

f = f |a(f )| 2 df |a(f )| 2 df (1)
This will show the dominant frequency range within the signal.

In the following Figure 7, it can be seen that the mean frequency starts very low and then increases with increasing energy in the high frequency range. As the channel network develops, we see that the mean frequency reaches to its maximum value.

Using cubic fitting, this mean frequency curve is estimated and compared with the optically obtained curves. An inflection point of the fitted curve is observed around time t 0 = 1.49 s after injection. This corresponds to the point where the finger ) is also presented.The area saturates to a maximum at a time close to t 0 . The length still increases slightly after t 0 , in particular, a small step in finger length (high slope of the length curve, corresponding to a large stick slip event connected to the finger) can be noticed around t 0 .

development stops and the regime inside the Hele-Shaw cell changes to a slow relaxation stage with slow fluid overpressure diffusion detailed in the discussion section.

Acoustic Events

It is possible to link the information received from the small scale lab experiments with large scale data, and vice-versa. In the experiments, we noted that the number of acoustic events occurring inside the Hele-Shaw cell is related to the empty channel area. However, after the fractured area reaches the final channel network shape, stress relaxation events are observed. These events are investigated further using event counting methods based on the ratio of Short Term Average over Long Term Average (STA/LTA) for event detection and compared with the evolution of the channel area inside the Hele-Shaw cell.

Event Detection

Detecting the number of events occurring during fracturing is a good indicator of the compaction level within the material. As long as there is motion of granular particles with the fast fluid flow it is very probable to detect some acoustic events in the experimental recording. However, it is very important to analyze and understand which part of the recording can be labeled an event and which part can not. STA/LTA threshold method is commonly used in seismic data interpretation [55][56][57][58]. If the parameters are selected carefully, it is very easy to use and is very robust [START_REF] Trnkoczy | Understanding & setting sta/lta trigger algorithm parameters for the k2[END_REF]. While LTA considers the average temporal noise to have an idea of the general behavior of the site, STA looks for intense changes in the signal in a small time window to detect acoustic events. Thus, it makes the ratio of those two parameters sensitive to the more complex events as well. When this STA/LTA ratio passes a pre-defined threshold, it is considered as an event (Figure 8), and the event counter is incremented. As long as the ratio stays above the threshold, it does not trigger again. Right after the ratio goes below the threshold, the algorithm can be triggered again for the next event to be counted. This algorithm (for one time window) can be generalized as follows:

(STA) (LTA) = 1 N s i+N s i (s(i) 2 ) 1 N l i+N l i (s(i) 2 ) (2) = (STA) (LTA) → > T → Event (3)
where s(i) 2 is the squared raw signal in the time domain (if necessary, a filtered signal s ′ (i) 2 can be used for different characteristic events), N s and N l are the length of the short (0.05 ms) and long time (1 ms) windows respectively and T is the predefined threshold for an event. The threshold to detect events may change between different datasets or different types of acoustic events.

Event Classification

One important thing that should be mentioned about the STA/LTA counting method is the frequency range used. After detailed analysis, it has become apparent that two different types of events exist within the experimental dataset. The first type of events are the non-impulsive low frequency (less than 10 kHz, similar to Figure 6) events (Type 1) which are related to the fluid flow and to the fluid-grain interactions rearranging the grains and producing major deformations and channeling. These events begin at the moment when the air injection starts and continue until a fully developed channel network is reached. It is possible to determine accurately this period from the optically acquired data, simply by calculating the area of the channel (estimated via the number of dark pixels in a binarized image) that saturates to a maximum value before the end of the acoustic emissions. Interestingly, in our findings, the number of Type 1 acoustic events follows a similar trend in time as the evolving emptied channel area within the Hele-Shaw cell (Figure 9). To enable the detection algorithm to distinguish between the different types of events, frequency filters are applied. To detect Type 1 events, a butterworth bandpass filter with corner frequencies 100 Hz and 10 kHz is applied to the raw signal s(i) before the STA/LTA event detection is applied.

The second type is the aftershock-like events, Type 2. Unlike the Type 1 events, in power spectrum their energy is spread over a wide frequency range similar to the one presented in Figure 10. These events are similar to the stick-slip relaxation events following a big earthquake in real scale. Compacted grains are rearranging their positions to have a more compacted state due to continuous air injection which results in some of their energy being released as acoustic emissions. To detect Type 2 events, a highpass filter for frequencies higher than 100 kHz is applied to the raw signal s(i) before the STA/LTA event detection is applied. In Figure 11 the evolution in the power spectrum with time is presented. In the figure, a Type 2 event occurred at t = 2.8 s is also presented. The raw signal is filtered by using a butterworth bandpass filter with corner frequencies 100 Hz and 10 kHz before computation. Bottom: The raw signal is filtered by using a highpass filter for frequencies higher than 100 kHz before computation. We found that some of the Type 2 events can be classified as Type 1 events since they have energy in the low frequency higher than the threshold. FIGURE 9 | Top: Signal recorded during the experiment. Bottom: Blue curve shows dark pixel counting over the image to compare compacted area with the number of different acoustic events occurred. Circles represent cumulative number of events normalized with maximum number of Type 1 events (i.e., 17 events). Cumulative number of Type 1 events (red points) are fitting well to the normalized channel area with time presented in the lower plot. First Type 2 (blue) event occurs very close to the inflection point of the mean frequency curve. Type 2 events are very impulsive and noticeable in the top panel as well.

After investigating further, we also noticed that the occurrence frequency of these events decay with time, similar to the Omori Law [START_REF] Omori | On after-shocks[END_REF][START_REF] Omori | On after-shocks of earthquakes[END_REF]. In the following section, curve fitting to the number of Type 2 acoustic emissions assuming an Omori Law decay (Figure 12) will be discussed.

Omori Law

Omori [START_REF] Omori | On after-shocks[END_REF][START_REF] Omori | On after-shocks of earthquakes[END_REF] worked on the half-day and monthly frequencies of aftershocks of the 1891 Nobi earthquake in Japan [START_REF] Omori | On after-shocks[END_REF][START_REF] Omori | On after-shocks of earthquakes[END_REF] . He found that the frequency of aftershocks n(t ′ ) at time t ′ can be expressed as:

n(t ′ ) = K(t ′ + c) -1 (4)
and the same equation for the cumulative number of aftershocks is given as

N(t ′ ) = t ′ 0 n(s)ds = Kln(t ′ /c + 1), ( 5
)
where c is a characteristic time, small and positive, K is the slope of the fit in the semi-logarithmic domain and N(t ′ ) is the number of cumulative aftershocks up to time t ′ . Following

Omori, [START_REF] Utsu | Magnitudes of earthquakes and occurrence of their aftershocks[END_REF] emphasized that the real aftershock activity decays with time differently than the originally derived Omori Law and proposed the equation for the frequency of aftershocks, with another fitting exponent p, as Utsu [START_REF] Utsu | Magnitudes of earthquakes and occurrence of their aftershocks[END_REF][START_REF] Utsu | The centenary of the Omori formula for a decay law of aftershock activity[END_REF]: 6) and called this the Modified Omori Law (MOL). The corresponding equation for the cumulative number of aftershocks is given as:

n(t ′ ) = K ′ (t ′ + c ′ ) -p (
N(t ′ ) = K ′ (p -1) c ′1-p -(t ′ + c ′ ) 1-p (7)
where c ′ is a characteristic time and K ′ is the slope having the dimension of time (p-1) . Using these approaches Omori Parameters K, c and p are estimated on our experimental catalog of Type 2 events, using time t ′ = tt 0 , where t 0 is the time defined in Section 4.1, that corresponds to the end of the channel formation and to the inflection point of the mean frequency recorded (Figure 12). A bin size to find frequency of occurrence of the events is selected as 0.2 s. However, it is more robust to use cumulative number of aftershocks (Equations 5 and 7) to avoid choosing an additional parameter for binning which can vary. Fitting is computed with 95% confidence bounds (i.e., two standard deviations away from the mean of the estimated value) to show the quality of fitting. Fitting parameters are calculated by a linear regression in the semi-logarithmic space of the cumulative number of events as K = 4.89 ± 1.33 and c = 0.13 ± 0.8 s. The Root Mean Square Error (RMSE) -average of the residual of the fit -is calculated to be 0.57 in this particular fit using the Equation (5).

Using the Modified Omori Law (i.e., Equation 7) fitting parameters are found as K ′ = 5.34 ± 1.52 s -0.45 , c ′ = 0.0066 ± 0.55 s, p = 0.55 ± 0.38 and RMSE is calculated to be 0.53 and the error bar is evaluated as the standard error [START_REF] Press | Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd Edn[END_REF]. The dashed lines in the Figure 12 represents two RMSE from the fitted curves.

Discussion

Discussion of the Experimental Results

The power spectrum of the mechanical signal at the different time windows show that the interactions inside the Hele-Shaw cell are evolving with continuous injection. In the beginning, there is a bulk movement of the grains due to the compaction caused by air injection. In Figure 6B it can be seen that the power spectrum follows a power law trend with an The circles represent the cumulative number of events and the solid lines represent the fitting curves given by the Equations ( 5) and (7). t ′ = 0 in this representation corresponds to time t 0 = 1.49 s after the start of the air injection.

exponent equal to -1.61 without any major peak. This indicates, although there are many different phenomena with various characteristic frequencies happening at the same time (e.g., bulk compaction, interactions between solid and fluid phase, stress redistribution among the medium etc.), that a scale-free mechanical phenomenon-probably related to the developing branching pattern-is dominating the emissions. Since the low frequency has higher energy, it can be seen that the mean frequency of the initial part t < 1.5 s is low as well (Figure 7).

As the channeling and fingering starts, the mean frequency is shifting from lower frequencies f < 10 kHz to the higher frequencies. Moreover, this evolution can be seen in Figure 11.

Rising mean frequency is also seen in the signals recorded during some volcanic processes, as e.g., the 2004 eruption process of the Arenal volcano, Costa Rica, and was interpreted as related to the stress increase with time due to the fluid pressure [START_REF] Almendros | Array analysis of the seismic wavefield of longperiod events and volcanic tremor at arenal volcano, Costa Rica[END_REF].

Additionally, the rate of this increase decreases as the channel network establishes. The trend which represents the channel area and also increasing permeability due to fracturing, is given in Figure 9 by checking the dark pixels in the images through time.

Occurrence frequency of low frequency earthquake-like events (described as Type 1) decreases with increased permeability due to fingering and fracturing indicated in the work of Frank et al. about real scale events seen in Mexico [START_REF] Frank | Along-fault pore-pressure evolution during a slow-slip event in Guerrero, Mexico[END_REF].

The inflection point (Figure 7 red dash dot line) of the mean frequency curve (t 0 = 1.49 s) shows the point where the power spectra evolves from a power law trend (Figure 6B) to stickslip aftershock like events (Figure 10, Bottom Right) whose power spectra consists of several peaks in the frequencies higher than 10 kHz (Figure 11) which causes the increase in the mean frequency. This inflection point occurs almost at the same time as the small step in the maximum finger length curve (high slope of the length curve, Figure 7 green) can be noticed around t 0 . This small step shows a stick slip motion in the most advanced finger tip. The inflection point also corresponds to the start of occurrence of the first Type 2 event which is presented in Figure 9. We can thus attribute it to a change in dynamics of the pneumatic fracture process: the main tree of channels stabilizes, and a slower stress relaxation process starts in the bulk around the channels, which leads to different type of acoustic event signatures, associated to the appearance of the impulsive Type 2 events.

The presence of this inflection point in the average frequency of the power spectrum, and the appearance of impulsive events (Type 2), can be suggested as the interesting signatures of the stabilization of the main channel during the monitoring of field scale hydrofracture or pneumatic fracture: these clear signatures can be detectable even in opaque field scale microseismic monitoring studies. This change of event types can also be seen in Figure 9 which shows the event detection using STA/LTA event detection procedure. As presented in the bottom panel of Figure 8, Type 2 events are detected starting from the inflection point of the mean frequency curve (t 0 = 1.49 s). Furthermore, we noticed that the cumulative number of Type 2 events are following an Omori Law trend which is shown in the Figure 12. This is due to the similarity between the aftershocks which is observed in real life and the Type 2 events in the medium which are due to the local stress relaxations. K = 4.89 ± 1.33 and c = 0.13 ± 0.8 s for the Omori Law and K ′ = 5.34 ± 1.52 s -0.45 , c ′ = 0.0066 ± 0.55 s, p = 0.55 ± 0.38 for the Modified Omori Law are computed as fitting parameters. It can be said that the Modified Omori Law-Equation ( 7)-gives slightly better fitting results to this dataset (RMSE is 0.57 for Omori Law and 0.53 for Modified Omori Law) due to the application of another parameter which is making the fitting procedure more sensitive.

The obtained p-value in this study (p = 0.55) corresponds to the p-value range for earthquakes studied by different scientists (i.e., 0.6-1.6 [START_REF] Utsu | The centenary of the Omori formula for a decay law of aftershock activity[END_REF], or 0.3-2.0 [START_REF] Helmstetter | Diffusion of epicenters of earthquake aftershocks, Omori's law, and generalized continuous-time random walk models[END_REF]). There are several studies for different types of earthquake mechanisms with varying p-values.

Utsu [START_REF] Utsu | The centenary of the Omori formula for a decay law of aftershock activity[END_REF] presented some of them: p = 0.9 ± 0.1 is found for the 1962 Westport earthquake in New Zealand, p = 1.4 is found for the 1965 Hindu Kush intermediate-depth (h = 214 km) earthquake, and p = 0.84 is found for the aftershocks following a rock fracture event in a mine causing a low magnitude earthquake (M = -2). This variability has been attributed to the variations in the state of stress, temperature, structural heterogeneities, material parameters, etc., in different tectonic regimes. However, a single, dominant factor controlling this parameter has not yet been identified [START_REF] Lindman | Physics of the Omori law: inferences from interevent time distributions and pore pressure diffusion modeling[END_REF].

In different examples, post-injection fluid-induced seismic events at various injection sites such as Fenton Hill, USA, and Soultz-sous-Forêts, France, appear to decay with an Omori-like fashion [START_REF] Langenbruch | Decay rate of fluid-induced seismicity after termination of reservoir stimulations[END_REF]. The rate of the post-injection seismic events in these cases has not been analyzed using the conventional MOL, thus a direct comparison between the obtained p-values cannot be carried out. In other examples, Nur and Booker [START_REF] Nur | Aftershocks caused by pore fluid flow?[END_REF] found that the decay of earthquakes following a step in the pore fluid pressure obeys a Modified Omori Law (MOL) with p = 0.5, which was consistent with a pore fluid diffusion law. Similarly, Yamashita worked on a model linking the fault slip and fluid flow in a system preloaded in shear. He found that the p-value is 0.48 for an early time period in a sequence of induced secondary aftershocks [START_REF] Yamashita | Regularity and complexity of aftershock occurrence due to mechanical interactions between fault slip and fluid flow[END_REF]. This link was also derived by and Rozhko et al. [75] for the events related to pressure changes in operation wells. Their formulation relates the rate of acoustic events to the temporal change in pore pressure. This is very similar to the p-value obtained by fitting the MOL to our experimental results.

Physical Explanation of the Experiments

The observed p-value in our experiments is also consistent with the derivations and approximations carried out in the models derived by Niebling et al. [20,24] or Johnsen et al. [25] for systems related to the one currently studied. The exponent of this Modified Omori Law and its prefactor can be directly related to the stress relaxation due to the diffusion of fluid pressure with the following approximations -that will be validated a posteriori by the agreement between the exponent and prefactor derived and the one that can be measured directly.

After the initial fast stage where the large empty channel is created, the compressibility of the fluid induces a slower seepage into the surrounding material. The boundary conditions can be considered as fixed for the granular material along the channel: owing to the large permeability of the channel, the boundary condition along the channel corresponds to an approximately homogeneous fluid pressure equal to the imposed inlet pressure, which also corresponds to the total stress since no grains are present in the channel, and the total pressure and pore pressure correspond on this boundary. The outlet boundary conditions can be considered as fixed pore pressure equal to the atmospheric one, and fixed displacement due to the semipermeable grid. The boundary conditions along the plates correspond to no flux for the fluid, and no normal displacement for the particles.

Since no large motion occurs in the following relaxation stage, the total stress field in the medium can be approximated as constant. This total stress σ T is the sum of the solid stress σ s arising from the forces transmitted in the grain-grain and grain-plate contacts, and the pore pressure p:

σ T = σ s + pId, ( 8
)
where Id is the identity matrix. This general formulation of summing the solid-bearing and fluid-bearing parts of the stress is for example demonstrated by Jackson [START_REF] Jackson | The Dynamics of Fluidized Particles[END_REF], and is valid as long as the large scale shear stress due to fluid-solid momentum exchange is negligible. It was shown to hold, for example, in sheared saturated granular layers [START_REF] Goren | A general criterion for liquefaction in granular layers with heterogeneous pore pressure[END_REF].

The total stress being fixed with fixed boundary conditions, the balance between solid-bearing and fluid-bearing stress changes as the fluid seeps in due to pore pressure relaxation. For every grain-grain contact and grain-plate contact, the stick/slip criterion can be described using a Coulomb friction criterion according to which slip occurs when the equality is reached:

σ s s ≥ µσ s n , (9)
where σ s n and σ s s are respectively the norms of the normal and shear (frictional) solid stress transmitted by the contact considered, and µ is the friction coefficient depending on the characteristics of the contacting elements (grains and plates).

The fluid can be explicitly considered by expressing this law in terms of total stress and fluid pressure, using Equations (8,9) which leads to a formulation corresponding to Terzaghi's 1936 effective stress formulation [START_REF] Goren | A general criterion for liquefaction in granular layers with heterogeneous pore pressure[END_REF][START_REF] Terzaghi | The shearing resistance of saturated soils and the angle between the planes of shear[END_REF]:

σ T s ≥ µ(σ T n -p), ( 10
)
where p is the pore fluid pressure. With fixed total stress, Equation ( 10) can be reformulated in terms of fluid pressure required for the grain pairs to slide:

p > σ T n - σ T s µ . ( 11
)
We assume a random distribution for the local total stress values, σ T s and σ T n , and call ξ the density of contacts per unit surface of the plates that will break with an acoustic emission due to sufficiently high pore fluid pressure. According to the above criterion, when p rises from the initial value to the final one a slip, during which micro acoustic events are produced, happens. The rate of events can be derived as follows: the system is considered to have a quasi one dimensional geometry where the grains ahead of the channel in the average flow direction correspond to x > 0. The initial value of the overpressure (i.e., pressure above the atmospheric pressure) is approximated as P = P max for x < 0 and P = 0 for x > 0. This holds if the initial channel creation was fast compared to pressure diffusion, so that the pressure skin depth stayed small with respect to the system size during channel creation, which was verified numerically by Niebling et al. [20].

The variation of the pore pressure over time is thus the key factor in this system to change stability of the porous medium. This pore pressure can be expressed using the fluid pressure diffusion into the porous medium. Considering that the grains are not moving, the pore velocity v d of the fluid (local velocity in the pores) and the Darcy velocity v D , can be computed from Darcy's law. As presented e.g., in Niebling's work [22] or Johnsen's [26]:

v d = v D φ = - κ φµ ▽P ( 12
)
where κ is the local permeability, φ is the local porosity of the granular medium, µ is the fluid viscosity and P is the pressure gradient driving the flow. Implementing this equation into conservation of mass we have:

∂ t (φρ(P)) + ▽(φρ(P)v d ) = 0 ( 13
)
where ρ(P) is the updated mass density of air, obtained assuming the state equation of a perfect gas which is following the relation ρ(P) = ρ 0 P P 0 which is valid for this type of setups since the pressure is not varying by orders of magnitude with no strong density changes over the cell. When the pressure dependence of the density is included in the analysis which is indeed more correct, the diffusion equation becomes nonlinear. It can be treated numerically, Niebling et al. [20]. As shown by Niebling et al., no qualitative changes are noted in this kind of regime though and the same kind of scaling law for the growth of the skin depth is observed. In this equation, ρ 0 is the density of air and P 0 is the atmospheric pressure. This puts Equation ( 13) into the form:

ρ 0 P 0 ∂ t (φP) -▽(ρ 0 κ µ ▽P) = 0. ( 14
)
Considering that porosity (φ = φ 0 ) and local permeability (κ = κ 0 ) does not change with time, after some simplification Equation ( 14) will become:

φ 0 P 0 ∂ t (P) - κ 0 µ ▽ 2 P = 0. ( 15
)
Then, from Equation ( 15) we can derive the diffusion equation as:

∂ t (P) = κ 0 P 0 φ 0 µ ▽ 2 P ( 16
)
where the part before the pressure gradient of Equation ( 16) can be described as diffusivity constant of overpressure in a porous medium, corresponding to the slow Biot wave [START_REF] Masson | Finite difference modeling of Biot's poroelastic equations at seismic frequencies[END_REF]:

D = κ 0 P 0 φ 0 µ . ( 17
)
The local permeability κ 0 can be computed using Carman-Kozeny equation [22]:

κ 0 = d 2 180 φ 3 (1 -φ) 2 (18)
where d = 80µm is the diameter of the grains and φ, the porosity of the medium, is around 48% which is typical for this kind of preparation. The results of Niebling et al. [20] show that in this kind of system, fast forming aerofractures are formed, followed by a slow diffusion of overpressure away from the large channelsour slow relaxation stage. He modeled the behavior of the aerofractures with two different stages, with initially a fast channel formation in the beginning due to the high fluid pressure and following, close to the injection point, a thick compaction front which satisfies Equation (11).

As the overpressure diffuses in the medium, the grains can slip and rearrange, possibly giving rise to the acoustic events. The fundamental solution of Equation ( 16) is derived with an approximated 1D boundary condition, having an initial condition P = 0 at t = t 1 describing the finger tip as a flat boundary at x = 0, and the space ahead of it as x > 0. Then, the initial condition for the overpressure is P = 0 for x > 0 and t = t 1 (when the finger stops and the overpressure has not yet penetrated far from the boundary). Furthermore, boundary conditions corresponding to the imposed overpressure are P = P max at x = 0 (the finger boundary), and P = 0 at x → +∞ (far ahead of the finger). The solution of the diffusion Equation (16) gives the overpressure at distance x ahead of the tip of the empty channel, and at time t > t1, after the end of the channel formation [25]:

P(x, t) = P max erfc( x 2 √ D(t -t 1 ) ) = P max 2 √ π ∞ x/2 √ D(t-t 1 )
e -y 2 dy (

) 19
where t 1 is the time when the end of channel growth stops, P max is the saturation level of pressure. As it is described in Equation ( 11) fluid pressure is directly effective on the number of failing contacts. This exponentially decreasing pressure field, for order of magnitude estimates, can be approximated as a zone of (skin) depth 2 √ D(tt 1 ) ahead of the finger tip experimenting a significant overpressure rise close to the maximum overpressure P max , and a negligible overpressure in the zone ahead of it. Since ξ contacts break-with an acoustic emission-per unit area when the overpressure rises from 0 to P max , the cumulative number of contacts failing (i.e., the cumulative number of acoustic emissions) can be approximated via the equation below:

N(t) = 2Lξ D(t -t 1 ) = (t -t 1 )/τ ( 20
)
where L is the width of the cell. Furthermore, there is a characteristic time related to the pore pressure diffusion and density of contacts breaking due to this pore pressure increase,

τ = 1/(4DL 2 ξ 2 ). ( 21
)
This time τ is different from the other characteristic time t 1 , that corresponds to the start of the diffusive behavior for the pore pressure.

Equating Equations ( 7) and ( 20) (and recalling that the time t ′ is measured with respect to t 0 in this last equation), we obtain the prediction:

((t -t 1 )/τ ) 0.5 = K ′ 1 -p ((t ′ + c ′ ) (1-p) -c ′(1-p) ), (22)
which predicts that 1 -p = 0.5, i.e., p = 0.5, t -

t 1 = t ′ + c ′ , c ′ = 0 which leads to t -t 1 = t ′ . Considering t ′ = t -t 0 ,
this becomes t 1 = t 0 and finally Equation ( 22) turns into 1-p) . This is indeed the case up to the error bars: t 1 = t 0 -c ′ , with c ′ = 0.0066 ± 0.55 s, i.e., the start of the diffusive regime t = t 1 corresponds to the end of the growth of the finger t = t 0 up to the error bar, and p is found to be p = 0.55 ± 0.38, equal to 0.5 up to the error bar. The prefactor of the MOL, found with the central p-value of 0.55 to be K ′ = 5.34 ± 1.52 s -0.45 , allows to evaluate the characteristic diffusion time τ = (K ′ /(1 -p)) 1/(p-1) = (5.34/0.45) (-1/0.45) = 0.004 s. (23) The frequency of occurrence of events is the time derivative:

K ′ /(1 -p) = τ -(
dN dt = 1 2 √ τ (t -t 1 ) . ( 24
)
The present study thus suggests that p = 0.5 in the Modified Omori Law is the signature of the slow stress relaxation due to the diffusion of the overpressure in the medium surrounding the cavities.Presumably, this p-value may decrease if the fluid injected slowly enough so that the diffusion skin depth becomes large with respect to the channel width during the injection, which is a contrary to the situation depicted here [20]. Also, large injection pressure may lead to subcritical crack growth and change p-value [40,42,[START_REF] Santucci | Slow crack growth: models and experiments[END_REF].

The prefactor of this law, corresponding to the characteristic time τ = 0.004 s for this diffusion, can thus be related to the porosity and permeability of the medium around the channels, the viscosity and compressibility of the fluid, and the density of failing sites leading to events.

Putting τ back into the Equation ( 21) it is possible to calculate ξ via the following equation:

ξ = 1 2L √ Dτ (25)
Recalling the expressions Equations ( 17) and ( 18), we obtain for this medium a diffusivity

D = d 2 φ 2 P 0 180(1 -φ) 2 µ = 0.20 m 2 /s (26)
where we have used the air viscosity, µ = 1.810 -5 Pa s, and the atmospheric pressure value P 0 = 10 5 Pa. These two last equations, Equations ( 25) and ( 26), allow to express the density of the triggered seismogenic contacts, with L = 0.4 m, as:

ξ = 1 2L √ Dτ ≃ 45/m 2 (27)
for this presented experimental setup. This indicates that around 45 contacts per square meter can give rise to strong events due to overpressure rise (which is far lower than the total number of contacts ≈ 2.3 × 10 10 for the grains and cell thickness considered). This shows that just a small subfraction of contacts give rise to strong events. It is also very consistent within the order of magnitude with the fact that 12 events of Type 2 were observed during this relaxation in the system of size around 50 cm ahead of the main finger by 40 cm width, i.e., of size 0.2 m 2 corresponding to a density around 50 m -2 microseismogenic events happening during this stress relaxation.

Conclusion

A purpose-built Hele-Shaw cell experiment was designed to enable both optical and acoustic recordings associated with controlled fracturing of a porous medium via air injection. The optical and acoustic recordings are analyzed together to obtain a further understanding of the type of deformation occurring within the cell. Based on the evolution of the power spectrum and mean frequency of the acoustic data with time, together with the growing channel network, it can be determined that the low frequency content f < 10 kHz of the acoustic emissions is directly related to the permeability state of the medium. Low frequencies dominate the power spectrum as long as the medium has not reached its final fractured state. After this state is reached, aftershock-like events seem to release the stress from the medium in the post-fractured phase. These aftershock-like events have a broadband spectrum where the energy is not focused on the low frequency but spread on a wide span of frequencies. These events occur right after reaching the final fractured state and their frequency of occurrence decays with time. Both power spectra evolution and diminishing frequency of occurrence are present in real scale microseismic data. In addition, it is possible to estimate the number of aftershocks by using Modified Omori Law in experimental and microseismic data. These signatures, inflection point in the average frequency, appearance of impulsive events of high frequency, and starting point of an Omori Law, can be used straightforwardly in large scale microseismic monitoring of fluid injection and well stimulation. The permeability of the medium can also be directly estimated from the prefactor of the Omori Law.

74

Chapter 3. Investigating Solid-Fluid Interactions inside the Hele-Shaw Cell

Draft Article: Numerical Studies of the Acoustic Emissions during Pneumatic Fracturing

In the previous section, we analyzed the experimental data of the system acquired using acoustic sensors and a fast speed camera. Even though these mentioned tools gave us a lot of information about the system there are still some questions that are not answered by these analyses (e.g. air vibrations inside the system, contributions to the different peaks of the spectrum etc.). It is possible to address some of these questions by using numerical approaches. Here in this section we explain how to have more information about these unknowns using a numerical model, with a physics based on the experimental design and observations. This study is mainly conducted by the main author under the supervision of Renaud Toussaint. Then, to develop the work further, we discussed several times with all coauthors. All co-authors have contributed with some suggestions on references, scientific methods on analysis and comparison, and presentation of the obtained results. This draft article is going to be submitted to Journal of Geophysical Research: Solid Earth in early 2017. Please consult the latest version.

To establish the numerical model, first we split the total stress into two contributions: (1) fluid pressure, (2) solid stress. (1) To simulate the fluid pressure field evolution, we used the conservation of mass equations for the solid and the fluid in the Hele-Shaw cell. Afterwards, we introduced momentum exchange between solid and fluid parts using Navier Stokes equation. Due to slow fluid flow in pores, Reynold's Number is smaller than 1. (2) For the solid stress at play between the grains and the surrounding plates, we used a Coulombian friction model, together with Jannsen's equation modified for horizontal setups. After combining these two stress maps, we obtained total stress evolution map inside the Hele-Shaw cell. Using how the stresses are evolving inside the cell, it is possible to calculate the vibrations. We convolved these stress variations with the wave propagators inside the glass plates constituting the Hele-Shaw cell sides, i.e. the Lamb wave Green function as derived by GoyderWhite. Then we compared the power spectrum of these synthetic vibrations with experimentally recorded ones.

After checking different contributions on the power spectrum, we found that the pore pressure variations due to air vibrations in the open channels are essentially governing the power spectrum. In the beginning (t < 0.05 s) air pressure bursts into the empty area in the cell, and eventually reaches and pushes the granular medium, which creates the major contribution in the low frequency region (f < 20 kHz). We found that the major peak of the power spectrum on experimental and numerical results are corresponding very well. Furthermore, we noticed a power law decay in the power spectrum of air vibrations (as mentioned in experimental results). This might be the influence of the channel network which displays a fractal shape, having no characteristic size in space.

Moreover, as the grains get compacted around the main channels, the solid stress starts to generate powerful signals. Since the solid stress has a step-like loading due to sudden compaction and sudden stress release (if carving occurs after compaction), this loading is exciting all types of frequencies on the glass plates of the Hele-Shaw cell, which eventually generates a power spectrum similar to the experimentally observed Type-2 events (as explained in the previous section). Comparing directly the simulations with the experimental results, we obtain that the rate of occurrence of these solid stress based events are similar to the experimentally observed Type-2 events.

Introduction

Fluid flow inside granular media is a very common phenomenon in nature and in industrial applications. Particularly if the fluid flow is pressurized, it can induce rapid porosity changes inside the porous medium which is in contact with the fluid [Gidaspow, 1994;Kunii and Levenspiel, 1991;Goren et al., 2010Goren et al., , 2011]]. This type of brittle deformation is visible in nature (e.g. volcanic activities, tremors), and in engineering (ground improvements, CO 2 sequestration, fracturing applications) [Charléty et al., 2007;[START_REF] Cuenot | Analysis of the microseismicity induced by fluid injections at the EGS site of soultz-sous-forłts (Alsace, France): implications for the characterization of the geothermal reservoir properties[END_REF]Dorbath et al., 2009;Aochi et al., 2011;Schuring et al., 1996;[START_REF] Agency | Accutech Pneumatic Fracturing Extraction and Hot Gas Injection, Phase One: Applications Analysis Report[END_REF]Gao et al., 2014]. Deformations based on momentum exchange between non-miscible solid and fluid phases can lead to different characteristic shapes and structures of deformation (bubbles, fractures, channels). These shapes change depending on the characteristics of the solid (density, particle size, permeability, etc.) and the fluid (viscosity, compressibility, etc.) and boundary conditions (applied pressure, shape and permeability of the boundaries, gravity orientation, etc.). The microseismicity monitoring is a routine protocol during well operations in the industry, but interpretation of this data requires further

-2-
research [START_REF] Valkó | Hydraulic Fracture Mechanics[END_REF]Cornet et al., 1998;Cornet, 2015]. In the experimental work [Turkaya et al., 2015], authors recorded the vibrations on the glass plates of a Hele-Shaw cell. These vibrations are generated during the changes in the total stress in the Hele-Shaw cell which are due to the variations of the applied load on the glass plate.

Investigating the experimental setup, it is possible to say that these changes in the applied load are caused by three major things [Turkaya et al., 2015;Niebling et al., 2012a,b]: Air

Vibrations Inside the Carved Empty Area, Pore Pressure Diffusion Into The Pores, and Solid Stress due to Compaction. We will here describe the basis of a numerical model taking these effects into account. First, we will explain the small simple regimes which we observe inside the Hele-Shaw cell and then combine these mentioned effects to obtain a more general poro-elasto-plastic formulation. Next, we will study the spectral characteristics of the waves emitted during the process, generated by this code, and compare these to the experimentally recorded ones. We will show that how the corner frequency in the power spectrum transforms to a power law decay as the channel network evolves, what are the contributions of the different stress sources on the peaks of the power spectrum and their characteristic properties (frequency band, corner frequency, slope etc.).

Air Vibrations Inside the Carved Empty Area

Air Vibration in an Empty Room can be explained and modeled using the classical wave equation for sound waves [START_REF] Royer | Elastic Waves in Solids I: Free and Guided Propagation[END_REF][START_REF] Griffiths | Linear and nonlinear waves[END_REF]:

∇ 2 P = 1 V 2 s ∂ 2 P ∂t 2 (1)
where P is the air pressure and V s is the sound wave velocity in the air. The same equation is valid for chambers with irregular shapes as well such as carved channels. For the geometry in which the air vibrates, photos taken during the experiment are used. First, at every time step, the image that is going to be used is binarized and rescaled to make the simulation computationally cheaper and stable. We could apply this equation to the determined empty space, with impermeable boundary conditions. In this study we are interested in, the air -porous medium interface are not sealed but permeable however, this equation is covering only the cases where all the boundaries are sealed (no leakage of the fluid into the pores). Thus, the pressure diffusion inside the porous medium needs to be included to make the simulation realistic.

-3-

Air Pressure Diffusion into the Granular Medium

As the air is pumped into the medium, it vibrates in the open space. When the overpressure reaches to the solid-fluid interface, as fluid penetrates into the porous medium, the stress state changes both in the fluid and the solid. If we assume that the solid part is permeable and rigid enough to not to be deformed by fluid stress, (which will be eventually included in chapter 2) we can use the classical diffusion equation to model this diffusion. The classical pore pressure diffusion equation is given as [Johnsen et al., 2008;Turkaya et al., 2015]:

∇ 2 P = 1 D ∂P ∂t , ( 2
)
where D is the diffusion constant. A model taking into account both regimes in a single framework, amounting to a model of poro-plasto-mechanics, will be developed in the next sections. We will show that in the two limits of large and low porosities, it is equivalent to the two behaviors described here (i.e. the diffusion one and the vibration one).

Granular Compaction

During fracturing, the state of stress in the solid part changes as well. In a Hele-Shaw cell like this where the solid medium cannot go outside the boundaries, compaction in the porous medium is the only way to create empty, carved areas inside the medium.

Starting from the solid-fluid interface, there is a region where the grains are fully compacted (compaction front). As presented in earlier works [START_REF] Knudsen | Granular labyrinth structures in confined geometries[END_REF]Sandnes et al., 2011] L = φR 1-φ is the thickness of this compaction front (see Fig. 1) with typically L = 3 cm that depends on the solid fraction φ = 0.48 and the radius of curvature of the fingers R = 6 cm which are obtained from experiments for this simulation:

σ(x) = ρ s g∆z 2ψ {ψγ + 1} e 2γψ x ∆z -1 , (3)
here in this equation ρ s is the initial density of the granular medium, x is the distance from the solid-fluid interface (see Fig. 1), ∆z = 1.5 mm is the aperture of the Hele-Shaw cell, ψ is the Janssen's coefficient [START_REF] Janssen | Versuche über getreidedruck in silozellen[END_REF][START_REF] Knudsen | Granular labyrinth structures in confined geometries[END_REF]Sandnes et al., 2011] taken here as 1, γ is the friction coefficient of the porous medium taken as 0.3.

Solving the equation 3, it is possible to determine the solid stress map of the medium at a specific time t in the cell.

-4- -5-

As explained in Figure 2, the stress variations on the glass plate generate Lamb Waves on the Hele-Shaw cell. Implementing this stress map to the far field Lamb wave approximation presented in [START_REF] Goyder | Vibrational power flow from machines into built-up structures, part i: Introduction and approximate analyses of beam and plate-like foundations[END_REF] we obtain the vibrations on the plate.

In this study, as presented in the experimental work, we expect the power spectrum and thus the mean frequency of this power spectrum to evolve as the medium is fractured.

The experimental and numerical results in various aspects will be compared in section 5.

Numerical Modeling of the Air Pressure

We start with the Navier-Stokes equation [START_REF] Guyon | Hydrodynamique Physique. Savoirs Actuels[END_REF] adapted to the Darcy flow [START_REF] Darcy | Les fontaines publiques de la ville de Dijon: exposition et application[END_REF] as given in [Niebling et al., 2010a;Jackson, 2000;Van der Hoef et al., 2006] as:

ρ f ∂ ∂t (v f ) = -∇P - v f -u κ µφ, (4)
where ρ f is the density of the injected air, φ is the porosity, v f is the velocity of the air, u is the granular velocity, µ is the viscosity of the air, κ is the permeability of the medium computed using the Carman-Kozeny equation for porous medium K por e = 3.5e -11 m 2 [Turkaya et al., 2015] and for the carved parts assuming a Poiseuille flow between two plates K empty = 1.59e -7 m 2 [Niebling et al., 2010a]. From equation 4 it is possible to express v f as:

v f = u - κ φµ ∇P + ρ f ∂(v f ) ∂t . (5)
During this coupled flow, the mass (both air and grains) is preserved. Thus, the stresses inside the Hele-Shaw cell can be explained with the mass conservation equations of the air and solid phases [Johnsen et al., 2008;Niebling et al., 2010a,b]. For the air phase mass conservation can be written as:

∂ ∂t (ρ f φ) + ∇ • (ρ f φv f ) = 0, (6)
Implementing v f in Eq. 6 we have:

∂ ∂t (ρ f φ) + ∇ • ρ f φ[u - κ φµ ∇P] - ρ 2 f κ µ ∂v f ∂t = 0, (7)
Taking into account that that the air is compressible, the change in the air mass density can be expressed as:

ρ f = ρ 0 + ρ 0 c(P -P 0 ) (8)
-6-where ρ 0 is the reference density of air at atmospheric pressure and room temperature, and c = 1/P 0 is the isentropic compressibility of the air -considering air as a perfect gas.

Then, for the grains the mass conservation equation is defined as:

∂ ∂t (1 -φ) + ∇ • ((1 -φ)u) = 0. (9)
which can be simplified into ∂φ ∂t = -u • ∇φ + (1φ)∇ • u. Inside the pores, introducing the compressibility of the air (eq. 8), the eq. 7 will transform to:

φρ 0 ∂P ∂t + ρ f ∂φ ∂t = ∇ • [ρ f κ µ ∇P] -∇ • (φρ f u) + ∇ • ( ρ 2 f κ µ ∂v f ∂t ), (10)
If the grains are moving with fluid, Darcy contribution in eq. 7 will vanish. Then, we can get:

ρ f ∂φ ∂t + ∇ • (φρ f u) = ρ f ∇ • [(1 -φ)u] + ∇ • [φρ f u], (11)
which simplifies to:

ρ f ∂φ ∂t + ∇ • (φρ f u) = ρ f ∇ • u + uφ∇ρ f ( 12
)
introducing the compressibility equation (Eq. 8) to the end of the right hand side (i.e.

∇ρ f ) we have:

ρ f ∂φ ∂t + ∇ • (φρ f u) = ρ f ∇ • u + φρ 0 cu • ∇P. (13)
Putting this into eq. 10 we get: 14) which is equal to:

φρ 0 c[ ∂P ∂t + u • ∇P] = ∇ • [ρ f κ µ ∇P] +∇ • [ ρ 2 f κ µ ∂v f ∂t ] -ρ f ∇ • u (
φ[ ∂P ∂t + u • ∇P] = ∇ • [ P κ µ ∇P] -P∇ • u +∇ • [ ρ 2 f κ ρ 0 cµ ∂v f ∂t ] (15)
where P = ρ f ρ 0 c = P -P 0 + 1c. It should be noted that inside the pores, the pressure difference is small compared to the background pressure, which leads to P ≈ P 0 [Niebling et al., 2010b]. It is possible to apply some simplifications to this equation. In a clear fluid we can say that ∇ •

ρ f κ ρ 0 cµ ∂v f ∂t reduces to ρ f κ ρ 0 cµ ∂∇•v f
∂t . Furthermore, for eq. 6 is simplified using three assumptions:

(i) ∇ • (ρ f φv f ) ∼ ρ f φ∇ • (v f )
) ∂∇•v f ∂t ∼ -cρ 0 /ρ f ∂ 2 P ∂t 2 ∼ -c ∂ 2 P
∂t 2 this approximation is valid inside the porous medium, obtained via implementing the Darcy velocity v f = κ φµ inside the diffusion equation in the porous medium (i.e. eq. 2).

Applying all these simplifications we eventually have the general formulation that accounts for the air vibrations in the empty area and the pressure diffusion in the porous medium:

[ ∂P ∂t + u • ∇P] = ∇ • [ P κ µφ ∇P] -P∇ • u - ρ 2 f κ ρ 0 µφ ∂ 2 P ∂t 2 . ( 16
)
When κ is large (i.e. in the empty channels) we have P∇ 2 Pρ f

∂ 2 P ∂t 2 = 0. Knowing that, P ρ f = 1
V s , we can obtain the air vibration equation in the empty channels (eq. 2) we defined in the introduction. Moreover, if κ is small (i.e. in the pores) last term vanishes and the eq. 16 converges to eq. 2.

To compute u we subtracted the experimental images that are taken with 8 ms timesteps.

The difference on the interface of the porous medium is computed as distance. To avoid any numerical instabilities, granular velocity in the empty area is set to zero. The velocity of the grains are computed on the interface and decayed proportional to the compaction zone.

Furthermore, by adding the solid stress computed via Eq. 3 one can compute the force map on the glass plates. As described in figure 2 first, the stress map on the plate in an instant in time is computed. Then, the signal is generated by convolving force variation on point with a green function. Following this, these signals are compared in time and in Fourier domain to see the evolution due to solid-fluid interactions inside the Hele-Shaw Cell.

Numerical Solution Scheme

We used a Forward Finite Difference Approximation Explicit Formulation (i.e. f ′ (x) ≈

f (x+h)-f (x) h
where f (x) is the function, h is the step along x [Levy and Equations, 1992]) is applied to solve numerically Eq. 16 inside the Hele-Shaw cell during air injection.

Since the equation 16 is made of combining several equations, forward difference method is the simplest (even though presence of the stability problems) to implement to this case.

-8- Defining parts of this equation we have:

∇P = P x+1,y i -P x,y i ∆x + P x,y+1 i -P x,y i ∆y , (17)
and

∇ 2 P = P x+1,y i -2P
x,y i + P

x-1,y i ∆x 2

+ P

x,y+1 i -2P

x,y i

+ P x,y-1 i ∆y 2 , (18)
can be given for the gradients in space. For the derivation in time we discretize as:

∂P ∂t = P x,y i+1 -P x,y i ∆t , (19)
and

∂ 2 P ∂t 2 = P x,y i+1 -2P x,y i + P x,y i-1 ∆t 2 . ( 20
)
Additionally, we define the diffusivity constant D = κ µcφ and wave velocity V s = 1 ρc . Putting equations 17, 18, 19, 20 back to Eq. 16 we will have:

P x,y i+1 -P x,y i ∆t + ρ 2 f κ ρ 0 µφ P x,y i+1 -2P x,y i + P x,y i-1 ∆t 2 = P κ µφ P x+1,y i -2P
x,y i

+ P x-1,y i ∆x 2 + P x,y+1 i -2P
x,y i

+ P x,y-1 i ∆y 2 - P φ u x+1,y i -u x,y i ∆x + u x,y+1 i -u x,y i ∆y - u φ P x+1,y i -P x,y i ∆x + P x,y+1 i -P x,y i ∆y (21)
-9-and simplifying eq. 21 we have the final equation:

P x,y i+1 = 1 ∆t+ ρ 2 f κ ρ 0 µφ ∆tP x,y i + ρ 2 f κ ρ 0 µφ 2P x,y i -P x,y i-1 + Γ (22)
where Γ is the right hand size of the eq. 21. Solving equation 22, we obtained the map of the fluid pressure evolution inside the Hele-Shaw cell. Forward difference scheme is very sensitive to the ratio of the time step to the grid size in space ∆t/∆x = 0.0013 s/m. The result of this equation with the corresponding image of the Hele-Shaw cell is given in figure 3 at t = 0.16 s starting from the time of injection. The corresponding colorbar shows the pressure applied on the plate in Pa. The pressure is applied at the right boundary from time zero as a line source and kept constant through time. For the cases with noise we sent the pressure with a Gaussian noise over time to reflect the pump fluctuations. The boundaries at the top and the bottom are sealed and the boundary on the left is fixed to atmospheric pressure (i.e. P = 0). As indicated in figure 2, having the pressure variation on the glass over time, it is possible to generate vibrations synthetically by convolving the forces with the Green function of Lamb waves, as will be explained in the next section.

Signal Generation

A far field approximation of the Lamb waves is given in [START_REF] Goyder | Vibrational power flow from machines into built-up structures, part i: Introduction and approximate analyses of beam and plate-like foundations[END_REF]Farin et al., 2016;Turkaya et al., 2016] as:

ξ(R) = -iF(k) 8Bk 2 2 Rkπ e -i(Rk-π/4) (23)
where ξ(R) is the displacement of the plate in Fourier domain in a distance R from the source, F(k) is the applied force in Fourier domain at a point like source, B is the bend-

ing stiffness B = h 3 E 12(1-ν 2 )
where h is the thickness of the plate, E is the Young's modulus, and ν is the Poisson's ratio of the plate material. k is the wavenumber where dispersion relation with angular frequency

ω is k = ω 1/2 ρh B 1/4
. If the frequency is higher than the tions simple, in the numerical studies, we did not include the reflections. Even though this may cause some differences in the amplitude on the power spectral distribution curves, the main structure which is the essential part of this study remains the same.

Results and Discussion

Temporal Evolution of Pressure inside the Cell

The pressure at a given point in a Hele-Shaw cell evolves with the evolution of the medium during injection process. Several mechanisms (location of the point selected, the material in this discrete volume, distance from the inlet etc.) are effective to change the pressure exerting on the glass from that finite volume. In figure 4, the evolution of stress at a point is shown. Two large jumps, due to the compaction zone of the solid passing through the point, are clearly visible. After the first carving (around t = 0.2 -0.3 s), again a second time migrating solid partition passes from the point, causing a sudden increase in the stress applied on the glass at t ≈ 0.3 s.

Since the only vibration source is the total stress variations inside the cell, the frequency of these changes are directly linked to the power spectrum of the Lamb waves generated on the glass plate.

In figure 5, we see that there is an increase in pressure at all points in varying slopes based on their distance to the outlet. In the second phase, as the solid-fluid interface reaches to the marked point on the bottom image (in red (triangle up), blue (circle), green (triangle down), and cyan (square) curves), pressure oscillations due to the rapid momentum exchange occurring between solid and fluid parts are observed. These vibrations are due to fact that the solid part is deforming as the air pressure penetrates into the medium, carving and compacting the granular medium around it.

Mean Frequency Evolution of Pressure Over Time

To compare the evolution of pressure inside the Hele-Shaw cell, we considered small time windows (5 ms) of the Lamb waves generated with different vibration sources all over the plate. After applying a Fast Fourier Transform to these windows, the mean frequencies are computed using the equation 24:

-

f n = ∫ f s( f , t n )df ∫ s( f , t n )df , (24)
where f n is the mean frequency of the time windows n and s( f , t n ) is the power spectrum of that specific time windows -similarly to the definition of Vinningland et al.

[ Vinningland et al., 2007aVinningland et al., ,b, 2010]. In figure 6, we compare the mean frequency evolution of the experiment with the different contributions of pressure in the Lamb wave simulations. Initially, the experimental mean frequency (top subfigure, blue curve) is following the simulations based on total stress. Then, the experimental average frequency decreases to the simulated air pressure average frequency (when % 50 of the channeling is completed t ≈ 0.21 s) and follows it until the end of experiment except the instants when the stick slip events are occurring (e.g. around t ≈ 0.9 s). In these cases, the mean frequency jumps up to the total stress -15-level and decreases back to the air pressure one in the end of the event. This experimental mean frequency is accurate as long as frequencies essentially under 50 kHz are involved.

To ensure this accuracy, a lowpass filter of 50 kHz is applied with sensor sensitivity correction.

Since the simulation results are not showing a mean frequency rise in the experimental results after channel formation, we tested these results without a lowpass filter to compare with the results obtained in [Turkaya et al., 2015] to see the mentioned characteristic properties of the mean frequency in simulation results. After recomputing the mean frequencies without filter we saw that the mean frequency rise due to stick slip events after the channel formation exists as well. In addition to the results mentioned in that article, it is possible to say that due to higher pressure (P ′ = 0.5 bar in the article, P ′ = 1.0 bar in this research) used in this experimental results, there are more oscillations in the mean frequency towards the end of the experiment.

Power Spectral Evolution of Pressure Over Time

The vibrations on the plate are due to the stress variations on the glass plates. It is possible to investigate the different signatures on the power spectrum by looking at a specific point to associate the change in the spectrum to a type of change in the stress.

In figure 7 (Top), spectrogram due to vibrations generated by air pressure in empty region is shown. (Middle) Spectrogram due to vibrations generated by total stress in the cell is shown. Colorbar of spectrograms are plot in logscale (log 10 ). (Bottom) Air pressure and total stress evolution over time on a single point is shown. It can be seen that initially with increasing pressure the power spectrum becomes more energetic. As it can be seen in the bottom figure, the only difference in these two sources are when the compaction zone passes through the point taken as the source. The spike in the loading corresponds to the compaction front passing on the point where the power spectrum has an abrupt increase in all frequencies as well. Following this, the point remains in the carved area where the air is bursting into the empty channel giving power to the low frequency range of the spectrum. As the air penetrates into the pores towards the boundary open to, the pressure and the vibrations diminish, leading to a decay in power spectrum.

Another parameter to study is the effect of the noise in the input pressure. In figure 8, two spectrogram belonging to the cases having flat and noisy input pressure are com- pared. Gaussian noise is added to the input pressure (1 dB noise is added to a flat unit overpressure before multiplying with the overpressure P ′ = 0.5). Both spectrograms have the same colorbar in (C * m 4 /s 2 ), corresponding to (in log 10 scale). In the top figure, the rise and the decay over time in the spectrogram is clearly visible, but in the bottom figure where the air input is noisy, this noise covers the decay in the power spectrum. However, the behavior of the general system does not significantly change with respect to the noisy inlet pressure.

It is possible to investigate further the evolution in power spectrum by looking to the entire spectrum.

In figure 9 the experimental and numerical spectrograms are compared. Even though the log scale colorbar are different in magnitude, it is possible to see the similar decay in energy in high frequency band in the system as injection continues and the system gets more permeable due to the fracturation process. Furthermore, using curve fitting and comparing the fitting parameters we can identify characteristic domain marked by power laws -18-in the spectral domain, and perform a quantitative comparison between experiments and simulations.

Looking at the power spectra in different snapshots of time belonging to different contributions in figures 10,11,12,13, and 14 we can say that initially (t < 0.05 s) air pressure is dominating the signal in total stress due to the resistance to flow in the porous medium. Carving empty areas will give space to the air pressure to penetrate in for air pressure and thus, it will generate vibrations on the plate. These vibrations are particularly in mid-range frequency of the power spectrum from 1-2 kHz up to 20 kHz which was also mentioned in the experimental study [Turkaya et al., 2015]. This indicates that in large times, as the pressure diffuses through the medium, air pressure in the pores become significant in the stress state. The developing channels in this system are fractal in shape, as can be seen in the optical data, and as analyzed in details in Eriksen et al. draft (they have a mass fractal dimension around 1.5 to 1.7, see [Eriksen et al., 2016]). These structures grow from an initially simple geometry and once they are established, the power spectrum of the acoustic emissions become close to power law. The vibration of the air in such cavities, or the solid stress evolving along such boundaries, or the air pressure diffusing from such boundaries, could be responsible for this power law decay.

Right after the initial opening, these peaks diminish and the system starts to be governed by the solid stress variations due to compaction up to the moment when the solid forces exerting on the plate is stable. On the plots, it is possible to see different power law regimes, fitted by straight curves. An exponent (slope) b1 is obtained after a low pass filter of 23 kHz is applied to the spectrogram, and for b2 a high pass filter of 23 kHz is applied. This filter frequency was identified as corresponding to the most marked change of slope. b3 shows the slope of the fit to the entire spectrum. In figures 10 and 11, it is possible to see a clear increase in corner frequency until t=0.25 s where the slope is similar to a power law decay with an exponent equal to b3 which is close to the low frequency slope b1. The corner frequency in figure 10 increases from 15 kHz up to 49 kHz. Divided by the air velocity in empty space this corner frequency evolution corresponds to an evolution in characteristic length (λ = c air f c where c air = 343 m/s) from 2.5 cm down to 7 mm. This can be related to the development of the channel network towards a fractal shape by having thinner outer fingers with time. The fractal shape is scale free therefore we do not see a characteristic peak in the power spectrum.

-19- In solid and total stress (fig. 12 and13), the initial burst during channeling forces solid pressure to have two different slopes. Again as in the cases with air pressure in the empty areas and pores, (even though it comes later due to difference in the dynamics) a power law decay is developed in the power spectrum. As defined in [Eriksen et al., 2016]this experiment belongs to a category in which there is a main dominating channel with some dendritic branches. In the article, the authors noted that the fully developed channel for this category forms t ≈ 1 s. This corresponds well with the time we obtained as a development of a power law decay in power spectrum (due to the fractal dimension of the channel network) in figures 12 and 13.

It is possible to see these evolutions during fracturing in experimental data as seen in figure 14. Analyzing part by part the different sources characteristics in the simulations, we see that both the air pressure vibrations in the cavity becoming fractal, and the solid stress arising from the surrounding grains, are the main source of noise -usually, the solid stress, and once the medium stabilizes, the air vibration. They both present more or less this power law character of the power spectrum, with an exponent around -1.6.

This value is compatible with the experimentally obtained b value -1.61 given in [Turkaya et al., 2015].

Conclusion

Evolution of different source contributions in the acoustic emissions are investigated using numerical simulations. Vibrations due to solid stress variations are found to be governing the power spectra. In the beginning (t < 0.05s) air pressure bursting into the empty area in the cell and eventually reaches/pushes the granular medium creates the major contribution in the low frequency region ( f < 20kHz). Afterward, as the grains getting compacted the solid stress starts to generate powerful signals. Since the solid stress has a step-like loading due to sudden compaction and sudden stress release (if carving occurs after compaction), this loading is stimulating all types of frequencies on the glass plate of the Hele-Shaw cell. In power spectrum, the corner frequency increases due to the initial burst then vanishes due to the fact that the power spectrum having a power law decay.

This power law decay is caused by the fractal dimension of the channel network, inducing all frequencies homogenously and leading a one main slope on the spectrum.

-25-

A: Response of the Channel

The vibrations recorded by the accelerometers are significantly influenced by the spatial variation of these loadings (due to changing channel size) and the response of the accelerometer used. Among these, accelerometer response is corrected as described in [Turkaya et al., 2015] Vibrations affected by the size of the porous apertures are significantly higher than the observable range, therefore not included in these analyses. For the rest of the mentioned sizes, we need to calculate the natural frequency and its harmonics to find the frequency range where the influence of size effect is significant. The natural frequency f n of a single degree of freedom object is defined in [START_REF] Chopra | Dynamics of structures[END_REF] as:

f n = 1 2π k m (A.1)
where k is the stiffness of the spring (in N/m) and m is the mass of the object (in kg). Adapting this formulation to body of a unit size we can say that, it is equal to the p wave velocity is v p = M ρ = k m and thus, the equation A.1 evolves into [START_REF] Royer | Elastic Waves in Solids I: Free and Guided Propagation[END_REF][START_REF] Aki | Quantitative Seismology[END_REF]:

f n = v p 2π (A.2) -26-
However it can be seen that eq. A.2 still does not take the effect of the harmonics of f n into account. Therefore a geometrical factor needs to be added in all dimensions which modifies eq. A.2 into:

f n (i, j, q) = v p 2π 2i -1 L x 2 + 2 j -1 L y 2 + 2q -1 L z 2 (A.3)
where the modes are defined with i, j, q = 1, 2, ..., and L x ,L y ,L z are the dimensions of the channel. Solutions of this equation for different channel sizes are presented in Fig.

A.1. fracturation. It is possible to say that, the bump in the low frequency range (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) is spreading towards a wider frequency band as the channel sizes are getting smaller and smaller. It can also be said that, the fundamental mode of the glass plate is only effective in the high frequency vibrations. Furthermore, we observed that none of the channel size based frequency ranges are corresponding to the drop in the power spectra (around 35 kHz), however the small peak right after 35 kHz can be explained with the peak in red and green curves (corresponding to the big gap and the main channel) in the bottom subfigure.

3.6. Draft Article: Explanation of Earthquake Types using Lab-scale Experiments 107

Draft Article: Explanation of Earthquake Types using Lab-scale Experiments

This study is mainly conducted by the main author under the supervision of Renaud Toussaint. Then, to develop the work further, we discussed several times with all co-authors. All co-authors have contributed with some suggestions on references, scientific methods on analysis and comparison, and presentation of the obtained results. This draft article is in preparation to be submitted to a journal of general scientific interest (e.g. Science, Nature) in early 2017. Please consult the latest version. In this work, we compared the characteristic properties of the acoustic emissions presented earlier in this chapter with the real scale earthquakes. The earthquakes are classified as Type-A and Type-B (there are several other classes as well but we are focusing on these two in this article) are corresponding to the Type-2 and Type-1 acoustic emissions respectively. In this article, we compared the mechanism of how these earthquakes (and the acoustic emissions) are generated and their characteristic features in the power spectrum. We concluded that the sudden movements of the earth generating high frequency earthquakes (Type-A) are similar to the Type-2 stick-slip motion that we observed in the Hele-Shaw cell.

In the power spectrum of these events, it is possible to see several peaks spread over wide frequency ranges. If we take into account the scale difference in these two phenomena it is possible to conclude that they are similar. Furthermore, we compared the Type-A earthquakes which are dominated by the low frequency energy with Type-1 events that are obtained in the Hele-Shaw cell. As pressurized fluid starts to interact with the porous medium (compaction, fracturing, channeling) these events start to emerge and decay as these interactions continue. Moreover, we observed that Type-1 events take longer time than Type-2 similar to the difference between Type-A and Type-B events.

used up to 1000 fps. Using optical and acoustic datasets and numerical simulations, the mechanics leading Type-A and Type-B earthquakes are explained and the results are shown to be compatible with the power spectral signature of the real earthquakes.

Introduction

The mechanics causing earthquakes are always a challenging issue to completely understand how the planet we live upon is evolving. From small magnitudes to large magnitudes the earthquakes are classified with respect to different characteristic properties (1)(2)(3)(4)(5)(6)(7). Among these In figure 1, it is possible to see the spectrograms belonging to Type A and Type B events obtained from volcanoes as presented in (8). The Type B event (top) has an exponential decay in spectrogram similar to pressure diffusion mechanics explained in (22). On the other hand there is a clear burst in the Type A event (bottom) starting and ending rapidly similar to a pulse shape.

2

In this research, we developed an experimental setup to generate Type A and Type B events in lab scale. Fluid flow into a porous medium in a Hele-Shaw cell is an easy way to reproduce source mechanics of small scale Type B events. This type of experiments are previously researched experimentally and numerically by several scientists (9)(10)(11)(12)(13)(14). To understand the percolation and fracturation due to bubble rise inside immersed porous medium several vertical Hele-Shaw cell studies are made (15,16). Moreover, flow rate and solid fractions causing different type of response in the porous medium are investigated in (17). However, instead of the optical response of the porous medium, in this research we focus on the acoustic emissions which are the direct representation of earthquakes at such small scales. Therefore, we used accelerometers (18) to record these emissions during the experiment in various locations (possibly leading a triangulation based source detection (19,20)). The generated signals are traveling through the glass plates and eventually recorded on the different sensors. Even though an ordinary Hele-Shaw cell as in described in ( 13) is largely enough to represent Type B events, it does not have enough confinement permits observing pressure diffusion that we would like to observe. Thus, a rectangular Hele-Shaw cell, having 3 boundaries are sealed and a semipermeable (fluid can escape the system while the grains are restrained inside) 4 th boundary made of a 50 µm steel mesh is designed. This design allows a relaxation / reorganization phase inside the porous medium due to air pressure diffusion into the system, changing the stability at normal stress and shear stress due to increasing pore pressure forcing the some group of grains into stick slip motion until they find a more stable form (18). These stick-slip motions is a good lab scale replica of high frequency, Type A, earthquakes. In figure 2, it is possible to see this experimental setup modified after Turkaya et al. 2015 (18). The Hele-Shaw cell is placed horizontally under the high speed camera. Below the Hele-Shaw cell, in different positions on the glass plate, accelerometers are stick using silica gel. The whole system is activated in the same instant to have synchronization between optical and acoustic recordings using a TTL sig- et al. 2015 (18). A high speed camera is recording the Hele-Shaw cell from above. The accelerometers are stick on the bottom glass of the Hele-Shaw cell. The recorded signal is amplified and recorded in a computer. A signal generator is used to trigger optical and acoustic recordings at the same time to ensure synchronization.

In figure 3, the evolution of the porous medium is presented. In (A), initial state of the porous medium is shown. In (B), after the injection starts, initially flat solid fluid interface starts moving due to the developing pressure. The porous medium inside the Hele-Shaw cell compacts 5 giving more empty area to air, eventually starts fracturing and channeling if the pressure is high enough (C). Different pressures leading different types of structures in this Hele-Shaw cell are investigated by Eriksen et al. 2016 (21). Furthermore, in a general 2D porous medium stress and stability analysis is made by Eriksen et al. 2015. (17). As injection continues, the channeling structure continues to penetrate inside the medium including thinner and thinner branches to the main channel. Eventually a fully developed fracture network having a fractal dimension (21,22) is carved inside the porous medium. 6 Equation 2, can be modified to show the critical pore pressure required to make grains slide as:

p > σ T n - σ T s µ . (3)
Right after a channel is stopped moving, the diffusion equation governs the pore pressure around the channel the details of this evolution is given in (18,32). The frequency of the aftershock occurrence of the Type-A events are found to be decaying with time. Therefore, an Omori-like decay was expected in the results we obtained as well. It has been observed that the lab-scale Type-A events are following an Modified Omori Law with parameters c ′ = 0.0066 ± 0.55 and p = 0.55 ± 0.38. To check the validity of the Modified Omori Law we can derive the parameters c ′ and p theoretically. The Omori parameters found in this experiment (p ≈ 0.5 and c ≈ 0) are compatible with the real scale studies. Nur and Booker found p = 0.5 in their study of earthquakes occurred due to a step in the pore fluid pressure. In the work of Yamashita about the early stage induced secondary aftershocks a p value equal to 0.48 is found (33). There are also other studies linking p value with the pore fluid induced aftershocks eventually finding similar p values to our experimental results (34)(35)(36)(37).

Experimentally reproduced high frequency event can be seen in figure 4. In top subfigure the signal recorded in the accelerometer 1 during the event is shown. The amplitude is rapidly developed as seen in Type A earthquakes. In the bottom, normalized power spectrum of the whole event is seen. It is possible to see that the energy is spread to the wide frequency range.

Therefore, these events are dominated by high frequencies, having a mean power spectral frequency for this event f mean = 224 kHz. All of these characteristic features are corresponding to the features of real earthquakes defined in Lin et al. 2007 (8).

Low Frequency Acoustic Emissions (Type-B)

As discussed in the introduction part, there are emergent earthquakes dominated by the low frequencies. In the experiments with a Hele-Shaw cell, it is possible to reproduce this kind of signals as well. As discussed in (18), after the channeling starts, it is possible to see a bump in the low frequency region (up to 20 kHz) of the power spectrum of the recorded signal. While the injection continues to carve newer and thinner branches and thus, compacts the material this bump spreads over the higher frequencies, eventually diminishes to the noise level when the fracturation is finished. This spreading is a very characteristic feature of the Type-B earthquakes as well. Another finding from the same experiment is that the frequency of occurrence of low frequency acoustic events (defined as Type 1 in ( 18)) is diminishing with increasing channeled area (permeability). This phenomenon is also noted during a slow-slip event in Guerrero, Mexico by Frank et al. (38). Figure 5 shows an example of a Type-B event reproduced experimentally. In the top subfigure, it is possible to see two consecutive emergent events starting at t 1 ≈ 0.55 and t 1 ≈ 0.8.

In the bottom, the power spectral distribution of this low frequency event is given. Particularly it is possible to see that there is a plateau in the low frequency containing significant energy in this spectrum. Therefore, these events are dominated by low frequency energy, leading a mean power spectral frequency for this event f mean = 14 kHz. These fluid based events are mainly like pressure bursts following one another.

A poro-elasto-plastic pressure model is combined with the solid stress inside the cell to investigate the pressure evolution inside the Hele-Shaw cell. Convolving the pressure evolution inside the Hele-Shaw cell (the variation of the thrust on the glass plates) with the far field Lamb wave approximation (39) it is possible to simulate lab scale Type-B earthquakes. After the analysis of power spectrum obtained experimental and numerical results, we see that the initially 11 there are two slopes during the burst of air inside the porous medium. This phase continues until the channel structure has a fractal dimension with all its branches. When this structure is formed, the corner frequency disappears and power spectrum obtains a power law decay due to the homogeneity of this fractal form (21,22). Finally the power spectrum diminishes back to the initial noise level as seen in the nature. Using numerical simulations, different characteristic features of pressure and porous medium and different source mechanics of the vibrations are investigated. Air pressure in the empty channels, pore pressure, solid stress and total stress in the cell are computed for vibrations. Initial air burst is generates very powerful signals until the compaction in the solid partition. As the solid part compacts, solid stress starts to govern the signature of the power spectrum. As seen in the experimental results, the corner frequency evolution and power law decay development in the power spectrum has observed in all different source contributions as well (22). Furthermore, results of the numerical simulations are shown that the noise in the input pressure does not effect these mechanics significantly.

Conclusion

In this research, we studied how Type-A (high frequency) and Type-B (low frequency) earthquakes are occurring using experimental and numerical models. Different source mechanics are investigated using computer simulations. Stick-slip nature of the Type-A events are investigated using the experimental setup where the slowly developing air pressure may break the stability of Coulomb friction. Therefore, a sudden slip occurs where the grains are moving as a group to a more stable state as happening in nature of faults. This sudden slip stimulates vibrations of a wide frequency band due to having a box-car shaped source function (40). The rise and decay in the spectograms of the Type-B events due to their emergent nature are explained with fluid bursting inside the emptied zone (e.g. due to compaction, channeling, or fracturing). As the channels are opening, fluid is flushing in vibrating with the characteristics of the size of the 12 emptied zone giving more energetic pulses in low frequencies. Therefore these events are not impulsive but slowly growing. Initially, there is a corner frequency in the power spectrum of lab scale Type-B events. Afterward, with developing channel network having a fractal nature, this corner frequency disappears and a power law decay is observed in the power spectrum. The shape of this channel network is self-similar meaning that small branches resemble the large main structure. The power spectrum of these events does not have a dominant peak but have a power law decay due to the size effect of this network. This evolution is seen both in experimental and numerical studies. Finally, in the steady state, the power spectrum of these events decay back to the initial noise level.

Conclusion and Future Work

The articles provided in this section discuss the acoustic emissions generated during pneumatic fracturing, due to solid-fluid interactions, inside a Hele-Shaw cell. The first objective of this chapter is to describe the experimental methods, i.e. to precise how the equipment is working, how the measuring devices are calibrated. In the first section, we developed a stepwise experimental procedure and sample preparation method. The way to assure the reproducibility of the experimental work are given early in the chapter. Following this section, the analysis of the acoustic emissions is detailed. The objective is to understand how the acoustic emissions are evolving related to the evolution in the porous medium with the effect of air flow inside. Using the optical data acquired inside the Hele-Shaw cell the evolution of the acoustic emissions (particularly in Fourier domain) are investigated. From the images we can observe the evolution of the emptied area, and the development of the channel network. The power spectrum and mean frequency of the acoustic emissions changes with the evolution in the channeled area. As the fracture network grows, the low frequency of the power spectrum diminishes. Therefore, it evolves towards certain characteristics (e.g. permeability), on a time scale similar to the development of the channel network. Furthermore, due to the diffusion based pressure migration into the pre-compacted areas it is possible to see some stress rearrangements inside the porous medium (due to the increased pore pressure stress state of the medium changes) eventually generating stick-slip motions. These acoustic emissions are different than the low frequency dominated motion. Due to the rapid movement of the grains with shearing motion it is possible to say that this particular motion causes pulse-like stress variations leading high frequency earthquake like acoustic emissions. Due to the similarity in the causing source, instead of a main shock these events more similar to aftershocks and thus they follow the statistical rules developed a-priori for the earthquake aftershocks (i.e. Modified Omori Law). Power spectra evolution and diminishing frequency of occurrence are present in real scale microseismic data as well. Appearance of aftershock-like events, mean frequency, and starting point of an Omori Law, can be used straightforwardly in large scale microseismic monitoring of fluid injection and well stimulation. The permeability of the medium can also be directly estimated from the prefactor of the Omori Law.

Using the numerical models the different contributions of the stress variations of the earthquake is explored. Initially, air fills the empty box like area before hitting to the solid interface. In this period (t < 0.05 s), vibrations due to the air pressure governs the power spectrum. As the solid starts deforming, this deformation creates a fully compacted zone around 3-4 cm (compaction zone) starting from the interface. In the numerical studies, we concluded that the main contribution of the vibrations are coming from the solid part compacting. In power spectrum, we saw that the corner frequency increases initially, but as the channel network develops the fractal structure of the channel has a direct impact on the power spectrum, creating a power law decay in the power spectrum of the emissions.

Introduction

Having two glass plates covering the porous medium which permits to see what is inside of the Hele-Shaw cell is one of the advantages we used in this project. A high speed camera (up to 1000 frames per second (fps)), seeing and recording what is going on inside during the fast fluid flow acquires very rich information about solidfluid interactions. The images are recorded from a top-down view over a horizontally placed cell. The camera may see the whole cell during the entire experiment enabling recording information with a high resolution (15 pixels ≈ 1 cm). Eventually, very high resolution data in time and in space are obtained. Interpreting the acquired experimental data first requires an analysis protocol having several steps. After a pretreatment of images (noise reduction, light correction, cropping etc.) the images are used for several purposes. Having a good contrast between the grains and air inside the cell, it is possible to convert images into binary images showing the structure inside. Furthermore, some of the grains were dyed a priori with Indian Ink to have contrast between grains. This step is essential to use image correlation tools such as particle image velocimetry to obtain displacement fields and their evolution over time. Also, for the simulations mentioned in chapter 2 the binary images are used to input the channel network to compute the stress map inside the Hele-Shaw cell. In this present chapter, the image processing methods are used to characterize and locate flow events both in space and in time. The evolution of these characteristic properties, development of fractal shapes inside the medium during fingering and pattern formation are investigated. Furthermore, it is possible to extract the information about how the porous medium is compacting, carving, eroding, and jamming in some of the locations. The different techniques we can use and the quantities possible to extract will be presented in this section. The two main topics in this chapter are:

• The flow regimes, pattern formation during air injection inside a deformable porous medium. The variations inside the pattern formation at different imposed entrance pressure levels are investigated. This is a draft article going to be submitted soon in Physical Review Letters E.

• The deformation and the compaction inside the porous medium. The characteristics of the fractures, (i.e. width, size, time of development) will be explained using the analysis of the experimental data and numerical simulations.

Draft

Article: Pneumatic fractures in confined granular media 127

Draft Article: Pneumatic fractures in confined granular media

The optical datasets acquired during the aerofracturing experiment which is described in section 3.2 are analyzed further using image processing methods. This study was mainly conducted by Fredrik K. Eriksen under the supervision of Renaud Toussaint. The experiments are co-operated with Semih Turkaya. Then, scientific discussions to develop the scientific work, before redaction and publication, were done with all co-authors. All co-authors contributed with some suggestions on references, scientific methods on analysis and comparison, and presentation of the obtained results. This draft article is submitted to Physical Review Letters E in November 2016. Please consult the latest version.

As a result of these analyses we observed that the channel formations seen in our experiments lead to patterns very similar to the patterns seen in nature or in similar experiments, e.g. viscous fingers in empty Hele-Shaw cells, manganese dendrites, or lightning bolts. In addition, we observed a feature of the channel growth, specifically the presence of an active growth zone outside a structure which is "frozen", due to screening of the pressure gradient by the most advanced parts of the structure. This phenomenon is typical for Laplacian growth systems.

The resulting channels in our experiments are different than channels formed in similar systems having open outer boundaries. In the open systems, the channels are generally much smoother and do not have more than 1-2 branches, while in the closed system we see ramified structures.

We noticed that these patterns have a well defined fractal dimension, and that the fractal dimension of the channels is found to be fairly stable along the main part of the structures. The local box-counting method indicates that when there is little or no erosion, the final structures end up with typical fractal dimensions in the range D = 1.53-1.60. This range is similar to the result of the observations of viscous fingers in saturated porous media. Moreover, the analyses show that even if the structures are growing differently with different pressure ranges, the roughness of the air-solid interface are very similar.

PACS numbers: 83.60.Wc, 81.05.Rm, 47.20.Ma

I. INTRODUCTION

Several processes in engineering, industry and earth sciences involve pneumatic (gas) or hydraulic (liquid) fracturing of the soil, which occurs when fluids in the ground are driven to high enough pressures to deform, fracture and generate porosity in the surrounding soil or rock. For example in environmental engineering, pneumatic or hydraulic fracturing is done to enhance the removal of hazardous contaminants in the vadose zone (soil remediation) [1,2], for soil stabilization injection to ensure a solid foundation for structures [3], or in packer tests for project planning, risk assessment and safe construction of dams and tunnels [4]. In industry, hydraulic fracturing is done to enhance oil and gas recovery [5][6][7], CO 2 sequestration [8], water well-and geothermal energy production [9][10][11]. Related natural processes such * Also at Department of Physics, University of Oslo, Oslo, Norway;

Correspondence: eriksen@unistra.fr as subsurface sediment mobilization are studied in earth sciences, where sand injectites, mud diapirs and mud volcanoes are formed due to pore-fluid overpressure [12][13][14][15][16][17].

For example, the Lusi mud volcano in Indonesia is the biggest and most damaging mud volcano in the world [18], having displaced 40 000 people from their homes, and has been active since May 2006. There is an ongoing debate about how it was triggered, i.e. whether it formed naturally by an earthquake or geothermal process [19][20][21][22], or that it is a man-made consequence of a nearby drilling operation by a company probing for natural gas [23].

Fluid injections into granular media has been extensively studied in laboratory experiments and simulations, where a common method to simplify the problem is to confine the experiment within a quasi-2-dimensional geometry, i.e. a Hele-Shaw cell. In [24,25], the decompaction, fluidization regimes, and coupling between air and granular flow was studied in dry granular media in open circular and rectangular cells during air injection at different overpressures. Similar behavior was seen for liquid saturated granular material injected with the same 2 liquid [26], so it is reasonable to assume that studies of pneumatic fracturing also has applications in hydraulic fracturing. The patterns formed during fluid injection into a granular medium, and evolution of the fluid-solid interface, have been found to resemble Diffusion Limited Aggregation patterns (DLA) and viscous fingering [27], a fingering instability that occurs when a less viscous liquid is injected into a porous medium containing a more viscous liquid with which it cannot mix [28]. As mentioned in [25], the main difference between the viscousand granular fingering instabilities is the absence of interfacial tension in the granular case. For example, the stabilizing forces in viscous fingering are surface forces and viscous drag, while in granular fingering it is the build-up of friction between particles and against the confinement. However, both instabilities are driven by the pressure gradient across the defending medium, which is largest on the longest finger tips, making more advanced fingers grow on expense of the less advanced ones. A notable difference between air injection into a dry granular medium and a saturated one is that the overpressure initially diffuses into the packing in the dry case, while it is already a steady-state Laplace field over the defending liquid in the saturated case.

Further, during air injection into liquid saturated granular media and suspensions, the characteristics of emerging patterns and behavior of the media depend on injection rate, and the competition between mobilized friction and surface forces [29][30][31][32][33][34][35][36][37][38][39][40][41][42]. For example, one observes flow regimes such as two phase flow in rigid porous media [38][39][40][41][42], capillary fracturing, stick-slip bubbles and labyrinth patterns [29][30][31][32][33][34][35][36][37]. In the opposite case, during liquid injection into dry granular media [43], the flow behavior goes from stable invasion towards viscous fingering for increasing flow rate and viscosity of the invading fluid. At intermediate conditions, fractures open up inside the invaded region. The same trend is shown in numerical studies for gas injection into granular media containing the same gas [44].

Granular fingering instabilities have also been studied in closed vertical cells, where gravity drives the flow as heavier beads fall down from a granular layer at the top of a lighter fluid layer [45][46][47][48][49][50]. When the beads detach at the front, they form fingers of falling granular material surrounding finger-like bubbles of rising fluid. These fingers are found to coarsen over time until they reach a typical wavelength depending on the interstitial fluid and bead size.

Typically, in all processes involving fluid injection into granular media, there are flow regimes where the medium has either solid-like behavior or fluid-like behavior. This is one of the special properties of granular materials, which also show gas-like behavior in some cases [51].

In this paper, we present an experimental study on flow regimes and pattern formation during air injection into confined granular media. More specifically, we inject air at constant overpressure into a dry granular medium inside a Hele-Shaw cell, where air escapes at the out-let while beads cannot. The motivation of this setup is to study the granular Raleigh-Taylor instability in compacting granular media, and the coupling between compaction and flow. A similar, but smaller system has been studied in numerical simulations by Niebling et al. [44,52]. During air injection at different overpressures they found that fractures grow faster, longer, as well as coarsen with increasing injection pressure. Here, the fracture propagation velocity is roughly constant initially and scales with the square root of the injection pressure. Further, by varying the interstitial fluid viscosity, two flow regimes were identified; one with finely dispersing bubbles and large scale collective motion of particles, the other one with build-up of a compaction front and fracturing. These flow regimes depend, respectively, on whether the particles are primarily accelerated by the imposed pressure gradient in the fluid, or interactions through particle contacts. This in turn depends on the diffusivity of the interstitial fluid pressure in the granular medium. We analyze our findings in light of these results, to investigate what is similar and what is different in our roughly 10 times larger cell with the same cell gap. As opposed to similar experiments, with open outer boundary conditions, after the flow compacts the medium there is no decompaction. We thus expect the material behavior (at high enough overpressure to displace beads) to have a transition from fluid-like to solid-like during experiments, and that eventual invasion patterns will initially resemble viscous fingering in the fluid-like regime, crossing over to stick-slip fracture propagation as the medium becomes more solid-like, until it reaches a final structure as the compacted medium has reached a completely solid-like behavior. What is less obvious, is how the flow patterns in this system change with the injection pressure. By varying the imposed overpressure, we identify and describe the different flow regimes. Due to the confined nature of the experiment, it is thought to be a laboratory analog to pneumatic and hydraulic fracturing of tight rock reservoirs where the free boundary at the surface is very distant from the injection zone. Therefore, new insight into this problem may have industrial applications in addition to increase the understanding of flow and transformations in porous media.

It is also worth to mention a closely related project [53], where acoustic emissions recorded during the experiments are analyzed. Here, it is shown that different stages of the invasion process can be identified acoustically in terms of characteristic frequencies and distinct microseismic events.

II. METHODS

A. Experimental setup

The experimental setup is a linear Hele-Shaw cell, partially filled with Ugelstad spheres, i.e. dry, non-expanded polystyrene beads with a diameter of 80 ➭m ± 1 %. The 3 cell is made out of two rectangular glass plates (80×40×1 cm in length, width and thickness respectively) clamped together on top of each other with an aluminum spacer controlled separation of 1 mm. A cell volume (76×32×0.1 cm) is formed between the plates by an impermeable sealing tape as shown in figure 1, with one of the short sides left open (outlet). Next, beads are filled into the cell by pouring them through the open side until they occupy about 90 % of the cell volume, followed by closing the open side with a semi-permeable filter (a 50 ➭m steel mesh) to keep beads inside the cell while allowing air to escape. The cell is then flipped vertically to place the granular medium against the semi-permeable outlet by using gravity, resulting in a volume packing fraction of approximately ρ s = 0.44 ± 0.04, assumed to be more or less uniform across the medium. This leaves a volume empty of beads on the sealed side of the cell, opposite to the semi-permeable outlet, with a linear air-solid interface. An inlet hole on the sealed side of the cell is connected to a pressurized air tank which lets us inject air at a constant and maintained overpressure, P in = P abs,in -P 0 (absolute pressure -atmospheric pressure), ranging from 5 to 250 kPa. This will force air to move through the granular medium, towards the semi-permeable outlet, where P abs,out = P 0 = 100 kPa, or in terms of overpressure above the atmospheric one, P out = 0.

During experiments, the prepared cell is positioned horizontally. A selected overpressure is set at the pressure tank outlet and verified by a Honeywell pressure sensor with an accuracy of ±4 kPa. The tubing between the pressure tank and the cell inlet is equipped with an electronic valve such that the air injection is started with a digital trigger signal. Positioned above, with a topdown view of the cell, a Photron SA5 high speed camera is started with the same trigger signal, recording the air invasion at a framerate of 1000 images/s and a resolution of 1024×1024 pixels (1 pixel ≈ 0.7 mm in the cell). Light from a 400 W Dedolight studio lamp provides uniform and flicker-free illumination onto the white beads of the medium. A small fraction of the beads (< 10 %) are dyed black with ink to create tracer particles that are used for tracking frame-to-frame deformations in the granular medium. The experiments are run for 10 s, but typically the fracturing and/or compaction of the granular medium takes less than 5 s.

B. Image processing

In analysis of the images from the high speed camera, we investigate flow regimes, characteristics and evolution of the fracturing patterns formed. We perform image processing with Matlab to obtain the information contained in the images. Quantitative analysis of the invasion patterns is done by converting the grayscale raw data into binary images, i.e. images with either black or white pixels, where the white pixels (value=1) represent The granular medium is confined inside the cell by three impermeable boundaries and an airpermeable boundary on the outlet side. The granular medium is placed against the outlet boundary such that it has a linear interface against a region empty of beads on the sealed inlet side, where pressurized air can be injected. B) Side view of the setup showing the high speed camera placed above. The glass plates are clamped together with aluminum framing while the cell gap is controlled with spacers, which are not shown here.

the pattern and the black pixels (value=0) represent the background, or the bead-filled region. We obtain such binary image sequences from the raw data by thresholding each frame with the initial image, such that the pixels having a value less than 30% of the corresponding initial value becomes white and the rest remain black, as shown in figure 2C. In addition, once a pixel is invaded (and made white) it will remain white for the rest of the experiment to correct for mobilized beads flowing inside channels, which occurs due to erosion and fingers merging in some experiments. From the binary images we extract information such as invasion depth over time, average finger thickness as function of depth, and fractal dimensions.

III. RESULTS

Depending on injection pressure (and friction related to packing fraction) we observe the granular medium to exhibit either a solid-like behavior, or a transition from fluid-like to solid-like behavior during air injection. In the solid-like regime, there is no apparent deformation (bead displacements), and the air is reaching the cell outlet by seeping through the network of pores between beads. On the other hand, in the fluid-like regime we observe significant deformation, where beads are displaced by an amount corresponding to several bead sizes, by momentum exchange between air molecules and beads. At sufficient overpressure the granular medium has a behavior much like a viscous liquid being invaded by air, as the flow opens up channels empty of beads, a Saffman-Taylor like instability [28] resulting in dendritic invasion patterns. A stabilizing mechanism of this flow instability is the build-up of friction and stress chains between beads during compaction, so due to the boundary conditions in our experiments we always end up with a solid-like medium towards the end of an experiment.

To explore the various flow regimes we have systematically varied the air injection pressure within the range of 5 kPa to 250 kPa, and we have identified 6 flow regimes occurring in this range. Typical final structures are shown for each category in figure 3, along with a visual representation of the associated range of injection pressures, which is also listed in table I. Videos of channels formed in categories 3 -6 are found in the supplementary material. We have grouped the experimental results into categories from 1 to 6, corresponding to flow regimes defined by final pattern characteristics and how the channels expand. When discussing the growth, new growth dA refers to the area added to the channel between two snapshots, and the active growth zone is the region of the pattern where dA is added, described in terms of extent d behind the most advanced channel tip.

• Category 1 -No deformation:

The porous medium appears rigid (solid-like), with no significant bead displacement or invasion pattern formed.

• Category 2 -Initial compaction without fingers:

The medium is initially compacted during a collective displacement of beads towards the outlet side, without the formation of invasion channels

• Category 3 -Short thin fingers: A short dendritic invasion pattern is formed, where the final length of the channel is less than 50% of L = 70 cm, the initial length of the granular medium, with an average finger thickness less than 2 cm. The channel growth is like fingering with an active growth zone typically within 5 cm from the most advanced tip, i.e. the channel grows primarily on the longest fingers, while it is frozen in the parts more than 5 cm behind. The new growth as function of distance behind the longest finger typically peaks halfway into the growth zone.

• Category 4 -Long thin fingers: A long dendritic invasion pattern is formed, where the final length of the channel is more than 50% of L, with an average finger thickness less than 2 cm. The channel growth is like fingering with an active growth zone typically within 10 cm of the most advanced finger tip. The amount of new growth as function of distance behind the longest finger typically peaks halfway into the growth zone.

• Category 5 -Long fingers with thicker main channel: A long dendritic invasion pattern is formed with a thicker main channel than the branch fingers, i.e. due to erosion inside the channel. The final length is more than 50% of L, with an average finger thickness between 2 cm and 4 cm. The channel growth is like fingering with an active growth zone typically within 15 cm of the most advanced finger tip. The amount of new growth is typically more spread out as function of distance behind the longest finger, where the growth between 10-15 cm behind the most advanced tip is mostly due to erosion. The erosion is not persistent, and rarely occurs more than 15 cm behind the most advanced tip.

• Category 6 -Long dominating channel: A long and thick main channel with some dendritic branches is formed. The width of the main channel increases by erosion (i.e. motion of beads nonperpendicular to the channel interface, but with a major tangential component along it) and with side channels and branches merging into it. The final length is more than 50% of L, with an average finger thickness more than 4 cm. The channel growth is more like fracturing, with an active growth zone spanning the entire length of the channel. Erosion and finger merging is observed to occur even at distances of 30-40 cm behind the most advanced tip. Again, erosion is not persistent and stops around the point when the channel reaches 90 % of its final length. When the erosion stops, we observe that the main channel walls and side fingers are slowly expanding perpendicularly to the average flow direction.

The pressure ranges for the categories overlap, indicating diffuse boundaries between flow regimes, probably in part due to randomness in the initial bead configuration. In addition, table I lists the typical values of P in we used for experiments within a category, and how many experiments were done. Plots of new growth between snapshots, i.e. area added to the channel between two snapshots, are shown for typical experiments in figure 4. In the plots, dA/W is the ratio of new growth to the cell width W plotted as function of distance behind the most advanced channel tip, where dA is found by image subtraction of binary images of the patterns, dA = A(t + ∆t) -A(t). The snapshots are separated by the durations the channels grow in 10 % increments of their final lengths. The categories 4 and 5 are not fundamentally different, but since category 5 has erosion and category 4 does not, we make separate categories for the purpose of averaging results among similar experiments. For categories 3-5 the growth typically extends in a region of 10-15 cm behind the most advanced tip, whereas for category 6 the growth extends along the whole channel length. The category 6 channels have a rapid and large volume change compared to the other ones, with a high air flow velocity resulting in erosion along the entire channel. In addition, the rapid expansion could let the channel pressure remain somewhat lower than P in during growth, such that when it stops the pressure rises to P in . If it is so, pressure gradients could form at the solid-air interface far behind the most advanced finger explaining the small growth there. The displacement of the granular medium outside the channels is similar for all categories;

In front of channels and branches, beads are pushed like a piston and compacted until the system is jammed. Beads are also pushed out from the side of fingers in the active growth zone, opening up for branches.

In the rest of the paper, we focus on the categories 3 -6, which are the flow regimes where invasion fingers are formed. For these flow regimes, we discuss the typical characteristics of the patterns formed and their growth dynamics. In figure 5 we see the average finger thickness w of the final structures as function of depth x into the granular medium, measured from the initial boundary position. The average finger thickness of a pattern at a given depth into the granular medium is found by intersecting the binary image of the pattern with a line perpendicular to the average flow direction at that depth, and count the total number of white pixels intersected. In other words, the total width of the pattern is measured followed by dividing it by the number of connected pixel groups, i.e. the number of fingers. This is done for each experiment and the results are averaged within each category. We see that the patterns in categories 3 and 4 have roughly the same thickness, averaging around 1.4 cm and is more or less constant as function of depth. We do not observe erosion inside these channels, so a finger thickness around 1.4 cm seems typical when there is no erosion. The category 5 patterns show initially thicker fingers along the first half of the structures, here with an average thickness of 3.3 cm, which then decreases down towards the typical thickness seen for categories 3 and 4. The larger initial finger thickness is due to erosion inside the main channel, meaning that beads near the channel walls are re-mobilized in the flow direction after the initial growth, making the channel grow perpendicular to it. As average thicknesses here suggest, this effect typically increase the channel width by an amount of around 2 cm. The category 6 patterns typically have thick fingers throughout their entire length up to the tips, with an average of 5.3 cm. This thickness is achieved both by erosion inside the main channel, and with fingers merging together. After fingers merge and trap clusters of beads, these beads are mobilized in the flow direction. In general we see an increase in finger thickness and length for increasing pressures. Some examples of both erosion and finger merging are shown in figure 6 for an experiment with P in = 200 kPa.

The plot in figure 7 shows the average number of fingers N f as function of depth x into the medium (found as the number of connected pixel groups on vertical crosssections, as explained for figure 5) for the final structures in each category. For patterns in the categories 3 and 4, the number of fingers increase at lower depths than for patterns in the categories 5 and 6, where there is a more pronounced initial region with fewer fingers, crossing over to more fingers during the last half of their length. This could indicate that channels with few branches form when the friction is low compared to the driving force ∝ ∇P , with a cross-over to the patterns branching out with more fingers as the friction is becoming comparable to the driving force due to compaction. The sudden drop in N f at the deepest parts of the patterns is due to the finite size of the structures, where only a few of the longest fingers reach.

To get a description of the shape of finger tips, we look at how the finger thickness w grow as function of distance d from the finger tip, which is found to follow a power law within the length of the tip region [27] (where the finger thickness increase with d before it saturates):

w(d) ∝ d β . ( 1
)
When w is plotted as a function of d in a log-log plot, the slope β for w(d) within the tip region reveal information about the shape of the tip; it has a more rounded shape for β < 1, while it has a pointy cusp-shape for β > 1. We have averaged w(d) over time for each experiment, i.e. obtaining the average finger tip shape during flow in each experiment. Figure 8A shows a log-log plot of the finger thickness as function of distance from the finger tip, averaged per category. The collapsed average slopes from all experiments seem to fall along the same line with β = 0.68, indicating more rounded finger tips with the same shape in all experiments. However, a purely round profile would give β = 0.5, so the fact that we measure a bigger value means that there still is a singularity in curvature towards the tip -even if it is not a spike shape. The singularity can be shown as follows, where

w ∼ d β ⇒ d ∼ w 1 β ⇒ dd dw ∼ w 1 β -1 , (2)
gives the slope of the front. From this we find the expression for the curvature as,

κ ∼ d 2 d dw 2 ∼ w 1 β -2 (or d 1-2β ). (3)
If β > 0.5, 1 β -2 < 0 and κ diverges at the tip where w (and d) → 0. At the same time, if β < 1, the slope does not diverge at d = 0, and yields a more rounded profile. The new growth is measured on intervals between successive snapshots where the channel length grows 10 % of the final channel length for each snapshot. The large initial new growth (0 % -10 %) in all plots corresponds to the initial bulge that forms before the channeling instability begins. The plots indicate the depth of the active growth zone, which goes from very narrow for category 3 (the growth is focused on the most advanced tips), to the whole channel length for category 6.

We have not observed any cusps as reported in [27] with β = 1.43±0.2, consistent with theory on fluid fingering in the zero-surface tension limit. The evolution of β over time is plotted in 8B, which shows that for all experiments, β fluctuates around the estimated value β = 0.68 over time, indicating that the tip shapes does not seem to change significantly during the different stages of the experiments.

In figure 9A, the finger tip position averaged per category is plotted as function of time, which is found by recording the maximum depth of the patterns in the binary images at each snapshot. We see that for all experiments, the channel length grows linearly with time initially before it decelerates towards a final length, and that both the growth velocity and final depth of the in-vasion structures increase with increasing injection pressure. The initial constant growth rate seems to scale with the injection pressure as P in to show that they fall along the same line initially. In the simulations done by Niebling et al. [44,52] of a similar, but smaller system they too observe an initial linear growth velocity crossing over to decay over time, however we here observe the initial growth velocity to scale with P in as v ∝ P 3 -6). We see that for the structures in Cat. 5 (erosion) and Cat. 6 (erosion and finger merging) the fingers are thicker -up to the finger tip for Cat. 6.

with time, v(t) ∝ t -α , where α = 2.5 ± 0.2. By looking at the invasion patterns, we notice that they have self-similar features, e.g. a smaller branch resemble the whole larger pattern. To characterize this feature we analyze fractal dimensions of the final patterns in three different ways. The box-counting dimension D B is found for each final structure by covering the binary image with boxes, i.e. dividing the image into equal squares of sides s, and count the number N of squares that contain a white pixel as function of box size s. For a fractal, the relationship between the number of boxes covering the pattern and their size follow a power law

N (s) ∝ s -DB , (4)
such that the box-counting dimension D B is found as the negative value of the slope of N (s) in a log-log plot [29,[54][55][56]. By obtaining box-counting data over a range of sizes s, we estimate D B from the slope of linear fits between an upper cutoff s = 32 cm (cell width) and a lower cutoff s = 1 cm (typical for thinner fingers). In addition, we obtain the box-counting dimension D F of the front (air-solid interface). This is found in the same way as the box dimension, but with binary images where only pixels on the perimeter of the patterns are white. Finally, we estimate local fractal dimensions D L as function of depth along the structures. To get the local fractal dimension at a given depth, we intersect the structure at that depth with a vertical line (perpendicular to the flow direction) and do a 1-dimensional box counting along that line, i.e. divide the line into pieces of equal length l and count the number N of line segments containing white pixels as function of l. Again, for a fractal we have the power law N (l) ∝ l -DL ,

and we find the local fractal dimensions at given depths from slopes of linear fits between l = 16 cm and l = 1 cm. To compare D L with D B we use one of Mandelbrot's rules of thumb [55,57]. It states that the codimension of an intersected set equals the sum of the codimensions of the individual intersecting sets, here given by

E 2 -D L = (E 2 -E 1 ) + (E 2 -D) ⇓ D = D L + 1, ( 6
)
where D is the fractal dimension of the pattern, E 1 = 1 is the dimension of the line intersecting it, and E 2 = 2 is the dimension of the image plane containing the sets. .53 -1.62 [29, 56, 58].

Figure 12 shows the results of box-counting the airsolid interface of the final structures, averaged per category. We see that categories 3, 4 and 5 have interface dimensions D F = 1.40, 1.52, 1.48 respectively which are close to their corresponding box-counting dimensions, while category 6 patterns have D F = 1.41 which is about 0.2 less than the corresponding box dimension. This can be interpreted in terms of finger thickness, i.e. the front becomes more similar to the structure for decreasing finger thickness.

Figure 13 shows the local fractal dimensions D L as function of depth into the granular medium. For all flow categories, the curves show a more or less constant behavior over most of the depths (the initial high dimension is due to a larger initial area empty of beads, and the final lower dimension is due to finite size effects of the patterns), which indicates that the patterns have consistent local fractal dimensions over a range of x. The categories 3, 4 and 5 patterns have typical local dimensions of 1.53, 1.58, 1.60, respectively. This is within the range of dimensions seen for viscous fingers in saturated

IV. SCALING OF GROWTH DYNAMICS

As we have shown in figures 9 and 10, the growth of the channel length is linear with time initially until it crosses over to decay towards a final channel length. More specifically, the growth velocity is constant initially, before it crosses over to a power law decay with time. The constant initial velocity v 0 is found to scale with injection pressure P in as

v 0 = C • P 3 2 in , (7)
where C ≈ 10 -3 2 cm/(s•kPa in , where Pin = 64, 110, 134 and 213 kPa, for category 3, 4, 5 and 6 respectively. The initial constant plateau correspond to the period with linear growth over time, which then crosses-over to a declining slope close to -2.5 (dashed line) for all the categories. The slopes show that the decay in growth velocity follows a power law equation with time. in along the ordinate axis in figure 10. By inserting the value for C into equation (7), we can instead collapse the curves along the ordinate axis with v(t)/v 0 resulting in dimensionless and normalized units for the velocity. In order to collapse the curves along the time axis, we define a critical time t = t c as the time when the growth velocity crosses over from being constant to follow a power law decay with time. At the critical time, the finger tip position can be described as

x c = x(t c ) = v 0 • t c , (8)
which gives

t c = x c v 0 = x c C • P 3 2 in . (9)
Further, we check if x c can be described as a typical fraction R of x f , the final length of the channel, such that

x c = R • x f ∝ x f , which gives t c = R • x f C • P 3 2 in ∝ x f P 3 2 in . ( 10
)
By plotting v(t * )/v 0 as function of

t * = R • (t/t c ) = R • ((t • C • P 3 2
in )/x f ) in a log-log plot, we do indeed find that the curves are collapsed along the time axis with a cross-over point corresponding to R ≈ 0.6 (see figure in supplementary material). Since we now have the constants C and R, we can calculate v 0 and t c for individual experiments by inserting the respective P in and x f into equations ( 7) and (10). Figure 14 shows log-log plots of v ′ = v(t ′ )/v 0 as function of t ′ = t/t c for individual experiments with injection pressures in the range of 50 -250 kPa, which follow the same dimensionless curve given by the function

f v (t ′ ) = 1, if t ′ ≤ 1 t ′-α , if t ′ > 1, where α = 2.5. ( 11
)
Similarly, we collapse the fingertip position x ′ = x(t ′ )/v 0 •t c and plot it as function of t ′ = t/t c for the same individual experiments, also shown in figure 14 (Note that the data is from single experiments, i.e. not the average values as in figure 10, explaining the increased amount of noise). Here, the collapsed data follow the curve given by the function

f x (t ′ ) = t ′ , if t ′ ≤ 1 1 1-α t ′1-α + α α-1 , if t ′ > 1, (12)
which is found by integrating each part of equation ( 11), requiring that f x (t ′ > 1) → 1 when t ′ → 1. Further, f x (∞) = α/(α -1) = x f /x c = 1/R, giving R = (α -1)/α = 0.6, which fits very well with our observations. By substituting α into equation (10) we get,

t c = α -1 α • x f C • P 3 2 in . ( 13
)
With equations (7 -13) we propose a description of the fundamental dynamics of the system. Since we have experimentally estimated α and equation (7) relating the initial growth velocity to injection pressure, we can estimate the evolution of the channel length during the invasion for a given injection pressure and final channel length.

V. DISCUSSION AND CONCLUSION

The channel formations seen in our experiments result in patterns very similar to other natural patterns arising from Laplacian growth, e.g. DLA clusters, viscous fingers in empty Hele-Shaw cells, manganese dendrites, or lightning bolts. In such systems, ramified structures expand at a rate proportional to the gradient of a Laplacian (∇ 2 φ = 0) potential field. This type of patterns is in the DLA universality class, where a fractal dimension of D = 1.71 is expected [40,56]. However, the fractal dimensions found for our patterns mainly take values of D ∈ [1. 53 -1.60], which is more similar to viscous fingers in porous Hele-Shaw cells. in and tc = ((α -1)/α) • (x f /v0(Pin)), where C = 10 -3 2 cm/(s•kPa (DBM), where the interfacial growth rate is proportional to the pressure gradient of a power η higher than 1, i.e. v ∝ (∇P ) η , where η = 2 for viscous fingers in porous media [40,56]. In our system the potential field is the overpressure in the medium, which diffuses into the medium initially (not Laplacian). However, the diffusing pressure field quickly approaches the Laplace solution outside finger tips (after 0.4-0.5 s) and reaches a global steady state on the order of a second. In addition, we observe a feature of the channel growth which is typical for Laplacian growth systems, i.e. there is an active growth zone outside a frozen structure, due to screening of the potential 13 gradient by the most advanced parts of the structure. The resulting channels in our experiments are different than channels formed in similar systems having open outer boundaries, such as in [24][25][26]. In the open systems, the channels are generally much smoother and does not have more than 1-2 branches, while in the closed system we see ramified structures. This is thought to be a consequence of the outer boundary conditions which directly influence how the medium can be deformed. In both systems, the beads are compacted and pushed in front of the growing channel like a piston, but in the open system there is also a decompaction front moving inwards from the open outlet, easing the further displacement of beads. On the other hand, in the closed system, after the compaction front hits the outer boundary (which happens between t = 150 -300 ms) further deformation is more difficult due to compaction. Thus, the local configuration of the pressure field on the tips could play a larger role on where the channels grow in this case. In addition, the overpressures imposed for channeling in our experiments (P in = 40-250 kPa) are one order of magnitude larger than the ones for these open systems (P in = 2-20 kPa). Although the cells in [24,26] are smaller than ours, the one in [25] is of comparable dimensions.

For the growth dynamics, the dependency of v 0 on P in , and the exponent α could be system dependent, i.e. the parameters could depend on e.g. bead size and confinement. Thus, it is of interest to investigate further how the dynamics are changing in various confined granular media, and see how the Eqs. (7 -13) apply in those cases. The dependency of x f on P in is also of interest to investigate further, such that if system dependent parameters are known, the expected result of an air injection could be estimated by the injection pressure alone. A brief discussion of what we could find from our experiments is presented in figure 15. Another interesting feature found to be common for all experiments in the system is the typical profile of the finger tips, corresponding to a growth in finger thickness w as function of distance d behind the tip as w ∝ d β where β = 0.68 on average.

As shown in figure 13, the fractal dimension of the channels is found to be fairly stable along the main part of the structures, even for the thick fingers in category 6. The local box-counting method indicates that when there is little or no erosion, as is the case for categories 3 -5, the final structures end up with typical fractal dimensions between D = 1.53 and 1.60. This range is similar to the one observed for viscous fingers in saturated porous media, i.e. D = 1.53 -1.62, suggesting that the channels formed in our experiments are in the DBM universality class of patterns where the interface expands as (∇P ) 2 . The thicker category 6 structures have a fractal dimension of D = 1.76, which is closer to DLA-patterns and viscous fingering in empty saturated Hele-Shaw cells (D = 1.71). However, the higher fractal dimension is probably also an effect of box-counting inside the main channel which is rather space filling compared to the system size, and this makes the box-counting slopes for category 6 more uncertain than for thinner fingers. Furthermore, the global box-counting dimensions D B are slightly lower than the locally estimated fractal dimensions, which could be a finite size effect since less developed outer parts are also taken into account. However, for the most developed and best preserved channels (Cat. 4 and 5), we find that D B = 1.53 and 1.54, which is still in good agreement with the fractal dimensions for viscous fingers in porous media. Finally, the fractal dimensions of the front is found to lie between D F = 1.41 and 1.52 for all categories. This indicates that even if the structures themselves appear to change with increasing injection pressure, the roughness of the air-solid interface seems to be more or less the same.

Conclusion and Future Work

In this section, the objective was to interpret the optical experimental results obtained via the high speed camera. The results show that there are two regimes in this experimental setup. The response of the porous medium vary depending on the pressure and packing fraction (and/or friction). ( 1) Solid regime where the porous medium is totally stable, where the fluid is seeping through the solid part to the atmosphere. ( 2) Fluid phase where we observe significant deformation, channeling, fingering and compaction. Due to the boundary conditions (3 sealed boundaries and 1 fluid-permeable boundary) as the fracture state of the medium reaches its final state, with the effect of compaction the system becomes more rigid (jamming of grains) eventually the behavior of the porous medium converges to the solid regime. Furthermore, 6 different categories (having ≈ 50 kPa pressure difference in each step) are defined depending on the size and thickness of the channels. The pattern of the channeled system was also investigated. It has been observed that the smaller channels resemble the whole larger pattern. The fracture patterns are found to follow a fractal dimension. This fracture pattern does not exist initially, which gradually sets in and defines the behavior of the system as mentioned in the previous chapter using the power spectral evolution in the recorded acoustic emissions.

There is also an on going work investigating the deformations of the porous medium inside the Hele-Shaw cell during aerofracturing, based on optical data analysis. Stress, strain and deformation maps of the porous medium inside the Hele-Shaw cell are compared with the associated pressure gradients. Optically obtained displacement maps are combined and compared with a numerical model based on fluid diffusion. Mechanics causing small and large scale displacements are investigated.

Chapter 5

Finding the Source of the Acoustic Emissions We measure and compare the accuracy of these techniques as a function of sampling rate, dynamic range, array geometry, signal to noise ratio and computational time. We show that this new and very versatile technique works better than the conventional techniques mentioned over a range of sampling rates from 8 kHz to 1 MHz. It is possible to have a decent resolution (3 cm mean error) using a very cheap equipment set. The numerical simulations allow us to track the contributions of different error sources in different methods. The effect of the reflections is also included in our simulation by using the imaginary sources outside the plate boundaries. Moreover, the possibility of extending this method for applications in three dimensional environments is discussed. The implementation in 3D will allow this particular localization method to be used in many other areas, such as monitoring industrial activities (e.g boreholes drilling/production activities) or natural systems (e.g earthquakes, volcanoes, avalanches).

Following this method development, in the second part, we applied this mentioned energy based wave localization method on the stick-slip motion based acoustic emissions during pneumatic fracturing in a Hele-Shaw cell. The result of this localization is compared with the optical localization using digital image correlation between the images right before the start and right after the end of the acoustic signal. The characteristics of the mechanics leading to these events (i.e. source mechanics) are compared using the displacement recordings and polarization of the signals. Shear and normal stress based events are identified and discussed.

Microseismic event monitoring is widely done in the industry to understand and to control the fracturing applications or in safety of the nuclear power plants. The location where the movement is happening, where the deformation is focused, or from where the energy is released have key effects to assess the risks associated with the seismicity. In this work, using the optical data acquired during the air injection experiments inside the Hele-Shaw cell it is possible to detect a map showing the displacement map of the deforming porous medium.

In figure 5.2, the incremental displacement between the images taken before and after an acoustic emission is shown. The displacements (red zones) are focused on the channel tips. Even though this map shows the area which has displaced, it does not specifically show the point which generated the acoustic emissions. For this, we applied the Arrival Time Localization (ATL), Time Reversal Localization (TRL), and Estimating Source Energy Homogeneity (ESEH) -which is particularly developed for the Lamb waves on plates which corresponds well with this type of signal on the glass plates of the Hele-Shaw cell -methods to detect the point (or points, depending on the estimation parameters) where the acoustic emission source is located.

Eventually, we obtained a very good correlation between the displacement maps obtained using digital image correlation and possible source locations (usually a point cloud due to uncertainities in the event duration and the arrival time of the reflected signals) obtained from acoustic localization methods for the events we recorded during the aerofracturing experiments. x + u 2 y where u x,y shows displacements in x and y axes. See the section 5.3 for more details. Source localization of acoustic signals is fundamental for a broad range of academic and industrial applications, from robotics to medicine and from telecommunications to geoscience. 1-7 Even though it is trivial to locate impulsive, non-dispersive signals without any reflections, 8 as the signal and propagating medium get more complex this task becomes more difficult. In this paper, a new approach of localization based on the signal energy is presented. It is possible to implement this estimated source energy homogeneity (ESEH) method into the authors' previous work 9 to localize microseismic events in aero-fracturing of a porous medium in a Hele-Shaw Cell, 10 or other types of 2D experiments such as fracture propagation experiments in plexiglass 11 or crack propagation in paper. 12 Furthermore, it is possible to implement this method into a 3D medium at a larger scale which then can be used in the industrial applications such as geothermal activities, 13-15 well stimulation, 16 etc.

We here present a set of experimental and numerical results that compare different methods of localization: Arrival Time Delay Localization (ATL), 17-19 Energy Based Localization (EBL), 20,21 Time Reversal Localization (TRL) 7,22 (see S1, the supplementary material), and ESEH.

A glass (or plexiglass) plate having dimensions 80 cm × 40 cm × 1 cm is used for experiments. A 4.5 mm diameter steel ball is dropped from 2-3 cm height above the plate on different locations. The signal generated by the impact is recorded with a set of 4 miniature piezoelectric accelerometers (Brüel and Kjaer -4374) with an optimal sensitivity in the (1 Hz-26 kHz) frequency range. The recorded signal is amplified and conditioned using a Brüel and Kjaer Nexus Charge Amplifier-Type 2692-A. Then, the signals are transmitted to the computer using a Ni-DAQ mx PCI-6133 acquisition card. This new method relies on the principle that the source energy, after correction from travel path-related attenuation, a) Electronic mail: turkaya@unistra.fr should be equal (or almost equal) at every sensor location. The energy spreading due to distance for plates of thickness h can be expressed, from energy conservation and isotropy along the plates, as E m = E s 2π Rh where E s is the total emitted source energy and E m is the energy density integrated over time and plate thickness, on a receiver m at a distant R. An additional attenuation factor due to the viscosity of the plate material can be incorporated as E m = E s 2π Rh e -α R m , where α is a damping factor of dimension of inverse length (when this one is lowly or non-dependent of frequency). The source energy can thus be expressed from the energy received at m as

E s = 2πRhE m e -α R m . (1)
Modifying this Eq. ( 1) for localization we have 23,24 E s (r s , r

n ) = ω N y q 0 2πR(n)ρhc(ω) |a(ω)| 2 ω 2 e α(ω)R(n) dω, (2)
where E s (r s , r n ) is the source energy seen from the sensor number n at position r n , the source is at position r s , ρ is the mass density of the plate, c(ω) is the group velocity over different ω (angular frequency) at which the energy is traveling, 24 a(ω) is the acceleration of the plate after Fourier transform, R(n) = r n -r s is the distance between the source and the receiver n, and ω Nyq is the Nyquist frequency (ω Nyq = π dt where dt is a discrete time step). The part e α(ω)R represents correction due to the material based attenuation which is equal to 1 in the case of glass plate. 24 For plexiglass plate for the frequencies lower than the cutoff frequency ω c = ω (1/h), this attenuation coefficient is equal to α(ω) = 1 m -1 , and for the frequencies higher than cutoff frequency it is α(ω) = 7.64 × 10 -4 ω 2/3 m -1 . 24 Defining an average operator as A = 1 N N n ′ =1 A(n ′ ) and using the normalized source energy Ψ s (r s , r n ) = E s (r s , r n )/max(E s (r s , r n )) where max(E s (r s , r n )) is the maximum over the N sensors, we can formulate the standard deviation of the energy recorded at N different sensors as σ(r indicate the position of the source. This method is based on the direct wave (signal without the arrival of reflected waves from the sides). Therefore, it is necessary to carefully define the time window prior to the arrival of reflection so as to achieve high quality estimates of the source position. We define a direct wave window ranging from the impact duration of the source by using Hertzian mechanics (30 µs and 96 µs for the glass and plexiglass plates, respectively) 25,26 up to the time of arrival of the reflected wave, i.e., τ source < L/V R , where L is the length of the raypath of the closest reflected wave (97 µs and 375 µs for the glass and plexiglass plates, respectively). We used window sizes 50 µs on glass and 200 µs on plexiglass for the experimental cases. The Rayleigh wave velocity is V R glass ≈ 3100 m/s and V R plexi ≈ 800 m/s on glass and plexiglass plates. 24 Tests on synthetic data using hertzian contact load during shocks 25,26 are conducted to check different methods (ATL, TRL, EBL, and ESEH) in different conditions before testing in physical experiments (see S2, the supplementary material). During the localization estimations, we used Rayleigh wave velocity V R as the wave velocity. However, it is possible to implement different group velocities for different frequencies in ESEH which decreases the error around ≈10 mm. In Table I, the average error of the different cases of signal generation is presented in this order: infinite and non-dispersive plate, infinite and dispersive plate by taking into account different velocities over the frequency spectrum, finite and nondispersive plate by taking into account the reflected waves from the boundaries, and finite and dispersive plate by taking into account both. Signals for 18 various source locations are generated and localized using different methods. From Table I, it can be seen that the effect of uncertainty due to the estimation of the source energy in the conventional method leads to some estimation error even in the simplest case. In the infinite cases where the reflections are not included, energy based methods are more precise. A clear definition of the primary signal means that ESEH performs well in all circumstances. ATL is comparable to ESEH when the medium is non-dispersive; however, with increased dispersivity the performance of ATL decreases due to the attendant uncertainty in the wave velocity. In Table II, it is possible to compare computational costs for each method in every signal generation context. This table is obtained by measuring the time spent during the estimation process with different methods. The window size of the signal is optimized for the cases to better include the direct wave. As a result, some of the less complex cases have longer signals and thus take more time to estimate. From Table II, it can be seen that ATL is the cheapest method while EBL is the most expensive.

s ) = 1 N N n=1 (Ψ s (r s , r n ) -Ψ s (r s , r n ) )
We tested experimentally the effect of the ball size (1 mm, 5 mm, and 10 mm) and stiffness of the plate (polyamide, glass, and steel) of the signals. We observed that these properties do not affect localization results significantly. The average estimation error for different sampling rates for a set of experiments using a 5 mm steel ball hitting on a plexiglass plate is presented in Figure 1.

On the other hand, due to the lower Rayleigh wave velocity (≈800 m/s) in plexiglass, direct wave is clearer. This allows energy based methods to work in this case with a resolution ≈3.3 ± 1.0 cm, which is calculated by averaging the difference between the real and estimated source position for all the experiments (including the human error on targeting the predefined source position). Analysis of the experiments with different spatial resolution (by changing the grid from 1 cm to 1 mm) has a minor effect on the resolution of the data, but it increases the computational cost 70 times (see S3, the supplementary material). The conventional energy based localization method can be optimized by having a linearly spaced grid of 10 points between the maximum and the minimum possible source energy (see S4, the supplementary material).

The sampling rate of the signal has a significant effect on the estimation process; the smaller the sampling rate, the higher the error on localization. However, EBL and ESEH are the most stable techniques and provide better results over the range of sampling rates tested.

To see the effect of signal to noise ratio (SNR) in localization, Gaussian white noise with different amplitude is added to the input data. In Figure 2, it can be seen that the average error from ATL is the largest when the signal quality is poor. Here, ESEH and EBL have very similar responses. Moreover, no dependence on the propagation medium (i.e., glass or plexiglass) has been observed for the decay of the localization error with respect to the SNR (see S5, the supplementary material).

To see the effect of array geometry on the resolution, we placed the sensors in different positions over the plate in the simulations. Based on our investigations, we found out that a "T" shaped sensor placement gives the best results (≈10% better estimation) (see S6, the supplementary material).

The theoretical formulation of ESEH can be used in 3D media by changing the cylindrical energy spreading (with a decay of the energy in 1/R) in Eq. ( 1) to spherical spreading and modifying the material based viscosity α for the selected 3D medium. Thus, we have for deep sources:

E m = E 4 3 π R 2 e -α R m and for superficial sources: E m = E 2 3 π R 2 e -α R m .
The two main advantages of employing a 3D medium are that (1) the signals reflected from the sides are generally clearly distinguishable (or not existing at all when boundaries are sufficiently distant) and ( 2) the waves are body waves, meaning that the medium is non-dispersive. We thus, anticipate improved resolution in a 3D medium compared to a system of plates.

The proposed method, ESEH, appears to be the best compromise with respect to location accuracy and computational cost among the set of methods tested in this study. Particularly, if the primary signal can clearly be defined, ESEH is very convenient to use. Despite being slower than ATL and TRL, it is more robust when the signal is noisy, dispersive, or fast propagating. This method can be easily applied in 2D environments such as touchscreens, or membranes, 12,27 or in 3D in the petroleum industry, seismology, or medicine.

See the supplementary material S1 for the explanation of the different localization methods used in this study, S2 for detailed information about the numerical simulation, S3 for the details of spatial resolution of EBL, S4 for the optimization of the source energy grid in EBL, S5 for details of the SNR comparison, and S6 for the details of the different array geometries compared.

the correct source position the signals should match giving the highest amplitude after superposition, while a false position causes signals to mismatch. Furthermore, this process is extremely sensitive to the reflected signals which can cause mismatch. A signal, reversed back in time to a grid point can be expressed as:

s ′ n (t, r g ) = s n (-t + |r n -r g | c ) (2)
where r n is the position of the sensor, r g is the position of the grid point, s n is the signal recorded at the sensor n, s ′ n is the signal reversed in time, and c is the group velocity of the signal. Maximum amplitude s n ′′ (r g ) in a grid point r g is computed by stacking the reversed signals using the expression:

s n ′′ (r g ) = max t ([ N n s ′ n (t, r g )] 2 ), ( 3
)
where N is the number of the sensors. The maximum of the stacked amplitude over the grid point gives the estimated source position r loc .

s n ′′ (r loc ) = max rg [s n ′′ (r g )] (4)
In this method, since the superposition of different recorded signals is essential to find the maximum amplitude, having the right propagation velocity is the key to have a higher quality estimation.

C. Energy Based Localization (EBL)

Energy based localization 10-12 uses the geometrical attenuation of the energy to model the signal energy at each receiver location. The ideal result is obtained, exploring a parameter space consisting of modeled source locations and source energies, to estimate a modeled energy as equal (or almost equal) to the received energy (computed using the recorded signal) at each sensor. First, the maximum and minimum possible source energy is calculated using the recorded signal. This is a very significant drawback of this method since it puts extra uncertainty to the system and thus, makes the application computationally expensive.

Finally, the grid point having the least difference between the modeled energy and the energy recorded in each sensor gives the estimated source location. To explain this method 2 mathematically, first the relation between the source energy E s and the modeled energy E m on the location of the sensor n can be expressed:

E m (n) = E s 2πR(n)h , (5)
where R(n) is the distance between the source and the receiver n, and h is the thickness of the plate. The recorded energy is given by:

E w (n) = 1 Z t 2 t 1 |s n (t) -s n | 2 dt, (6)
where Z is a proportionality constant depending on the sensor type and plate thickness.

t 1 and t 2 are the limits of the time window defining the primary signal, selected considering the impact duration of the source and the first arrival of reflections. Since ideally modeled energy E m should be equal to the recorded energy E w at each sensor, the location minimizing the RMS of the difference between the modeled and received energy will give the estimated source position. Therefore the minimum of the RMS error β(x, y, E s ) in the equation:

β(x, y, E s ) = 1 N N n=1 E w (n) -E m (x, y, E s , n) E m (x, y, E s , n) 2 , (7)
will give the estimated position of the source. However, this is a computationally costly approach which checks three parameters, instead of two which is the case in the two first methods introduced.

SUPPLEMENTARY MATERIAL S4 -SOURCE ENERGY RESOLUTION IN EBL

Separately from spatial grid (x,y), the uncertainty in the third dimension (E s ) needs to be considered for EBL. Using the possible maximum and minimum values of R(n) (considering the size of the plate) the bounds for the source energy E s are calculated using the Eq. 5

given in appendix A. Using linearly spaced grid between E smin and E smax the best fitting result is found. As seen in the figures 1 and 2 we found that the grid for E s does not have a significant effect on the estimation result.

SUPPLEMENTARY MATERIAL S5 -SIGNAL TO NOISE RATIO

There is no significant change, except the material based amplitude difference between signals, between glass and plexiglass plates (Fig. 1 and2). Therefore, the response of the localization methods to the added noise is similar (Fig. 3). Using the "awgn" function on matlab, generated white Gaussian noise based on the selected dB is added to the signal.

Since we are using SNR in dB, added noise depend on the power -and thus the amplitude -of the recorded signal. Moreover, in the analysis variations of SNR between sensors also tested and seemed not to have a significant effect on the results. Having a random SNR ratio between (6 -40 dB) for each sensor, after 20 random runs with 18 experimental dataset on glass plate the average error obtained on the signals ESEH 6.4 cm, ATL 9.1 cm, EBL 5.2 cm, and TRL 3.9 cm. In comparison with the signals having any additional noise -ESEH 5.9 cm, ATL 2.4 cm, EBL 4.9 cm, and TRL 3.4 cm -it is possible to say that TRL is the least affected by the variation of noise between sensor and ATL is the most influenced by this change. As we go down in the noise level to SNR 0-6 dB for 5 random runs the mean error values are ESEH 11.8 cm,ATL 25.4 cm,EBL 10.3 cm,TRL 10.9 cm all of them are higher (or very close to) than the random picking number (11 cm) for a grid system of this size.

SUPPLEMENTARY MATERIAL S6 -ARRAY GEOMETRY

Array geometry is efficient to see the difference in the signal in different positions over the plate. In figure 1, four different geometries are presented. Then, in table I the average error corresponding to that particular setup is shown. "T" shape geometry (Geometry 1) provides slightly better result (≈ 10%) than rest of the array geometries. -4-data obtained in sensor 1. An example of application of the STA/LTA procedure on the acoustic signal recorded during on experiment is visible in Figure 2.

Estimation of Source Energy Homogeneity (ESEH)

Estimation of Source Energy Homogeneity considers that the source energy calculated from different recordings should be the same after the correction of energy loss due to the travel-path attenuation (based on material and distance). The geometrical energy spreading can be given with the equation E m = E s 2π Rh where E s is the source energy, E m is the recorded energy at receiver m at a distant R where h is the plate thickness. E s can be estimated following [START_REF] Hibert | Slope instabilities in Dolomieu crater, Réunion Island: From seismic signals to rockfall characteristics[END_REF]Farin et al., 2016;Turkaya et al., 2016] by:

E s (r s , r n ) = ∫ ω N y q 0 2πR(n)ρhc(ω) |a(ω)| 2 ω 2 dω, (2)
where E s (r s , r n ) is the computed energy coming from a source at r s recorded by a sensor at r n . ρ is the unit density of the plate, c(ω) is the group velocity of the waves over different angular frequencies ω up to Nyquist frequency ω N yq , a(ω) is the Fourier transform of the accelerometric recordings. In the experimental setup, 4 sensors are placed on the glass plate on different locations (see figure 1B andC) to have a good spatial coverage. Standard deviation σ(r s ) for source energy E s on different sensors are calculated.

Minimum of this standard deviation indicates the source location [Turkaya et al., 2016].

We used 5 mm grid spacing for the 45 cm × 30 cm area covered by the porous medium inside the Hele-Shaw cell.

Image Processing for Deformation Localization via Optical Data

The optical images captured via a high speed camera (FastCam SA5 -Photron, recorded with 125 fps with a resolution 1024 × 1024 pixels) during the experiments.

These images are investigated to localize the deformation corresponding to the time interval when the acoustic signal is emitted. The deformation of the medium during the acoustic emission is analyzed based on the image analysis techniques eventually producing frame-to-frame displacement fields of the porous medium.

First, two images corresponding to the start and the end of the acoustic event is taken. Then, a pixel coarsening method is applied to reduce the noise in the final data as explained in [Johnsen et al., 2008]. Following this pre-treatment, a Digital Image Cor- LAB software) is applied to obtain frame-to-frame displacement fields. Ncorr cross-correlates subwindows between two images to find displacement fields in between. A detailed description of Ncorr is given in [START_REF] Blaber | Ncorr: open-source 2d digital image correlation matlab software[END_REF].

In figure 3, an example of the displacement magnitude field obtained is given. This map shows the magnitude of the absolute displacement of the medium (without any defined direction in the format of ||u|| = u 2 x + u 2 y where u x,y shows the displacement in directions x and y) occurred in between two images. The blue to red colormap shows the intensity of the displacement over the Hele-Shaw cell. As it can be seen for this snapshot, the displacements are focused on the fingertips. However, it is still possible to see displacement (of lower amplitudes) farther into the porous medium.

Mechanics of the Source

In this section we compare the different source mechanics of the acoustic events using the polarization of the signal. For example, a force occurred on the glass towards a sensor will create positive displacement on the recordings of that particular sensor. Having we can deduce the event type. Since we do not know exact duration of the primary signal, the time window of estimation is increased stepwise to have multiple estimations to be compared with DIC and polarization results.

Results and Discussion

In figures 5 and 6 the localization for two different events are given. In top subfigure, displacement magnitude maps (in gray) are overlaid by acoustic localization results.

Here, amplitude of this displacement is not essential for comparison since a small displacement may also generate an event strong enough to be detected. Therefore, we used a grayscale displacement field to see the compatibility of displacement between optical and acoustic localization. Compared with the Arrival Time Localization (ATL) and Time

Reversal Localization (TRL) method, the results obtained from ESEH method fit better to the DIC results than the ATL and TRL. Each colored dot is calculated for different time window length to define primary signal (see bottom subfigure for window lengths). The size of the time window is set to extend stepwise to give multiple estimations for a single event for a better comparison with DIC results. However since the size of the time windows does not affect the results of TRL and ATL, it is not possible to have multiple estimation results for each event.

To interpret the obtained results seen in fig. 5, we can say that it is possible to see that all the sensors are indicating similar polarizations. This can happen either due to the normal forces applying inside the area covered by the sensors (like a steel ball hit on a glass plate), or shear forces acting outside the sensor coverage area where the red arrow is placed. Fracture opening generates polarization both in positive (compaction) and negative (extension) due to shear forces. In figure 6, the polarizations are indicating shear movement that is generating opposite polarization in E-W axis. Instead of having similar polarizations, as seen in a compaction type event, we observed opposite sign of polarizations in E-W axes due to the shear forces acting on the plate.

Conclusion

In this article, we show that the digital images and acoustic signals are very good -12-image correlation method are compared with the acoustic signal based source localization results. Size of the time windows set to extend freely so that the method would give multiple results for a single event for a better comparison. Depending on the polarization found on different sensors, it is possible to define the type of the source leading acoustic emissions. Localization results obtained from ESEH seem to fit well with displacement map obtained with optical analysis. Different from the traditional methods, ESEH works on emergent signals. This method is the first application of the recently developed ESEH method in complex signals.

Conclusion and Future Work

The objective of this chapter was to find the source of the acoustic emission signals in the Hele-Shaw cell. Since it is a complex structure made of glass plates and steel clamps on 3 edges, we simplified, for the localization method, the structure to a basic glass plate having 4 free edges. To characterize the properties of the plate, we conducted several experiments based on launching a ball towards the plate from 2-3 cm distance. The generated vibrations are recorded at the accelerometers placed on the different positions on the plate. Following the signal generation and recording, we used several localization methods (including a newly developed energy based localization method, ESEH) on the recorded datasets to estimate the source positions which are known a priori. The results obtained with different methods are compared with the real source location. The distance between the real and the estimated source positions (error) obtained for each method are compared as well. A Lamb wave based numerical model is used to compare different cases (i.e. dispersivity and finite size of the medium) as well. In the experimental and numerical results we observed that ESEH works better than the other methods. This difference is clearer when the primary signal (signal without reflections) can be clearly defined. This method works faster than the conventional energy based localization methods.

In the second part of this chapter, we applied ESEH (and other methods) to the aftershock-like acoustic events obtained from the Hele-Shaw cell experiments defined in chapter 3. We compared the estimation results with the displacement magnitude fields (see chapter 3) to see if the both systems are compatible. We investigated different source mechanics (based on shear and normal forces) leading different polarizations (up-down displacements recorded on the sensors) on different receivers. In the events based on the variations of normal forces, all the receivers show similar displacement (if the source is in the area covered by the sensors) whereas in the events based on the variations of the shear forces, a clear distinction in the sign of the displacement recordings (and thus, the polarizations) are present. For future work, the possibilities to extend the application of ESEH in 3D media are presented. It is possible to extend it via changing the cylindrical energy spreading to spherical spreading (for deep sources) or to hemi-spherical spreading for superficial sources. Since in 3D medium the primary signals are clearly distinguishable (usually) and the medium is not dispersive we expect that ESEH will have even better resolution when applied in 3D media.

Chapter 6

Conclusion and Perspectives

During earthquakes the elastic waves are released from a source inside the earth. Like other waves, they travel through the medium, eventually reaching a receiver which can be a human feeling the shaking or a seismometer recording the motion of the ground. However it is not possible to say that all the earthquakes are identical. They are unpredictable events having several types of source mechanics. However, it is possible to understand their source mechanics in a fully controlled (and monitored) laboratory environment by creating lab-scale earthquakes. Replicating small scale earthquakes enables analyzing the mechanics causing seismicity with a very high resolution in time and in space. As explained in chapter 3, in the experiments we used 1 MHz sampling rate in accelerometers and we took images up to 1000 frames per second using the high speed camera.

The experimental work in the Hele-Shaw cell had two parts. In the part focused on the acoustic emissions, we saw that two types of acoustic events exists depending on the solid-fluid interactions inside the plate. ( 1) During the initial compaction phase, we have low frequency emergent acoustic events similar to low frequency earthquakes. (2) Following the main compaction/channeling we have the granular rearrangements due to slowly advancing air pressure diffusion. As diffusion weakens the stabilizing effect of normal stress, several granular groups show stick-slip motions towards a more stable form. These stick-slip motions are generating wide frequency band acoustic emissions similar to high frequency earthquakes. In the following section, these experiments are explained in more detail.

Combining a poro-elasticity based numerical model computing the fluid pressure map of the Hele-Shaw cell with the granular friction model computing solid stress of the porous medium, it is possible to obtain the evolution of the total stress exerting pressure on the glass plate. Using a Green function to generate Lamb waves on the plate, we can synthesize the acoustic emissions. In these simulations we saw that the initial burst of air is effective in sending acoustic waves until the solid deformation dominates stress changes and eventually governs the vibrations. Particularly, the fact that the channel network progressively acquires a fractal structure is found to be reflected in the fact that the power spectrum of the emissions also adopts a scale free (power law) structure. As this fractal dimension sets in, we start to see a power law decay in the power spectrum.

We compared these acoustic findings with digital image correlation results. The images acquired during the experiments are used to locate the incremental displacement maps, stress states etc.

Furthermore, these acoustic events are investigated to find the sources. First, 184 Chapter 6. Conclusion and Perspectives the characteristic properties (e.g. wave velocity, dispersion curve) of the glass plate of the Hele-Shaw cell are determined. A simple experimental setup consisting in a ball hitting on a plate to localize the source position acoustically is designed. Different source localization methods including a newly developed method (Localization Based on Estimation of Source Energy Homogeneity, ESEH) are compared with respect to the distances between the real position of the source and the estimated position using different methods. Following the simple experiments with a steel ball hitting on the glass plate, acoustic events during aerofracturing are localized. However, since the real source of the acoustic events are not known a priori, calibration and validation of the results require some more information. For this, we used the digital image correlation results. Displacement magnitude maps showing where the porous medium was moving during the acoustic emission is used to compare localization results with the optically obtained displacement maps. Finally, we found a good correlation between the acoustic and optical localization.

In the future, can apply the developed methods to investigate the following questions:

• System with different loading: By placing the Hele-Shaw cell vertical (or inclined) it is possible to apply the effect of gravity into the system. By varying the angle the effect of the gravity can be manipulated and investigated further.

• Mobile Top Plate: In the current setup both plates are clamped and unable to move except some small bending. By having the top plate mobile (see appendix) it is possible to study the effect of fluid pressure on the faults.

• Prestress and Shear Loading: The more complex fault conditions can be investigated by applying a prestress and shear loading into the system.

• Fluid extraction: In the experiment we have explained we only have injection experiments. The response of the porous medium to a negative pressure gradient is close to some industrial applications (e.g. petroleum extraction).

In 3D, use the method developed on large scale experiments or phenomena. Notably, vulcanologist are interested in emissions associated with fluid flow coming out of a crystal/fluid sponge in a reservoir, which is the current representation of the conduit flow dynamics.

Furthermore, the acoustic signal localization method (ESEH) can be implemented to 3D media to localize microseismic events. The energy based nature of the method is a very strong point when the reflections are negligible. start capturing images with 1000 fps). The air pressure (1 -2 bars) will be constant during a very short period (0.01 -0.02 sec) then will go back to zero.

There are several things that can be explored in this experiment. Firstly, we check the stress state inside the media using the Particle Imaging Velocimetry. Then, we compare it with the probable acoustic events and the signature in the power spectra. We are planning to play with the duration of the injection (1st exp: 0.01 seconds) and the pressure (1st exp: 1.5 bars) to see the effect on the porous medium.

In addition, at the initial stages, thickness of the steel spacers (1st exp: 0.625 cm), initial compaction of the media via hammering during filling phase (initial exp: no compaction), and function of injection pressure over time (1st exp: no rise time, 0.005 seconds in following experiments) can be varied to see and understand the response of the porous media. At the following stages, when a control over those parameters are achieved, a tilt on Hele-Shaw cell can be introduced with a systematic increase to understand how gravity affects the response.

A.2 Dispersion Curve for Lamb Wave Simulations

In this section we will give some more details about calibration of the Lamb wave simulation in two parts. First, we will give information about the dispersivity calculations of the medium and we will compare this with theoretical approaches. Second, we will discuss the effect of the shear forces in simulations.

A.2.1 Dispersion Curve

Surface waves are dispersive, meaning that different frequencies travel with different wave velocities. In the simulations that we made, we used a theoretical approach explained in chapter 2. In an experimental setup with a glass plate, using a signal generator we sent a signal covering a wide range of frequencies using a piezoelectrical sensor stick on the plate. To cover the effect of different distances, we placed an accelerometer right next to the source (distance smaller than 1 mm) then, we shifted the receiver farther from the source with 1 mm steps until 7 cm. To correct the shooting time, we also recorded the trigger signal in a different channel. After synchronizing the received signals using the trigger signal the time versus distance map (similar to the principle of a seismogram) is obtained in figure A.3: Following this arrival time map, we convert this to the Fourier domain in time and in space to obtain the wavenumber versus the angular frequency plot. In figure A.4 the two theoretically obtained dispersion curves are compared with experimentally obtained one.

Since the non-dispersive regime is not included in the theoretical study [Goyder 1980] we extended the slope of the last point on the dispersive regime to account for the non-dispersive frequency range. Recently, we saw that the experimental dispersion curve is somewhere in between. For the future work, using the experimental dispersion curve will improve the resolution of numerical simulation of acoustic emissions.

Résumé en anglais

Lamb waves, acoustic emissions, fracturing, porous medium, signal localization, earthquake types Fluid induced brittle deformation of porous medium is a phenomenon commonly present in everyday life. From an espresso machine to volcanoes it is possible to see traces of this phenomenon. In a rectangular Hele-Shaw cell we inject air into a loose porous medium. Then, we monitor this system using optical imaging using a high speed camera (1000 fps) and 4 high frequency resolution accelerometers. Using the numerical and experimental acoustic emissions, different sources of the recorded signal (vibrations due to air, changes in the effective stress due to fluid-solid interactions) are analyzed. We found that, the peaks in the low frequency range (f < 20 kHz) diminishes while the medium fractures. Furthermore, we propose a new signal localization method based on energy amplitude attenuation and inversed source amplitude comparison. Furthermore, using optical and acoustic datasets and numerical simulations, the mechanics leading Type-A and Type-B earthquakes are explained.

Figure 1 . 1 :

 11 Figure 1.1: Représentation de la cellule de Hele-Shaw (horizontale) et de l'état initial du milieu poreux, avec une chambre de compression sur la droite où la pression d'air est rapidement augmentée, puis maintenue.

Figure 1

 1 Figure 1.2: L'évolution des canaux dans la cellule de Hele-Shaw (horizontale). (A) On voit l'état initial du milieu poreux. (B) Début de la compaction, l'interface entre les parties solide et fluide est quasi-linéaire. (C) Début de la chenalisation. En plusieurs points des canaux apparaissent. (D) Les canaux se propagent jusqu'au centre du milieu poreux, ils deviennent plus fins avec l'influence de la compaction.

Figure 1 . 3 :

 13 Figure 1.3: Les émissions acoustiques pendant l'aérofracturation. La première photo montre la cellule de Hele-Shaw et présente l'état de fracturation du milieu poreux. L'écoulement du fluide a engendré une forte chenalisation et les canaux créés présentent de nombreuses ramifications. La deuxième figure montre l'amplitude du signal enregistré pendant l'expérience. Il nous renseigne sur l'allure du signal tout au long de l'expérience (durée totale 4 secondes). La troisième figure montre un zoom sur la fenêtre temporelle choisie et qui correspond à la photo de la cellule. Le signal enregistré correspond à un événement de type « stick-slip ». En utilisant ce signal, on peut calculer la puissance spectrale associée. La dernière figure montre la puissance spectrale du signal dans la fenêtre temporelle.

  Article: Aerofracture through a double looking glass, mixing optics and acoustics Cette étude est menée à parts égales par Semih Turkaya (analyse et interprétation des données acoustiques) et Fredrik K. Eriksen (analyse et interprétation des données optiques) sous la supervision de Renaud Toussaint.

Figure 1 . 4 :Figure 1 . 5 :

 1415 Figure 1.4: Les émissions acoustiques avant l'injection. La photo en haut montre la cellule de Hele Shaw à son état initial avant l'injection de l'air. La seconde figure montre le signal reçu par le accéléromètre lors de l'expérience. La troisième figure correspond au signal correspondant à la fenêtre rouge de la seconde figure. Cet enregistrement montre le bruit ambiant présent au sein du laboratoire. La dernière figure présente ce même signal. Il est représenté par son spectre de puissance.

Figure 1 . 6 :

 16 Figure 1.6: Les puissances spectrale des différents snapshots. Les courbes présentées de différentes couleurs correspondent aux puissances spectrales à différents moments de l'expérience. Elles indiquent que la puissance dans les basses fréquences diminuent avec la fracturation. La courbe observée et labélisée « type 2 event » correspond à la courbe de la puissance spectrale caractéristique d'un événement « stick-slip ».

Figure 1 . 8 :

 18 Figure 1.8: Les champs de pression de fluide et pression totale à l'instant t=0.536 s. La photo montre l'état de l'expérience à cet instant et indique les conditions aux limites associées. La seconde image au même instant représente le champ de pression de fluide calculé en utilisant le modèle numérique. La dernière figure montre le champ de pression totale (pression de fluide et contrainte solide) calculé grâce au modèle numérique toujours à cet instant donné.

Figure 1 . 10 :

 110 Figure 1.10: Les images de la cellule de Hele-Shaw après la binarisation des photos. Cela nous permet d'étudier uniquement les structures issues de la chenalisation et ainsi obtenir leurs dimensions fractales.

Figure 2 . 1 :

 21 Figure 2.1: "La Trahison des Images" ("The Treachery of Images") by René Magritte, 1928

Figure 2 . 2 :

 22 Figure 2.2: An image of the experimental setup for signal localization. 4 accelerometers are placed on the glass plate having 80 cm × 40 cm dimensions. A steel ball is fired at the red points to be estimated afterwards using different localization techniques.

Figure 2 . 3 :

 23 Figure 2.3: Example of aerofracturing during the injection experiment. Different snapshots on different times show how the medium is evolving as injection continues.

Figure 2 . 4 :•

 24 Figure 2.4: A drawing of the experimental setup used to investigate acoustic emissions during aerofracturing. Hele-Shaw cell is monitored using accelerometers (under the cell) and the high speed camera (looking from above) during continuous injection of constant air pressure.

Figure 1 :

 1 Figure 1: a) Hele-Shaw Cell used in the experiments. b) Material used in the experiment c) Channeling in the Hele-Shaw Cell

Figure 2 :

 2 Figure 2: Flowchart explaining two different localization processes.

Figure 3 :

 3 Figure 3: Validation of the Localization Methods using a Steel Ball. Figures on the left column shows localization result using the arrival time delay. Figures on the right column shows localization result using the signal energy.

Figure 4 :

 4 Figure 4: Localization of the signal received during an experiment. The figure of the left side shows localization result using the arrival time delay. The figure on the right shows localization result using the signal energy.

Figure 5 :

 5 Figure 5: Image b) shows the intensity difference between the granular packing (grey) and the airfilled region of the cell (black). The top image a) shows the image of the initial packing, and c) shows the emptied pattern (grey) as the difference a) minus b). The bottom image d) shows the segmented binary image from c)

Figure 7 :

 7 Figure 7: Example of a PIV field between the initial packing and the pattern image in figure5. Blue represents unmoved particles and red is the highest velocity magnitude in the frame. Because of the relatively large time separation of the selected frames, the velocity field close to the emptied region (white) is poorly defined. However, the banana-shaped part to the right of the velocity field gives a good visualization of the extent of the compaction front.

Figure 8 :

 8 Figure 8: Signal Recordings and Time Windows. The uppermost image shows the Hele-Shaw cell on the instant of the signal received inside the red rectangle. The signal received during the whole experiment is given in the first plot. The instantaneous signal is transferred to the power spectrum which is shown in the lowermost plot.

Figure 9 :

 9 Figure 9: Signal Received during Air Injection

Figure 10 :

 10 Figure 10: Signal Received after Channel Formation

FIGURE 1 |

 1 FIGURE 1 | Schematic diagram of the Hele-Shaw cell including the dimensions. (A) Top view of the Hele-Shaw cell including its dimensions. The cell is made of two glass plates which are placed on top of the other with 1.5 mm spacing between them. A porous medium is placed inside the cell. Three sides of the cell are sealed and the 4th boundary is sealed with a semi-permeable filter. (B) Side view of the Hele-Shaw cell showing clamp details. To protect the glass plates from point loading, a rubber sheet is placed in between the screw and the plate. Aluminum spacers are placed in between the plates to provide equivalent spacing everywhere in the cell.

FIGURE 2 |

 2 FIGURE 2 | The acquisition chain of the aero-fracturing experiments with a Hele-Shaw cell. The signal acquisition card, camera and the electrovalve connected to the air pump are triggered at the same time via a TTL sent from the signal generator to have synchronized optical and acoustic data. The sensors are placed on the bottom glass plate of the Hele-Shaw cell.

FIGURE 3 |

 3 FIGURE 3 | Image of the Hele-Shaw Cell prior to injection (A), during quasi-homogenous compaction (B) and during channeling (C-D). Red dot shows the position of the acoustic sensor (accelerometer). Yellow circles in (C) represent the locations of the first finger-like carved formations.

FIGURE 4 |

 4 FIGURE 4 | Flowchart showing the acoustic-emission analysis procedure. Snapshots of the experimental signal is taken. Then, they are first converted to the fourier domain to obtain the power spectra. Afterwards, power spectra are compared with each other to understand differences.

FIGURE 5 |

 5 FIGURE 5 | Ambient Lab Noise recorded prior to experiment. The red strip in the top left panel (time window of 0.005 s) and its zoomed version (Bottom left panel) corresponds to the snapshot that the power spectra (Right) occurred at. Right panel: The uncorrected (blue) and corrected with sensor response (red) power spectra of the ambient lab noise.

FIGURE 6 |

 6 FIGURE 6 | Optical (A) and acoustic (B) experimental recordings just before channeling. Yellow circles on the left image shows the initiation points of the channeling. The red dot on the left image illustrates the sensor location where the signal on the right is recorded. Top Right panel: Acoustic signal recorded throughout the entire experiment. The red strip (time window of 0.005 s) corresponds to the snapshot that the image (A) and power spectrum (Bottom Right) occurred at.

FIGURE 7 |

 7 FIGURE 7 | Mean frequency evolution (normalized with the maximum value 200 kHz) during injection. As the channel network develops progressively the energy percentage in the low frequency range f <10 kHz diminishes, and the mean frequency gets higher. A cubic polynomial curve (red) is fitted to the mean frequency data to simplify the comparison and discussions. The location of the inflection point of this curve is indicated by the (red) dash dot strip (t 0 = 1.49 s). In this curve the normalized maximum finger length (green, maximum value is approximately 35 cm) and the carved area (black, maximum value is approximately 270 cm 2 ) is also presented.The area saturates to a maximum at a time close to t 0 . The length still increases slightly after t 0 , in particular, a small step in finger length (high slope of the length curve, corresponding to a large stick slip event connected to the finger) can be noticed around t 0 .

FIGURE 8 |

 8 FIGURE 8 | STA/LTA detection plots for Type 1 (upper figure) and Type 2 (lower figure) events. Top:The raw signal is filtered by using a butterworth bandpass filter with corner frequencies 100 Hz and 10 kHz before computation. Bottom: The raw signal is filtered by using a highpass filter for frequencies higher than 100 kHz before computation. We found that some of the Type 2 events can be classified as Type 1 events since they have energy in the low frequency higher than the threshold.

FIGURE 10 |FIGURE 11 |

 1011 FIGURE 10 | Optical (Left) and acoustic (Right) experimental recordings just before (Top) and during (Bottom) an aftershock like event. The red dot on the left image illustrates the sensor location where the signal on the right is recorded. Top Right: Acoustic signal recorded throughout the entire experiment. The red strip (time window of 0.005 s) corresponds to the snapshot that the image (Left) and power spectrum (Right) occurred at.

FIGURE 12 |

 12 FIGURE 12 | Conventional and Modified Omori Law Fits of the cumulative number of events. The circles represent the cumulative number of events and the solid lines represent the fitting curves given by the Equations (5) and(7). t ′ = 0 in this representation corresponds to time t 0 = 1.49 s after the start of the air injection.

Frontiers

  in Physics | www.frontiersin.org 13 September 2015 | Volume 3 | Article 70

3. 5 .

 5 Draft Article: Numerical Studies of the Acoustic Emissions during Pneumatic Fracturing 75

Figure 1 .

 1 Figure 1. Compaction zone in the Hele-Shaw cell. As the air pushes grains, compacted zone with a length L occurs. Following this zone exponential decrease along the length x on the density of the grains (down to the undeformed phase) is presented.

  considering that the local variations in density of the fluid and porosity are negligible compared to the fluid velocity variations. (ii) -∂ ∂t (ρ f φ) ∼ -φ ∂ ∂t (ρ f )φcρ 0 ∂ ∂t (P), considering porosity is constant over time, variations over time in the density of the fluid is expressed using the compressibility equation -7-(eq. 8). (iii

Figure 2 .

 2 Figure 2. The sequence of the simulation analysis. Signals are generated using the forces exerted on the plate. Then, the evolution of these signals are compared with the evolving state of solid-fluid interactions inside the Hele-Shaw cell.

Figure 3 .

 3 Figure 3. Hele-Shaw Cell during injection t = 0.16 s. Colorbar shows the pressure applied on the glass plate in Pa.

  the dispersion relation (which satisfies k h ≃ 1) then the dispersion relation will change to (i.e. ω k = dω dk = c) where c is the wave velocity. To keep the vibra-

Figure 4 .

 4 Figure 4. Fracture state of the porous medium is compared with the total stress at a given time. Top: White color shows the porous medium, black color shows the carved area, and blue dot on the image shows the point where the loading curve is taken from. Bottom: Bottom: Red star shows the instantaneous pressure. Blue and red curves show the air and total pressure over time respectively.

Figure 5 .

 5 Figure 5. Pressures inside the Hele-Shaw cell, pressure evolution on different points marked on the cell are presented

Figure 6 .

 6 Figure 6. The comparison of the mean frequency evolution inside the cell. Top: Mean frequency vs time for different stress contributions (simulations) in comparison with the mean frequency vs time obtained from experimental data. Bottom: The signal experimentally recorded on the plate during the air injection experiment.

Figure 7 .

 7 Figure 7. Power spectral evolution of the Lamb waves on the plate compared with the different types of loading.

Figure 8 .

 8 Figure 8. Comparison of two different spectrogram. Top: Spectrogram when the input pressure is flat.Bottom: Spectrogram when Gaussian noise is applied to the air input pressure.

Figure 9 .

 9 Figure 9. Power spectral distribution of the signals in different time steps. The contributions of different types of sources inside the Hele-Shaw Cell are indicated on the plot.

Figure 10 .

 10 Figure 10. Power law evolution and comparison of power spectra due to the air pressure in empty channels (ax b = y).

Figure 11 .

 11 Figure 11. Power law evolution and comparison of power spectra due to the air pressure in pores (ax b = y).

Figure 12 .Figure 13 .Figure 14 .

 121314 Figure 12. Power law evolution and comparison of power spectra due to solid stress (ax b = y).

  before the analyses. Influence of the channel size can be simulated by checking the harmonics of the natural frequency of an empty room of a size that particular channel. Here we take into account five different channel sizes:1. Main rectangular empty zone between the inlet and the porous medium ( 30 cm × 10 cm × 1 mm) 2. Main channel linking the rectangular empty zone and the small channels ( 5 cm × 25 cm × 1 mm) 3. Centimetric channels following the main channel and penetrating inside the porous medium ( 1 cm × 1.5 cm × 1 mm) 4. Milimetric aperture following small channels and diffusing the pressure inside the medium ( 1 mm × 1 mm × 1 mm) 5. Porous apertures having roughly size around a grain diameter(80 µm)

Figure A. 1 .

 1 Figure A.1. Figure showing the frequency ranges where the different size effects are present. Top: Several power spectral distributions recorded during the experiment. Bottom: Curves represent the different size effects and the response of the plate modes.

  types mentioned, we focus on (based on the definitions byMcNutt et al. 1996 (4)) the high frequency (Type A) and low frequency (Type B) earthquakes. Type A earthquakes are (mostly) due to the shearing on faults. The slip occurring between two faults generates an energy release, creating waves on a wide range of frequencies. These events have clear P and S waves having a dominant frequency band 5 -15 Hz. Type B earthquakes however, are triggered by the fluid pressure, examples of this type are bubble formation, fluid flow into faults in very shallow depths, hydrolicfracturing. Type B events can be caused by shear or tensile forces. These events are slowly developing (emergent) with respect to Type A events. This causes a lack of clear S waves. Dominant frequency band of Type B events are 1-5 Hz.

Figure 1 :

 1 Figure 1: The Type A and Type B earthquakes obtained from seismology recordings modified after Lin et al. 2007. Top: A long-period seismogram generated by a Type B Earthquake (M w = 3.5). Bottom: High-Frequency seismogram with clear P wave and S wave arrivals due to a Type A earthquake (M w = 3.6).

4 111Figure 2 :

 42 Figure 2: Experimental setup of air injection into the Hele-Shaw cell modified after Turkayaet al. 2015 (18). A high speed camera is recording the Hele-Shaw cell from above. The accelerometers are stick on the bottom glass of the Hele-Shaw cell. The recorded signal is amplified and recorded in a computer. A signal generator is used to trigger optical and acoustic recordings at the same time to ensure synchronization.

Figure 3 :

 3 Figure 3: Optical data showing the evolution of the porous medium in a Hele-Shaw cell. A: Initial state of the porous medium before injection starts. B: Initial pressure development in solid-fluid interface prior to fracturing. C: Fluid pressure is high enough to initiate channeling in the points indicated with yellow circles. D: Channel network is being developed inside the porous medium.

Figure 4 :

 4 Figure 4: Experimentally reproduced Type A event. Top: Signal recorded during the event. Bottom: Power spectrum of the signal recorded during the event, normalized with the peak amplitude.

Figure 5 :

 5 Figure 5: Experimentally obtained Type-B event. Top: Signal recorded during the event. Bottom: Power spectrum of the signal recorded during the event, normalized with the peak amplitude.

124Chapter 3 .

 3 Investigating Solid-Fluid Interactions inside the Hele-Shaw Cell

126Chapter 4 .

 4 Analysis of the Deformation of the Porous Medium Using Optical Data

FIG. 1 .

 1 FIG.1. Sketch of the experimental setup. A) Top-down view of the prepared cell: The granular medium is confined inside the cell by three impermeable boundaries and an airpermeable boundary on the outlet side. The granular medium is placed against the outlet boundary such that it has a linear interface against a region empty of beads on the sealed inlet side, where pressurized air can be injected. B) Side view of the setup showing the high speed camera placed above. The glass plates are clamped together with aluminum framing while the cell gap is controlled with spacers, which are not shown here.

FIG. 2 .

 2 FIG. 2. Example of image processing. A: The initial granular medium, B: A snapshot during the experiment, C: The resulting binary image where the pattern has been segmented. The segmentation is done by setting all pixels having 30 % less grayscale intensity in a snapshot compared to the initial image to the value 1 while the other pixels are set to the value 0.

FIG. 3 .

 3 FIG. 3. Visual representation of the pressure ranges for the observed flow regimes, and typical final structures for each category. In the left plot: The ordinate axis indicate the category, for which the filled/thick symbols with error bars show the mean pressure within two standard deviations Pin ± 2σ. The lines with arrows just below indicate the range of injection pressures for experiments performed in the category, and the open/thin symbols indicate the pressure for the individual experiments. Open/thin symbols stacked on top of each other indicate the number of experiments done at a certain pressure. In the photos to the right: Typical final structures per category, where the granular medium is gray and areas empty of beads are black, showing that for increasing Pin the fingers become longer and thicker. A full-size figure of the patterns shown here is found in the supplementary material.

3 2 in

 32 . In figure9B, the finger tip positions per category is divided by P 3 2

in 8 DepthFIG. 5 .

 85 FIG.5. Average finger thickness w as function of depth x into the medium for the final structures, averaged per categories where fingers occur (Cat.3 -6). We see that for the structures in Cat. 5 (erosion) and Cat. 6 (erosion and finger merging) the fingers are thicker -up to the finger tip for Cat. 6.

FIG. 6 .

 6 FIG. 6. Snapshots from an experiment with Pin = 200 kPa, showing examples of erosion and finger merging. The green rectangle indicates the area cut out and enlarged in the bottom right inset: Here we see erosion inside the channel as beads (gray area) are removed over time. The red arrows points at fingers that will merge into the main channel in the following snapshot, carving out trapped clusters of beads which are then mobilized in the flow direction. A full-size version of this figure is found in the supplementary material.

9 DepthFIG. 7 .

 97 FIG.7. The average number of fingers N f as function of depth x into the granular medium, averaged over the final patterns in the same category. The plot indicates that for increasing injection pressure, an initial region with fewer fingers extends deeper into the medium before an increase in N f towards the final depth. This suggests more branching of the channels as the beads become harder to displace. The sudden drop in N f close to the tips is due to the finite size of the patterns.

FIG. 8 . 10 timeFIG. 9 .

 8109 FIG. 8. A)The finger thickness w plotted in a log-log plot as function of the distance d from the finger tip. The data shown is the average of about 50 measurements in time for each experiment, averaged again for each category. The data is collapsed by the finger thickness at the base of the tips w b found as the average finger thickness where w has reached a more or less constant level, which seems to be a typical length scale for the tip region. The finger thickness at the base of the tips is found as. The average slope at different times during several experiments is β = 0.68, indicating a more rounded shape of the tips. In the inset: Snapshot of typical finger tips together with profiles corresponding to β = 0.68 (red lines). B) The slope β averaged per category as function of time during experiments. The plot shows that the slope for the tip profile fluctuates around β = 0.68 over time during the experiments, suggesting that the tips are more rounded than cusp-like at all times.

FIG. 10 .

 10 FIG. 10. Log-log plot of the finger tip velocity v(t) as function of time, averaged per category. The data is collapsed along the ordinate axis by dividing v(t) by P 3 2

FIG. 11 .FIG. 12 .

 1112 FIG. 11. Average box-counting slopes per category, the lower and upper limits are indicated with the vertical lines and mark the typical thickness of the thinnest fingers (1 cm) up to the cell width (32 cm). In the inset: Box-counting fractal dimensions DB averaged per category, the error bars show 2 standard deviations. Patterns in categories 4 and 5 with DB = 1.53 and DB = 1.54 are within the range of earlier measured values for viscous fingers in saturated porous media (between the dashed lines). The less developed patterns in category 3 is below this range, with DB = 1.41, and the thicker category 6 patterns are close above with DB = 1.63. The solid line indicates DB = 1.71 for Diffusion Limited Aggregation patterns.

FIG. 14 .

 14 FIG. 14. Left: The collapsed finger tip positions x ′ = x(t ′ )/v0 • tc as function of t ′ = t/tc, for individual experiments with various injection pressure. The dotted vertical line indicates the cross-over at t ′ = 1, and the dashed curve shows the proposed function fx(t ′ ). Right: Log-log plot of the collapsed finger tip velocities v ′ = v(t ′ )/v0 as function of t ′ = t/tc for the same individual experiments. The dotted vertical line indicates the cross-over at t ′ = 1, and the dashed curve shows the proposed function fv(t ′ ). Recall that the parameters used for the collapse are given by the injection pressure, as v0 = C • P 3 2

3 2 )

 2 and α = -2.5.

5. 1 Introduction

 1 Noisy acoustic signal localization is a difficult problem having a wide range of applications. In this chapter, first we propose a new localization method applicable for thin plates, based on energy amplitude attenuation and inversed source amplitude comparison. This inversion is tested on synthetic data using a direct model of Lamb wave propagation and on experimental dataset (recorded with 4 Brüel & Kjae r Type 4374 miniature piezoelectric shock accelerometers, 1 -26 kHz frequency range). We compare the performance of this technique with classical source localization algorithms, arrival time localization (ATL), time reversal localization (TRL), localization based on energy amplitude (EBL).The experimental setup in figure5.1 consists of a glass / plexiglass plate of dimensions of 80 cm × 40 cm × 1 cm equipped with four accelerometers and an acquisition card. Signals are generated using a steel, glass or polyamide ball (of different sizes) hitting the plate quasi perpendicularly from a height of 2-3 cm. Signals are captured by sensors placed on the plate on different locations.

Figure 5 . 1 :

 51 Figure 5.1: Experimental setup with steel ball hitting on thin plate (glass or plexiglass). In an area covered by the sensors placed as a rectangular geometry, steel ball is launched to the several hitting locations.

Figure 5 . 2 :

 52 Figure 5.2: Incremental displacement map between the snapshots before and after an acoustic emission. The displacement magnitude map (for a time step of 8 ms between two images) is obtained by following the expression ||u|| = u 2x + u 2 y where u x,y shows displacements in x and y axes. See the section 5.3 for more details.

FIG. 1 .FIG. 2 .

 12 FIG.1. Effect of the source energy estimation grid on the resolution. The result is identically the same for 10 and 100 points in the grid.

FIG. 1 .

 1 FIG. 1. Comparison of the recorded signal with the noise added signal on glass plate.

FIG. 2 .

 2 FIG. 2.Comparison of the recorded signal with the noise added signal on plexiglass plate.

FIG. 3 .

 3 FIG.3. Effect of the signal to noise ratio (added artificially) of the signal over the estimation quality in different methods on plexiglass plate.

FIG. 1 .

 1 FIG. 1. Different sensor geometries used in numerical simulations. White area represents the plate and blue dots show the sensor positions.

Figure 1 .

 1 Figure 1. (A) The acquisition chain of the aero-fracturing experiments with a Hele-Shaw cell modified after [Turkaya et al., 2015]. The signal acquisition card, camera and the electrovalve connected to the air pump are triggered at the same time via a TTL signal sent from the signal generator to have synchronized optical and acoustic data. The sensors are placed on the bottom glass plate of the Hele-Shaw cell. (B) the image showing the accelerometers stick under the Hele-Shaw cell is given. In (C), accelerometers are sketched and numbered. Red dots show the positions of the accelerometers with indicated distances. With red numbers each accelerometer (and thus, the recorded signal) is numbered for future reference.

  relation (DIC) procedure, called Ncorr (an open source 2D digital image correlation MAT--5-

Figure 3 .

 3 Figure 3. Displacement magnitude field obtained using the images at the start and at the end of an acoustic event. The red bar shows a distance of 5 cm. The colorbar shows the displacement occurred during a time step of 8 ms

(

  and coherent) monitoring tools to detect and localize fracturing and channeling in a porous medium. Different acoustic events having different source natures (and locations) are investigated optically and acoustically. Deformation magnitude maps obtained from digital -9-

Figure 4 . 10 -Figure 5 . 11 -Figure 6 .

 4105116 Figure 4. Type of different forces creating different polarization on the recorded signals. These differences are detected based on the sign of the recorded displacements. Left: Inplane forces generating compression (or relaxation) at all sensors. Right: Shear forces generating compression and relaxation in different sides.

Figure A. 1 :

 1 Figure A.1: The preparation procedure of the experiment. (A) The accelerometers are stick to the bottom plate and steel ball spacers are placed on the sides. Then, using a cardboard a mold is prepared to give shape to the porous medium. (B) The cell is placed vertically and the porous medium is filled with polystyrene beads. (C) The cell is placed horizontally and the cardboard is removed. (D) Finally, the spacers are removed and the top plate stays on the porous medium.

Figure A. 4 :

 4 Figure A.4: Comparison of experimental (green) and theoretical (red, magenta) dispersion curves.

  1.1 Problématique et Enjeux . . . . . . . . . . . . . . . . . . . . . 2 1.2 La Fracturation dans un Milieu Poreux . . . . . . . . . . . . 2 1.3 Localisation des Emissions Acoustiques . . . . . . . . . . . .

5 1.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5 Sommaire des Articles . . . . . . . . . . . . . . . . . . . . . . 8 1.5.1 L'Article: Aerofracture through a double looking glass, mixing optics and acoustics . . . . . . . . . . . . . . . . . . . . . . . 8 1.5.2 L'Article: Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium . . . . . . 11 1.5.3 L'Article: Numerical Studies of the Acoustic Emissions during Pneumatic Fracturing . . . . . . . . . . . . . . . . . . . . . . 13 1.5.4 L'Article: Pneumatic fractures in confined granular media . . 17

  Introduction Contents 2.1 Brief History of Research on Porous Medium . . . . . . . . 21 2.2 Brief History of the Research on Acoustic Emissions . . . . 26 2.2.1 Acoustic Instrumentation . . . . . . . . . . . . . . . . . . . . 26

	Chapter 2

2.3 Brief History of the Signal Localization . . . . . . . . . . . . 29 2.3.1 Time Difference of Arrival . . . . . . . . . . . . . . . . . . . . 30 2.3.2 Time Reversal Localization . . . . . . . . . . . . . . . . . . . 30 2.3.3 Energy Based Localization . . . . . . . . . . . . . . . . . . . 31 2.4 Motivations and Objectives . . . . . . . . . . . . . . . . . . . 32 20 Chapter 2. Introduction

  1, 2 , Megan Zecevic 3 , Guillaume Daniel 3 , Eirik G. Flekkøy 2 and Knut J.Måløy 2

TABLE I .

 I The flow categories with related pressure values; mean overpressure within two standard deviations Pin ± 2σ, as well as the number of experiments N analyzed per category.

	Category	Pin ± 2σ [kPa]	N
	1	9 ± 5	3
	2	26 ± 18	5
	3	64 ± 50	5
	4	110 ± 44	6
	5	134 ± 43	10
	6	213 ± 50	4

  FIG. 4.The ratio of new growth dA to the cell width W is plotted as function of distance d behind the longest finger tip for typical experiments in categories 3 to 6.

															7
		1							1
		0.9			Cat. 3		0 -10 % 10 -20 %		0.9			Cat. 4		0 -10 % 10 -20 %
							20 -30 %							20 -30 %
		0.8					30 -40 %		0.8					30 -40 %
							40 -50 %							40 -50 %
	New growth, dA / W	0.4 0.5 0.6 0.7					50 -60 % 60 -70 % 70 -80 % 80 -90 % 90 -100 %	New growth, dA / W	0.4 0.5 0.6 0.7					50 -60 % 60 -70 % 70 -80 % 80 -90 % 90 -100 %
		0.3							0.3
		0.2							0.2
		0.1							0.1
		0							0
		0	10	20	30	40	50	60	0	10	20	30	40	50	60
				Distance behind channel tip, d [cm]					Distance behind channel tip, d [cm]
		1							1
		0.9			Cat. 5		0 -10 % 10 -20 %		0.9			Cat. 6		0 -10 % 10 -20 %
							20 -30 %							20 -30 %
		0.8					30 -40 %		0.8					30 -40 %
							40 -50 %							40 -50 %
	New growth, dA / W	0.4 0.5 0.6 0.7					50 -60 % 60 -70 % 70 -80 % 80 -90 % 90 -100 %	New growth, dA / W	0.4 0.5 0.6 0.7					50 -60 % 60 -70 % 70 -80 % 80 -90 % 90 -100 %
		0.3							0.3
		0.2							0.2
		0.1							0.1
		0							0
		0	10	20	30	40	50	60	0	10	20	30	40	50	60
				Distance behind channel tip, d [cm]					Distance behind channel tip, d [cm]

  It has been established that flow in porous media is better described by another Laplacian model, i.e. the Dielectric Breakdown Model 12 FIG. 13. Left: Local fractal dimensions (DL + 1) as function of depth into the medium for the final structures averaged per category. We see that the patterns in all categories have more or less constant and well defined local dimensions as functions of depth. The main parts of the patterns in categories 3-5 have local dimensions corresponding to established values for viscous fingers in saturated porous media (between the dashed lines). For patterns in category 6, the local dimensions correspond to values closer to the fractal dimension for DLA patterns and flow patterns in a saturated Hele-Shaw cell (solid line). Right: Examples of local box-counting slopes at different depths of the same experiment.

											2.5
			2												x=10 cm x=20 cm
															x=30 cm
	L + 1 Local fractal dimension, D	1.8 1.2 1.4 1.6	Cat. 3						log 10 ( Box count )	2 0.5 1 1.5				x=40 cm Slope at 10 cm Slope at 20 cm Slope at 30 cm Slope at 40 cm
				Cat. 4
			1	Cat. 5 Cat. 6							0
			0	10	20	30	40	50	60	70	-1	-0.5	0	0.5	1	1.5
					Depth into granular medium, x [cm]					log 10 ( Box size ) [cm]
			1.6
			1.4
			1.2
	x'	0.8 1							log 10 v'
			0.6
									50 kPa
			0.4						100 kPa 150 kPa
									200 kPa
			0.2						250 kPa
									f x (t/t c )
			0
			0		5		10			15	-1	-0.5	0	0.5	1	1.5
						t'							log 10 t'

  The final channel length x f is plotted as function of Pin -Pcat,1, where Pcat,1 = 9 kPa is the average injection pressure for category 1, where x f = 0. We see that there is a rather scattered (due to randomness in the initial packings), but increasing trend in x f for increasing Pin, but also that the final channel length for the highest injection pressures is limited by the closed outer boundary. The fit (Pin-Pcat,1) 0.68 was found from the averaged x f per pressure for Pin ≤ 150 kPa and suggest a power-law relationship. The average limiting line xmax,avg was found by averaging x f over the experiments with Pin ≥ 200 kPa, to indicate the influence of the cell boundary. Hopefully, future experiments with much more repetitions could be conducted to reveal a clearer relationship between x f and Pin.

		70
		60
	f [cm] Final channel length, x	20 30 40 50			Individual experiments Average per pressure (P in -P Cat,1 ) 0.68
		10			x max,avg
					Cell boundary
		0
		0	50	100	150	200	250
				Injection pressure, P in -P Cat,1 [kPa]
	FIG. 15.

  Contents 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 5.2 Article: Localization Based On Estimated Source Energy Homogeneity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 5.3 Draft Article: Source Localization of Acoustic Emissions during Pneumatic Fracturing . . . . . . . . . . . . . . . . . . 166 5.4 Conclusion and Future Work . . . . . . . . . . . . . . . . . . 182 146 Chapter 5. Finding the Source of the Acoustic Emissions

Localization based on estimated source energy homogeneity

  This study is mainly conducted by the main author under the supervision of Renaud Toussaint. Then, to develop the work further, we discussed several times with all coauthors. All co-authors have contributed with some suggestions on references, scientific methods on analysis and comparison, and presentation of the obtained results. This article was published in Review of Scientific Instruments in September 2016.Semih Turkaya,1,a) Renaud Toussaint, 1 Fredrik Kvalheim Eriksen, 1,2 Olivier Lengliné, 1 Guillaume Daniel,3 Eirik G. Flekkøy, 2 and Knut JørgenMåløy 2 Acoustic signal localization is a complex problem with a wide range of industrial and academic applications. Herein, we propose a localization method based on energy attenuation and inverted source amplitude comparison (termed estimated source energy homogeneity, or ESEH). This inversion is tested on both synthetic (numerical) data using a Lamb wave propagation model and experimental 2D plate data (recorded with 4 accelerometers sensitive up to 26 kHz). We compare the performance of this technique with classic source localization algorithms: arrival time localization, time reversal localization, and localization based on energy amplitude. Our technique is highly versatile and out-performs the conventional techniques in terms of error minimization and cost (both computational and financial). Published by AIP Publishing. [http://dx.doi.org/10.1063/1.4962407]

	5.2. Article: Localization Based On Estimated Source Energy
	Homogeneity	149
	5.2 Article: Localization Based On Estimated Source
	Energy Homogeneity
	3 Magnitude, Sainte Tulle, France
	(Received 7 June 2016; accepted 25 August 2016; published online 8 September 2016)

REVIEW OF SCIENTIFIC INSTRUMENTS 87, 096101 (2016) Note: 1 IPG Strasbourg, CNRS, Université de Strasbourg, Strasbourg, France 2 Department of Physics, University of Oslo, Oslo, Norway

  2 . The minimum of this standard deviation, σ(r s ) over the plate will

	096101-2	Turkaya et al.		Rev. Sci. Instrum. 87, 096101 (2016)
	0034-6748/2016/87(9)/096101/3/$30.00	87, 096101-1	Published by AIP Publishing.

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Download to IP: 77.154.224.126 On: Thu, 08 Sep 2016 14:00:26

TABLE I .

 I Table showing the average distance between the real source position and the estimated position (i.e., error) in mm for different cases in the signal generation.

	Mean error, glass plate (mm)	ATL	TRL	EBL	ESEH
	Infinite, non-dispersive	0	0	7	0
	Infinite, dispersive	45	26	35	4
	Finite, non-dispersive	2	2	20	5
	Finite, dispersive	35	27	50	36
	Experimental	21	29	48	33
	Mean error, plexiglass plate (mm)	ATL	TRL	EBL	ESEH
	Infinite, non-dispersive	0	0	7	0
	Infinite, dispersive	30	61	27	11
	Finite, non-dispersive	0	0	21	3
	Finite, dispersive	36	37	48	40
	Experimental	48	34	34	30

TABLE II .

 II Table showing the computational time in milliseconds for different cases.

	Mean computational cost (ms)	ATL	TRL	EBL	ESEH
	Infinite, non-dispersive	6	31	5567	125
	Infinite, dispersive	6	31	5567	125
	Finite, non-dispersive	6	43	6983	257
	Finite, dispersive	6	22	4760	90

  FIG. 1. Figure showing the decay of the average error on plexiglass plate with different sampling rates in different methods. Signals with lower sampling rates are obtained by decimating the experimentally recorded signal with 1 MHz sampling rate. Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Download to IP: 77.154.224.126 On: Thu, 08 Sep 2016 14:00:26 FIG. 2. Figure showing the effect of the SNR (added artificially) of the signal over the estimation quality in different methods on glass plate.

	096101-3	Turkaya et al.	Rev. Sci. Instrum. 87, 096101 (2016)

TABLE I .

 I Effect of the array geometry on the localization error. Draft Article: Source Localization of Acoustic Emissions during Pneumatic FracturingThis study is mainly conducted by the main author under the supervision of Renaud Toussaint. The results obtained from optical dataset are done by Fredrik K. Eriksen. Then, to develop the work further, we discussed several times with all coauthors. All co-authors have contributed with some suggestions on references, scientific methods on analysis and comparison, and presentation of the obtained results. This draft article is going to be submitted to Journal of Geophysical Research: Solid Earth in early 2017. Please consult the latest version. This draft article is going to be submitted to Geophysical Research Letters in late 2016.

	5.3
	Geometry / Error (mm) ATL TRL EBL ESEH
	Square	36 37 48 40
	Top-Left	26 52 39 23
	Top-Right	34 43 49 32
	Bottom-Left	35 39 51 30
	Bottom-Right	30 48 38 23
	1
		165

Institut de Physique du Globe de Strasbourg, UMR7516, CNRS, University of Strasbourg/EOST, 5 rue Descartes, 67000 Strasbourg, France

Physics Department, University of Oslo, Oslo, Norway

Centre National de la Recherche Scientifique, Institut de Physique du Globe de Strasbourg, Université de Strasbourg, Strasbourg, France,

Department of Physics, University of Oslo, Oslo, Norway,

Magnitude, Sainte Tulle, France

Frontiers in Physics | www.frontiersin.org 2 September 2015 | Volume 3 | Article 70

Frontiers in Physics | www.frontiersin.org

September 2015 | Volume 3 | Article 70

Frontiers in Physics | www.frontiersin.org 4 September 2015 | Volume 3 | Article 70

in , while at a critical time tc there is a cross-over to the behavior v(t) ∝ t -α where α is close to 2.5 for all experiments. Finally, we explore the fractal dimension of the fully developed patterns. For example for patterns resulting from intermediate Pin around 100-150 kPa, we find that the box-counting dimensions lie within the range DB ∈[1.53, 1.62], similar to viscous fingering fractals in porous media.

Acknowledgments

Acknowledgments

We would like to thank Alain Steyer and Miloud Talib for their technical support during experimental campaigns. We would also like to thank Maxime Farin for stimulating discussions during this work. This project has received funding from the European Union's Seventh Framework Programme for research, technological development and demonstration under grant agreement no 316889, from the REALISE program of the Alsatian research network, from the Universities of Oslo and Strasbourg via a gjesteforsker program and an IDEX Espoirs award. Acknowledgments We would like to thank Alain Steyer and Miloud Talib for the technical support and Maxime Farin for discussions during this work. This project has received funding from the European Union's Seventh Framework Programme for research under grant agreement no 316889, from the REALISE program of the Alsatian research network, from the Universities of Oslo and Strasbourg via a gjesteforsker program and an IDEX Espoirs award. ACKNOWLEDGMENTS This project has received funding from the European Unions Seventh Framework Programme for research, technological development and demonstration under grant agreement no. 316889, ITN FlowTrans. We thank Alain Steyer and Miloud Talib for their technical support with experimental equipment. We thank A. Steyer and M. Talib for the technical support and M. Farin for discussions during this work. This project has received funding from the EU's 7th Framework Programme for research under Grant Agreement No. 316889, from the REALISE program of the Alsatian research network, from the Universities of Oslo and Strasbourg via a gjesteforsker program, and an IDEX Espoirs award. We also thank J. Farquharson for his help in language to make the article easier to understand. Acknowledgments We would like to thank Alain Steyer and Miloud Talib for the technical support. This project has received funding from the European Union's Seventh Framework Programme for research under grant agreement no 316889, from the REALISE program of the Alsatian research network, from the Universities of Oslo and Strasbourg via a gjesteforsker program and an IDEX Espoirs award.

Source Localization of Acoustic Emissions during Pneumatic

Fracturing

Conflict of Interest Statement:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Abstract

Fluid injection into a deformable porous medium is a challenging area of physics having a wide range of application. The food, construction and petroleum industries can be given as example areas of application. In this research, we investigate experimentally and numerically, in a Hele-Shaw cell, pneumatic fracturing the porous medium. In the experimental phase, we inject air inside the porous medium at a loose state to create compaction, channeling and fracturing while monitoring the cell with 4 accelerometers and a high speed camera. Furthermore, we developed a numerical model in two steps (1) poroelasto-plasticity based model to explain dynamic fluid pressure variations (2) solid stress model based on Janssen's silo theory. The contributions of the different pressure sources (air pressure in channels, solid stress etc.) are compared with respect to amplitude and frequency. Afterwards, the pressure variations are convolved with a Lamb Wave green function to generate acoustic emissions numerically. The physics behind the experimentally observed corner frequency evolving into a power law decay in the power spectra are explained using numerical models. The frequency bands (in the power spectrum) influenced by the size of the carved channels and the Hele-Shaw cell are computed are found to be compatible with the peaks of the acoustic emissions.

emissions in a porous medium, Frontiers in Physics, 3, 70. Turkaya, S., R. Toussaint, F. K. Eriksen, O. Lengliné, G. Daniel, E. G. Flekkøy, and K. J. Måløy (2016), Note: Localization based on estimated source energy homogeneity, Review of Scientific Instruments, 87 (9), 096,101. Valkó, P., and M. Economides (1995), Hydraulic fracture mechanics, Wiley.

Van der Hoef, M., M. Ye, M. van Sint Annaland, A. Andrews, S. Sundaresan, and J. Kuipers (2006), Multiscale modeling of gas-fluidized beds, Advances in chemical engineering, 31, 65-149.

Vinningland, J. L., Ø. Johnsen, E. G. Flekkøy, R. Toussaint, and K. J. Måløy (2007a), Granular rayleigh-taylor instability: Experiments and simulations, Physical Review Letters, 99, 048,001, doi:10.1103/PhysRevLett.99.048001. Vinningland, J. L., Ø. Johnsen, E. G. Flekkøy, R. Toussaint, andK. J. Måløy (2007b), Experiments and simulations of a gravitational granular flow instability, Physical Review Nonlinear,and Soft Matter Physics,76(5).

Vinningland, J. L., Ø. Johnsen, E. G. Flekkøy, R. Toussaint, and K. J. Måløy (2010), Size invariance of the granular rayleigh-taylor instability, Physical Review E -Statistical, Nonlinear,and Soft Matter Physics,81(4).

-31-Explaining type A (shear) and B (low frequency) earthquakes with lab scale experiments and numerical simulations Semih Turkaya teristics. The events are detected using a Short Time Average / Long Time Averaging (23)(24)(25)(26)(27).

When this ratio passes a predefined threshold, the subsignal is considered an event. After a pre-treatment of the signal using a butterworth filter having different frequency spans for high frequency and low frequency events. For low frequency events we used a detection criterion in the frequency span 0.1 -10 kHz and for high frequency events we changed this span to 100 -200 kHz (18,19).

High Frequency Acoustic Emissions (Type-A)

After the channel formation is stopped, there is a pressure diffusion phase which continues inside the plate. In a sub-group of grains susceptible to move eventually due to pore pressure the stability equation can be expressed as ( 28):

in equation 1, Id is the identity matrix, σ T is the total stress of the continuum, the solid stress is σ s , and the pore pressure p. The Coulomb condition (i.e. σ s s ≥ µσ s n where σ s n is the norm of the normal stress and σ s s is the norm of the shear (frictional) solid stress in this subgroup of grains. µ is the friction coefficient representing the characteristics of the contacting material solid phase and glass plates) (29). leading to develop stick-slip acoustic emissions combined with the pore fluid pressure in the Hele-Shaw cell can be derived similar to the Terzaghi's 1936 effective stress formulation (30,31): We perform experiments where air is injected at a constant overpressure Pin, ranging from 5 to 250 kPa, into a dry granular medium confined within a horizontal linear Hele-Shaw cell. The setup allows us to explore compacted configurations by preventing decompaction at the outer boundary, i.e. the cell outlet has a semi-permeable filter such that beads are stopped while air can pass. We study the emerging patterns and dynamic growth of channels in the granular media due to fluid flow, by analyzing images captured with a high speed camera (1000 images/s). We identify 6 typical flow regimes, depending on the imposed overpressure, ranging from no channel formation at low Pin around 10 kPa, to large thick channels formed by erosion and fingers merging at high Pin around 200 kPa. The flow regimes where channels form are characterized by typical finger thickness, final depth into the medium and growth dynamics. The shape of the finger tips during growth is studied by looking at the finger width w as function of distance d from the tip. The tip profile is found to follow w(d) ∝ d β , where β = 0.68 is a typical value for all experiments, also over time. This indicates a singularity in the curvature d 2 d/dw 2 ∼ κ ∼ d 1-2β , but not of the slope dw/dd ∼ d β-1 , i.e. more rounded tips rather than pointy cusps, as they would be for the case β > 1. For increasing Pin, the channels generally grow faster and deeper into the medium. We show that the channel length along the flow direction has a linear growth with time initially, followed by a power law decay of growth velocity with time as the channel approaches its final length. A closer look reveal that the initial growth velocity v0 is found to scale with injection pressure as v0 ∝ P

SUPPLEMENTARY MATERIAL S1 -LOCALIZATION METHODS

A. Arrival Time Delay Localization (ATL)

Arrival time delay localization is based on the difference in arrival time occurring when the signal is propagating from one source to the different receivers at different positions 1-6 . This travel time is defined by the absolute time passing between the start of the experiment and the first arrival of the signal computed via thresholding the amplitude of the signal by 3 to 6 times (depending on the signal to noise ratio) of its standard deviation. This method requires the wave velocity to be known a-priori to calculate the position which is fitting best to the acquired delay between the time at each sensor and the earliest time for all sensors. A grid search (i.e. an exhaustive search for all the possible positions over the plate or the medium) is conducted comparing the modeled arrival time delay with the delay obtained from the real data. The delay is used instead of absolute travel time since the time of causing event (e.g. rupture/hit, slip etc.) is unknown. The position corresponding to the minimum of the square root of the average of the squared difference (RMS) between these observed and modeled travel time difference gives the estimated source position. A simple mathematical explanation for this method is given below:

where n and N represent the sensor considered and the maximum number of sensors (taken as 4 in this study) respectively. d d (n) is the delay between sensors obtained by subtracting the minimum arrival time from arrival time set in real data, t m (n, x, y) is the time of flight modeled with the pre-defined (or approximated) wave velocity for a specific grid point (x, y) and γ m (x, y) is the minimum time of flight in the sensor array for that grid point. Eventually, the argument of the minimum of RMS, min

gives the estimated position of the source (x s , y s ).

B. Time Reversal Localization (TRL)

Time reversal localization 7-9 is based on sending signals back in time to a grid point, then superposing them to find the maximum possible amplitude for each grid point. For A. Non-Dispersive Signals

A finite element scheme is made to solve wave equation 1 numerically. To simplify the 2D equation into 1D, polar coordinates are used. In this equation c is the wave velocity and ξ is the displacement of the plate.

B. Dispersive Signals

In this equation 2,3 , ξ(R) is the displacement in fourier domain in the distance R from the source, k is the wavenumber having a dispersion relation

where E is the Young's modulus, ν is the Poisson's ratio, h is the thickness of the plate. Furthermore, for the frequencies greater than the cutoff frequency (which satisfies kh ≃ 1 dispersion relation) will be changed into non-dispersive one (i.e. ω k = dω dk = c) where c is the wave velocity. However, non-dispersive component of the signal corresponds to only ≈ 10% of the total energy which is negligible.

SUPPLEMENTARY MATERIAL S3 -SPATIAL RESOLUTION

The overall precision may increase with increasing spatial resolution of the grid. However, as the spatial resolution increases 10 times the computational time increases 70 times. In the figures 1 and 2 the grid sizes 10 mm and 1 mm are compared. We can see that, here for this setup, increasing the spatial resolution does not have a significant impact on the estimation results. Colorbar in each subfigure corresponds to the normalized results of the functions described in the corresponding sections. Excluding the time reversal method where the maximum amplitude is desired, all of the methods are looking for the minimum value in their corresponding grid search results.

Introduction

Acoustic signal localization is applied in many different areas of science [Gershman et al., 1995;Valin et al., 2003;Elnahrawy et al., 2004;Malioutov et al., 2005;Zhu et al., 2007;Fink, 2015;[START_REF] Garnier | Super-resolution in time-reversal focusing on a moving source[END_REF]. Locating the source of a received signal is an essential task in many different branches of science. In robotics, speech-source tracking is done to automate cameras to follow the speaker [Brandstein et al., 1997]. In electronics, touchscreens are a very popular example of a signal localization. It is necessary for the system to locate the touch of a user to transmit information to process a simple message [START_REF] Terlizzi | Multi-touch display screen with localized tactile feedback[END_REF]. Finding the epicenter of an earthquake is necessary to understand how it was generated, where exactly the movement came from, which region is more risky for the aftershocks [START_REF] Aki | Quantitative Seismology[END_REF]. Recently, some of the signal localization types are reviewed in the work [Turkaya et al., 2016]. Here in this article, we

show the application of Estimation of Source Energy Homogeneity (ESEH) to locate the source of the acoustic emissions during aerofracturing experiments [Turkaya et al., 2015].

After this application, these localization results are compared with the Digital Image Correlation (details can be found in [Eriksen et al., 2016]) results obtained from the optical recordings via a high speed camera during the experiment.

-2-

Experiments

The aerofracturing experiments here are done in a Hele-Shaw cell. As explained in a previous article [Turkaya et al., 2015], we are using two glass plates 80 cm × 40 cm × 1 cm with an aperture of 1 mm. Experimental setup is shown in figure 1 A). The acquisition chains of optical and acoustic data are presented. The system is triggered via a signal generator to have synchronization between optical and acoustic data. For optical data we recorded with 125 images per second and for acoustic data we used 1 MHz sampling rate. Three boundaries of this cell is sealed with a double sided tape, while the fourth side is covered with a mesh filter which lets the air pass to the atmosphere while keeping the grains inside of the cell. In this particular experiment, we used 1 bar overpressure from the inlet and we used 80 µm grains having material density 1.005 g/cm 3 in loose state inside the cell. The cell is filled vertically before sealing with a semi-permeable mesh.

Then, the cell is rotated vertically until a homogenous loose packing is reached. Finally, the cell is placed horizontally without changing this loose state. An empty space is provided between the air inlet and the solid-air interface to provide a homogenous pressure over the width of the cell.

During this experiment, acoustic signals are recorded using 4 accelerometers attached to the bottom plate (see Figure 1B). The sensors are placed towards the outlet (where most of the stick slip events are happening) to have a better signal resolution.

Event Detection

The acoustic events are detected using a well known Short Time Average / Long Term Average (STA/LTA) algorithm, which is particularly robust in complex cases like fracturing if calibrated properly [Allen, 1978;Baer and Kradolfer, 1987;Earle and Shearer, 1994;Wong et al., 2009;Trnkoczy, 1998]. This method is used in [Turkaya et al., 2015] as:

where N s,l are the size of the short (0.05 ms) and long (1 ms) time windows respectively, s(i) is the value of the signal (pre-filtered if necessary) at time i, and τ(i) is the ratio of these two averages. If τ(i), passes a threshold (τ th = 4.3) defined a priori, algorithm considers that short time window as an event and records it. For event detection we used the The detected event in the red window is presented. The experiments done in the Hele-Shaw cell presented in Chapter 4.1 was completely confined for grains and the plates were clamped on the corners to take into account the cases where only the porous medium are moving. To take into account the cases where the plate may move under the influence of fluid flow, we modified the experimental setup as follows:

1. The corners are entirely permeable, both air and porous medium can flow outside the plate freely.

2. The air is injected in the middle of the porous medium to have axisymmetric deformation on the plate.

3. Glass plates are replaced with plexiglass plates due to their light and deformable nature.

Using open boundaries on the Hele-Shaw cell makes the preparation significantly difficult. Since the boundaries are open, it is not possible to make a decompaction phase after filling to obtain homogenous loose solid fraction in the porous medium due to leaking problem. However, by using a rectangular filling shape we tried to give a preferred way of deformation to the porous medium. We filled simply can deform parallel to the shorter edge which have the same (approximately) solid fraction along this axis. The protocol we followed to prepare the Hele-Shaw cell for the experiment is presented in figure 5.1.

As seen in figure A.1A, first the acoustic sensors are placed on the bottom plate (80 cm × 40 cm × 1 cm). We are using 4 Brüel & Kjaer Type 4374 miniature piezoelectric shock accelerometers (up to 26 kHz frequency range) and 4 micro 80 PAC piezoelectric transducers (200 kHz to 900 kHz frequency range) close to the corners of the granular fill. Then, we put the plate horizontally (sensors are stick below the plate) and we place the steel ball placers (0.625 cm diameter) between on the bottom plate where they are restrained in one dimension. A steel ball can only move parallel to the closest plate edge. When the placers are rigid enough, a "U" shaped cardboard is placed between two plates to use as a mold for the porous medium. (B) Then we place the Hele-Shaw cell vertically (the opening is on top) and fill the opening up to the limit of the cardboard mold. (C) Afterwards, we place the cell horizontally and place the air pump into the inlet, remove the cardboard and steel placers so that the top plate only rests on the granular medium. (D) The air injection is controlled via an electro valve. The air pressure is constant when the valve is open (like a step function). The whole system is triggered together as explained in figure A.2 below:

The initial noise is recorded via two types (mentioned above) of acoustic sensors with 2.5 MHz sampling rate (1.25 MHz Nyquist Frequency) and amplified. After the pre-defined noise recording period (e.g. 0.5 sec) as recording continues the electrovalve is triggered (to let the air in) along with the high speed camera (to

A.2.2 The Effect of Shear Forces

To simulate the acoustic emissions inside the Hele-Shaw cell, we combined the fluid pressure computed using the poroelasticity with the solid stress computations using Janssen's theory [START_REF] Janssen | Versuche über getreidedruck in silozellen[END_REF]]. In this work, we only took into account the variation of the normal forces on the plate, but not the tangential forces. Since the receivers that are used are recording all the normal motion (both due to normal and shear forces) on the Hele-Shaw cell, the contribution of these forces are not separately known. Considering these vibrations are purely elastic, it is possible to say that the normal vibrations are due to shear forces are equal to the shear vibrations that due to the normal forces. Therefore, using the theory of elasticity it is possible to compare the effect of normal and shear forces by recording normal and shear vibrations due to a normal force.

In figure A.5, it is possible to see normal (blue) and shear (green) vibrations due to a steel ball hitting vertically on the plate (normal force). The ratio of the amplitude of normal to shear motions are more than 3 times. Furthermore, the signature of these two motions on the power spectrum are very similar. Based on these experimental findings, we found that neglecting tangential forces for acoustic emission simulations would not significantly effect our results.