Professor Jari Kaipio
email: j.kaipio@auckland.ac.nz

Professor Daniela Calvetti
email: daniela.calvetti@case.edu

Professor Stuart Bolton
email: bolton@purdue.edu

Professor Nuutti Hyvönen
email: nuutti.hyvonen@aalto.fi

Maître Jean-Daniel Chazot
email: jean-daniel.chazot@utc.friv

Matti Niskanen

BAYESIAN INVERSION METHODS FOR ACOUSTICAL CHARACTERISATION OF POROVISCOELASTIC MEDIA

Keywords: Bayesian methods, Poroelastic materials, Parameter estimation

Natural and artificial porous materials are encountered and used in many fields of engineering, medicine and physics. One of the current challenges is to reliably characterize the pore structure of these materials as it greatly affects their physical properties. Characterization is usually performed using deterministic inversion methods that are based on minimizing a cost function in the least squares sense. These approaches can provide reliable parameter estimates when the number of parameters is small, but with an increasing number of parameters the least squares estimates and, in particular, error estimates often become unreliable.

In this thesis, we study the characterization of porous materials in the Bayesian framework for inverse problems, where uncertainties can be treated and quantified naturally. First, we apply this approach to the case of air-saturated porous materials whose frame can be modelled as rigid. The case of porous materials

Titre : METHODES D'INVERSION BAYESIENNE POUR LA CARACTERISATION ACOUSTIQUE DES MATERIAUX POROVISCOELASTIQUES

Mots clés : Méthodes bayesiennes, Matériaux poroélastiques, Estimation de paramètres Résumé : Les matériaux poreux naturels et artificiels sont rencontrés et utilisés dans de nombreux domaines de l'ingénierie, de la médecine et de la physique. L'un des défis actuels consiste à caractériser de manière fiable la structure des pores de ces matériaux car elle influe grandement sur leurs propriétés physiques. La caractérisation est généralement effectuée à l'aide de méthodes d'inversion déterministe basées sur la minimisation d'une fonction coût par les moindres carrés. Ces approches fournissent des informations relativement fiables lorsque le nombre de paramètres recherché est faible. Lorsque ce nombre devient important, ces approches montrent une limite qui se traduit par une mauvaise estimation des paramètres de minimisation par les moindres carrés.

Dans cette thèse, la caractérisation des matériaux poreux est abordée par une approche bayésienne qui permet de pouvoir considérer naturellement et quantitativement les incertitudes sur les paramètres recherchées. Dans un premier temps, cette approche est appliquée au cas de matériaux poreux à squelette rigide et saturés d'air. L'application au cas de matériaux poreux saturés par un fluide lourd est abordée dans une deuxième partie. Les méthodes ont été testées sur des mesures, effectuées dans le cadre de ce travail de thèse, en particulier au tube à impédance ou dans une cuve à eau. Les méthodes et résultats développés dans ce travail mettent en évidence que les techniques d'inversion bayesiennes permettent de caractériser de manière fiable les matériaux poro-visco-élastiques

R ÉSUM É

Les matériaux poreux naturels et artificiels sont rencontrés et utilisés dans de nombreux domaines de l'ingénierie, de la médecine et de la physique. L'un des défis actuel consiste à caractériser de manière fiable la structure des pores de ces matériaux car elle influe grandement sur leurs propriétés physiques. Dans ce contexte, la caractérisation des matériaux poreux à l'aide de méthodes acoustiques est devenue une procédure courante dans de nombreux laboratoires. Cela se justifie en particulier en raison du caractère non destructif de ces méthodes associé à un faible coût de l'équipement nécessaire. La caractérisation est généralement effectuée à l'aide de méthodes d'inversion déterministe basées sur la minimisation d'une fonction coût par les moindres carrés. Ces approches fournissent des informations relativement fiables lorsque le nombre de paramètres recherché est faible. Lorsque ce nombre devient important, ces approches montrent une limite qui se traduit par une mauvaise estimation des paramètres de minimisation par les moindres carrés. Dans cette thèse, la caractérisation des matériaux poreux est abordée par une approche bayésienne du problème inverse. Cette approche présente l'intérêt de pouvoir considérer quantitativement les incertitudes sur les paramètres recherchées. Dans les approches bayésiennes proposées, l'estimation des paramètres et l'évaluation des incertitudes sont obtenues en utilisant des algorithmes de Monte-Carlo avec des chaînes de Markov. Dans un premier temps, cette approche est appliquée au cas de matériaux poreux saturés d'air. Dans ce cas, la modélisation est réalisée en supposant le squelette comme étant rigide. Des expériences sont ensuite réalisées dans un tube d'impédance à des fréquences audibles. Les mesures ainsi réalisées permettent de tester cette première application. L'application au cas de matériaux poreux saturés par un fluide lourd est abordée dans une deuxième partie où la modélisation permet de considérer l'élasticité du squelette et son effet sur la propagation. De nouvelles mesures sont réalisées dans une cuve à eau à des fréquences ultrasonores afin de tester cette deuxième application de l'approche bayésienne. Le développement de méthodes d'inversion permettant de caractériser de manière fiable ce type de matériaux poro-visco-élastiques constitue une contribution majeure de ces travaux de thèses. Les résultats de cette thèse indiquent que l'approche bayésienne pour la caractérisation des matériaux poreux fournie des estimations robustes pour les paramètres des modèles tout en y associant leur incertitude.

Universal Decimal Classification: 519.217.2, 534.6, 539.217, 539.3, 620.179.16 INSCPEC Theasaurus: porous materials; porosity; viscoelasticity; acoustic analysis; acoustic measurement; ultrasonic measurement; acoustic impedance; nondestructive testing; inverse problems; parameter estimation; Bayes methods; Markov processes; Monte Carlo methods Yleinen suomalainen ontologia: materiaalit; materiaalitutkimus; huokoisuus; joustavuus; kimmoisuus; akustinen analyysi; ultraäänitutkimus; rikkomaton aineenkoetus; inversioongelmat; parametrit; estimointi; bayesilainen menetelmä; Markovin ketjut; Monte Carlo -menetelmät vi ACKNOWLEDGEMENTS This thesis was carried out as a jointly supervised thesis between the Department of Applied Physics in the University of Eastern Finland and the Laboratoire d'Acoustique de l'Université du Mans in the Le Mans University during the years 2015-2019. The study received financial support from the Academy of Finland, Jenny and Antti Wihuri Foundation, Saastamoinen Foundation, the RFI Le Mans Acoustique, and the DENORMS COST action. I would first like to thank my supervisor in Finland, Docent Timo Lähivaara, for your guidance, encouragement, and help over these years. I would also like to thank Docent Tomi Huttunen for guiding me in the beginning of the project. It has been a privilege to work with both of you. In addition, I am grateful for the patience and understanding you showed as my plans were changing and I went on a rather long research visit to New Zealand. I wish to express my deepest gratitude to my supervisors in France, Professor Olivier Dazel, HDR Jean-Philippe Groby, and MCF Aroune Duclos. You made me feel very welcome to Le Mans, and your enthusiasm towards acoustics and porous materials provided an exciting atmosphere to work in. During my 1.5 years in Le Mans I learned a lot about acoustics, but more importantly about the French way of life, and I will never forget the raclette nights we shared together. I would also like to express my sincerest thanks to my New Zealand supervisor Professor Jari Kaipio, who welcomed me into his group in the University of Auckland when life was taking me to the other side of the world. My research visit there turned into almost two years, during which you were always there to offer me help, guidance, encouragement, and discouragement, in the correct proportions, when it was needed. It truly was an honour to work with you and your group.

In addition, I would like to thank the official pre-examiners Professor Daniela Calvetti and Professor J. Stuart Bolton for taking the time to review my thesis and providing excellent feedback.

During these four years I have had the pleasure to live and work in three different countries. I have met lots of amazing people, but moving around has meant I have had to say a lot of goodbyes as well. First I would like to thank my good friends Elie and Valentin for always being willing to act as my French-English translators, even when we had just barely met. When I moved to France I did not speak a word of French, but you helped me navigate through the French administration and most importantly, showed me which ones of the hundreds of different cheeses were worth trying (all of them). Thanks for all the gigs and all the evenings at Le Lézard. Also thanks to Charlotte, Mathieu, Quentin, Nico, Dimitri, Weichun, and many others for making my stay in Le Mans unforgettable. I want to also thank my Kuopio friends Antti, Timo, Petteri, Arttu, Johanna, Simo, and Pekka, for the friendship, housing assistance, and discgolf company. Further, I want to thank all my Auckland office mates, Ru, Owen, Pascal, Hwan, Susana, Vincent, Scarlett, and Annabelle, for the laughs and evenings we spent together. You were always down (up?) for some bananagrams to liven up the work day.

Of course, the longest lasting support for my scientific endeavours comes from my family. I owe a huge thanks to my parents Sirpa and Jorma, and my brothers Timo and Tuomas for always encouraging and supporting me.

LIST OF PUBLICATIONS

This thesis consists of the present review of the author's work in the field of acoustics and the following selection of the author's publications: I M. Niskanen, J.-P. Groby, A. Duclos, O. Dazel, J. Le Roux, N. Poulain, T. Huttunen, and T. Lähivaara, "Deterministic and statistical characterization of rigid frame porous materials from impedance tube measurements," J. Acoust. Soc. Am., vol. 142, no. 4, 2407-2418 (2017).

II M. Niskanen, O. Dazel, J.-P. Groby, A. Duclos, and T. Lähivaara, "Characterising poroelastic materials in the ultrasonic range -A Bayesian approach," J. Sound Vib., vol. 456, 30-48 (2019).

III M. Niskanen, A. Duclos, O. Dazel, J.-P. Groby, J. Kaipio, and T. Lähivaara, "Estimating the material parameters of an inhomogeneous poroelastic plate from ultrasonic measurements in water," J. Acoust. Soc. Am., vol. 146, no. 4, 2596-2607 (2019).

Throughout the overview, these papers will be referred to by Roman numerals.

AUTHOR'S CONTRIBUTION

The publications selected in this dissertation are original research papers on statistical characterisation of porous materials using acoustical methods. All publications are results of joint work with the author of this thesis, the supervisors, and the co-authors. The author of this thesis has been the principal writer and the corresponding author in all of the publications but the help, guidance, and contribution of the thesis supervisors has been significant. The author has also been responsible for conducting the experiments, developing the codes for the forward models and the inversion, and producing the numerical results.

ix Nearly all materials are porous to some degree. In some, the porosity occurs naturally, others are manufactured to be porous, and in fact it is rather difficult to find or manufacture a genuinely non-porous solid [START_REF] Rouquerol | Recommendations for the characterization of porous solids (Technical Report)[END_REF]. In this thesis, we examine materials whose pores form a fluid-filled interconnected network inside a solid frame. Such materials commonly around us include the human body with its bones and tissues, plants, rocks, soils, concrete, plaster, fabrics, wools, and so on [START_REF] Coussy | Mechanics and Physics of Porous Solids[END_REF]. It is well known that the structure of the pores greatly affects the materials' physical properties, such as density, thermal conductivity, and strength. In addition, the connected pore structure allows fluid to flow through the materials, which gives them special properties and enables countless applications. These include insulation and dissipation of thermal energy, absorption of acoustic and mechanical energy, filtration, diffusion, and mixing of flow, which is why porous materials find uses in fields such as civil engineering, geophysical exploration, food industry, and chemical industry to name a few [START_REF] Coussy | Mechanics and Physics of Porous Solids[END_REF]. It is therefore essential that we have reliable methods to characterise the pore structure of materials. Due to the practical importance of porous material characterisation, numerous methods have been proposed for it [START_REF] Rouquerol | Recommendations for the characterization of porous solids (Technical Report)[END_REF][START_REF] Bonfiglio | Inversion problems for determining physical parameters of porous materials: Overview and comparison between different methods[END_REF][START_REF] Liu | Porous Materials[END_REF][START_REF] Horoshenkov | A review of acoustical methods for porous material characterisation[END_REF]. These include measurements of fluid flow, gas adsorption, and X-ray and neutron scattering. One class of methods in particular, based on acoustic probing of the medium, has received a lot of interest in the recent years. Acoustic characterisation methods are based on the fact that pore parameters have a strong influence on the response of porous media to sound waves [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF][START_REF] Cox | Acoustic absorbers and diffusers: theory, design and application[END_REF]. In reverse, the acoustic response carries lots of information on the pore structure, which we can attempt to extract using inversion methods. Acoustical characterisation methods are non-destructive, relatively low-cost, and have the potential to provide information on a multitude of pore structure parameters from a single measurement.

The first step in acoustical characterisation is to solve the forward problem, that is, develop methods to model the acoustic response of porous materials that we can then compare to measurements. Mathematically, we can represent porous media to consist of two phases, a compressible fluid and an elastic solid phase, which both contribute to wave propagation as a coupled fluid-solid system. We call these kind of materials poroelastic. A generally accepted theory to model the deformations in this bi-phasic media is the Biot theory [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range[END_REF][START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range[END_REF][START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF][START_REF] Biot | Mechanics of deformation and acoustic propagation in porous media[END_REF][START_REF] Biot | Generalized Theory of Acoustic Propagation in Porous Dissipative Media[END_REF], which was published by M. A. Biot in 1956 [8, 9]. The Biot theory predicts the existence of three distinct waves in an isotropic poroelastic medium, all of which have been experimentally observed [START_REF] Plona | Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies[END_REF][START_REF] Kelder | Observation of the Biot slow wave in water-saturated Nivelsteiner sandstone[END_REF]. In some cases, such as with many air-saturated porous media, the coupling between the solid and the fluid phases is weak. In such a case, an incident acoustic wave may not transfer enough energy to the solid phase to cause it to vibrate, and we can model the porous medium as having a rigid frame. Only one type of wave, the longitudinal acoustic wave, can propagate in this medium. We can then simplify our models by replacing the poroelastic material with an equivalent fluid model whose properties are defined by the saturating fluid and the solid matrix.

The inverse problem in acoustical characterisation has been traditionally carried 1 out in the deterministic sense. The most straightforward approach is to minimise a cost functional, such as the least squares functional, of the model residuals. Studies that use a least squares minimisation approach to estimate the pore structure parameters include [START_REF] Atalla | Inverse acoustical characterization of open cell porous media using impedance tube measurements[END_REF][START_REF] Sebaa | Measuring flow resistivity of porous material via acoustic reflected waves[END_REF][START_REF] Sebaa | Ultrasonic characterization of human cancellous bone using the Biot theory: Inverse problem[END_REF][START_REF] Fellah | Ultrasonic characterization of porous absorbing materials: Inverse problem[END_REF][START_REF] Ogam | Recovery of elastic parameters of cellular materials by inversion of vibrational data[END_REF][START_REF] Ogam | Non-ambiguous recovery of Biot poroelastic parameters of cellular panels using ultrasonic waves[END_REF][START_REF] Fellah | Simultaneous determination of porosity, tortuosity, viscous and thermal characteristic lengths of rigid porous materials[END_REF][START_REF] Sadouki | Measuring static viscous permeability of porous absorbing materials[END_REF][START_REF] Zieli Ński | Normalized inverse characterization of sound absorbing rigid porous media[END_REF][START_REF] Liu | Inverse characterization of the frequency-dependent acoustic and elastic parameters of porous materials by surface-normal impedance method[END_REF][START_REF] Vanhuyse | Global optimisation methods for poroelastic material characterisation using a clamped sample in a Kundt tube setup[END_REF][START_REF] Verdiere | Inverse Poroelastic Characterization of Open-Cell Porous Materials Using an Impedance Tube[END_REF][START_REF] Lei | Inverse method for elastic properties estimation of a poroelastic material within a multilayered structure[END_REF]. Approaches where the parameters are connected to the data using an analytical inverse mapping have also been proposed [START_REF] Panneton | Acoustical determination of the parameters governing viscous dissipation in porous media[END_REF][START_REF] Olny | Acoustical determination of the parameters governing thermal dissipation in porous media[END_REF][START_REF] Doutres | Evaluation of the acoustic and non-acoustic properties of sound absorbing materials using a threemicrophone impedance tube[END_REF][START_REF] Groby | Analytical method for the ultrasonic characterization of homogeneous rigid porous materials from transmitted and reflected coefficients[END_REF]. When the number of estimated parameters is small, deterministic approaches can provide reliable parameter estimates for porous media, as shown by the review articles [START_REF] Bonfiglio | Inversion problems for determining physical parameters of porous materials: Overview and comparison between different methods[END_REF][START_REF] Horoshenkov | A review of acoustical methods for porous material characterisation[END_REF]. However, with a growing number of estimated parameters, the parameter estimates and especially possible uncertainty estimates often become unreliable.

Recent studies [START_REF] Pompoli | How reproducible is the acoustical characterization of porous media?[END_REF][START_REF] Bonfiglio | How reproducible are methods to measure the dynamic viscoelastic properties of poroelastic media?[END_REF] have shown that the reproducibility of deterministic characterisation methods between different acoustic laboratories is quite poor, with the exception of a few parameters. Some of the poor reproducibility can be explained by inhomogeneities in the measured samples and by a lack of standardised measurement and calibration procedures. On the other hand, acoustic models for porous media include nonlinear interactions between the model parameters, the pore parameters may be strongly correlated with each other, and the cost functions often exhibit several minima. These factors present considerable difficulties for minimising the cost function, and the point estimate provided by the optimisation may not summarise the complex parameter space structures very well. In addition, there are underlying uncertainties associated with the measurements that can be difficult to model, and are therefore often neglected. In the context of inverse problems, however, neglecting the modelling errors may yield meaningless estimates [START_REF] Calvetti | An introduction to Bayesian scientific computing: ten lectures on subjective computing[END_REF][START_REF] Kaipio | Approximate marginalization over modeling errors and uncertainties in inverse problems[END_REF][START_REF] Nicholson | Estimation of the Robin coefficient field in a Poisson problem with uncertain conductivity field[END_REF].

Another approach to the inverse problem is to solve it in the Bayesian framework [START_REF] Calvetti | An introduction to Bayesian scientific computing: ten lectures on subjective computing[END_REF][START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF][START_REF] Kaipio | Statistical and Computational Inverse Problems[END_REF][START_REF] Gelman | Bayesian Data Analysis[END_REF], which has several pragmatic advantages compared to the deterministic solutions. For example, it allows us to be explicit in defining prior information, and provides a way to directly quantify uncertainty. In addition, the Bayesian approach can naturally handle inaccuracies between the models and actual experiments. In the Bayesian framework, all unknowns are modelled as random variables, and the inverse problem takes the form of statistical inference. The solution of the inverse problem is then the posterior probability density (ppd), which combines our knowledge on the measurement and model uncertainties, possible prior information, and the measurement data. Once the posterior model has been derived, we can compute point and interval estimates to summarise the information in an easily understandable form. The maximum a posteriori (MAP) and conditional mean (CM) estimates, defined in Chapter 3, are examples of Bayesian point estimates. The MAP estimate can be solved with an optimisation approach, but finding the CM estimate, as well as computing the interval estimates, is an integration problem. In the recent years, Bayesian methods have started to gain popularity for characterising porous materials, especially ones that can be modelled as an equivalent fluid [START_REF] Dettmer | Full wave-field reflection coefficient inversion[END_REF][START_REF] Dettmer | Trans-dimensional geoacoustic inversion[END_REF][START_REF] Holland | In situ sediment dispersion estimates in the presence of discrete layers and gradients[END_REF][START_REF] Fackler | Bayesian acoustic analysis of multilayer porous media[END_REF][START_REF] Bonomo | A comparison of three geoacoustic models using Bayesian inversion and selection techniques applied to wave speed and attenuation measurements[END_REF], but methods for poroelastic inversion have also been proposed [START_REF] Bonomo | A comparison of three geoacoustic models using Bayesian inversion and selection techniques applied to wave speed and attenuation measurements[END_REF][START_REF] Chazot | Acoustical and mechanical characterization of poroelastic materials using a Bayesian approach[END_REF].

In this thesis, the aim is to develop robust acoustical characterisation methods for both rigid and elastic frame porous materials, using measurements done with typical existing instruments. In particular, we are interested in quantifying the uncertainty in the parameter estimates so that the repeatability and reproducibility of the characterisation methods could be improved. To this end, we adopt the Bayesian framework for inverse problems. To explore the true posterior density, we then use Markov chain Monte Carlo (MCMC) [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] methods, which are often the only practical way to draw samples from a high-dimensional and possibly multi-modal posterior. However, MCMC methods can be computationally very demanding, because often to obtain a representative sample of the posterior the forward problem needs to be solved hundreds of thousands of times, possibly even more. This is why we pay close attention to the speed of the forward models used, and use fast matrix methods to solve the Biot equations. In addition, we have been able to take advantage of the numerous advances to MCMC in the recent years [START_REF] Haario | An adaptive Metropolis algorithm[END_REF][START_REF] Earl | Parallel tempering: Theory, applications, and new perspectives[END_REF][START_REF] Andrieu | A tutorial on adaptive MCMC[END_REF][START_REF] Vousden | Dynamic temperature selection for parallel tempering in Markov chain Monte Carlo simulations[END_REF], which have made it possible to construct efficient samplers for poroelastic wave propagation problems.

The results of this thesis have been published in I-III.

In publication I, we examine rigid frame porous media that can be measured in a standard transmission impedance tube [START_REF] Song | A transfer-matrix approach for estimating the characteristic impedance and wave numbers of limp and rigid porous materials[END_REF]. These materials include acoustic foams, automotive felt materials, and different kinds of wools. We solve the inverse problem in the Bayesian framework, but in addition compare the results to inversion from a least-squares minimisation approach. We make the comparison to highlight the differences between the traditional way of solving the characterisation inverse problem, and the Bayesian approach.

In publications II and III, we focus on porous materials with an elastic frame. First, publication II provides a numerical study into the feasibility of estimating basically all Biot input parameters, just from acoustic reflection and transmission data in the ultrasonic frequencies. We consider some stability issues of the forward model, and then construct an efficient MCMC sampler to explore the posterior. We use methods such as Adaptive Metropolis (AM) [START_REF] Haario | An adaptive Metropolis algorithm[END_REF] and Parallel Tempering (PT) for the MCMC sampler and find that it can explore the peaky posterior density reasonably efficiently.

Publication III is a study into the characterisation of a poroelastic ceramic material, submerged in water so that we can excite the solid waves. We examine numerous measurements from several parts of a object made from the porous ceramic, and carry out the inversion separately in each measurement spot. This highlights the heterogeneity of the sample, and provides a way of quantifying it.

The introduction to this thesis is organised as follows. Chapter 2 presents briefly the basics of acoustic wave propagation in rigid and elastic frame porous media, and matrix methods for their numerical simulation. Chapter 3 introduces the Bayesian framework for inverse problems, along with the sampling methods used in this thesis. Chapter 4 reviews and studies the results of the publications I-III, and a conclusion is given in Chapter 5.

Wave propagation in porous media

In this chapter, we will review the basics of sound propagation in rigid and elastic frame porous media. Rigid frame porous media were considered in publication I, and elastic frame materials in publications II and III. We first consider mathematical models that describe the acoustically relevant properties of porous media, and define the parameters that are related to the pore properties. We also state the acoustic and poroelastic wave equations. At the end of the chapter we discuss the construction of the forward models and their computational implementation, and especially some fast matrix methods for solving the wave equations.

Porous materials consist of a continuous solid skeleton that forms the frame of the material and of a fluid phase that fills the pores in the solid, such as in the porous rock of Figure 2.1. For most applications, we are interested in materials whose pores form an interconnected network, which allows acoustic propagation within the porous frame. We call these kind of materials open-cell, as opposed to closed-cell porous materials where each pore is isolated from their neighbour. In addition, we can define a saturated porous medium, where all the pores are filled with the same fluid. This simplifies their modelling considerably. An example of an unsaturated porous medium would be a porous rock submerged in water where air is trapped inside some of the pores. In this thesis, we will consider saturated open-cell porous media, which we can model as a superimposition of two continua (the frame and the fluid) [START_REF] Coussy | Poromechanics[END_REF].

Figure 2.1:

An example of a naturally occurring porous material, scoria, which is volcanic rock solidified from foaming lava. Scoria typically has a high porosity and pore connectivity [START_REF] Saar | Permeability-porosity relationship in vesicular basalts[END_REF]. The rock is from the Maungauika volcano in Auckland, New Zealand.

Modelling attenuation mechanisms is of prime importance when modelling wave propagation in porous media. Acoustic waves in porous media are mainly attenuated through viscous and thermal losses that result from the interaction between the solid frame and the fluid. In the free field and at short distances, attenuation can be usually neglected without introducing large errors. However, when waves are propagating in small ducts (which porous media closely resemble), the effects of viscosity and heat conduction become important. These effects were understood already in 1868, when Kirchhoff published a general theory for sound propagation in tubes with an arbitrary cross section [START_REF] Kirchhoff | Ueber den Einfluss der Wärmeleitung in einem Gase auf die Schallbewegung[END_REF]. In his work, Kirchhoff considered the viscous and thermal effects together in the same equations [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF]. In the case of cylindrical tubes, Zwikker and Kosten [START_REF] Zwikker | Sound absorbing materials[END_REF] worked out simpler expressions, where the viscous and thermal losses are treated separately. Stinson [START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF] later generalised these equations for an arbitrary cross section, and the use of separate expressions for viscous and thermal losses for sound propagation in ducts is now commonly used.

Analytical description of sound propagation in porous media is, in general, impossible due to the complex geometry of their microstructure. Therefore, we have to model porous media in the macroscopic sense, relying on homogenised representations of the material. Such simplifications are possible when the long-wavelength condition is met, that is, acoustic wavelengths are much larger than the homogenised representative volume of the material, which in turn is much larger than the average pore size [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF].

RIGID FRAME POROUS MEDIA

Often in cases where a porous material is saturated with air and excited by an airborne wave, we can assume the solid frame is motionless. In general, the frame can be assumed motionless when the solid-fluid coupling is weak, and little energy is transferred to the solid. The strength of this coupling is a function of frequency, and it gets weaker as the frequency increases [START_REF] Zwikker | Sound absorbing materials[END_REF]. A porous medium with a motionless solid phase can be modelled as a fluid whose apparent properties are a function of the solid properties [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF]. We call such a medium an equivalent fluid.

Time-harmonic wave propagation in an isotropic equivalent fluid medium can be described by the Helmholtz-like equation

∇ 2 + k 2 eq p = 0, (2.1)
where ∇ 2 is the Laplace operator, p denotes the acoustic pressure, and keq is the equivalent wave number in the fluid. We adopt the exp(-iωt) time convention, where i = √ -1, ω = 2π f is the angular frequency, and f is the frequency. The wave number can be defined with the equivalent density ρeq and the equivalent adiabatic bulk modulus Keq as keq = ω ρeq Keq .

(2.

2)

The tilde over the quantities denotes that they are complex-valued (due to attenuation) and frequency-dependent (because the governing physics is different at low and high frequencies). The acoustic behaviour of a fluid medium is completely described by ρeq and Keq , and from these parameters we can derive interesting properties such as the fluid's equivalent speed of sound, ceq , and the equivalent characteristic impedance Zeq as ceq = Keq ρeq , and Zeq = ρeq Keq .

(2.3)

Any two of the above parameters can be used to describe the medium's acoustic behaviour.

Several models exist to express the acoustic properties ρeq and Keq of an equivalent fluid. The simplest of these are empirical. They depend on very few macroscopic parameters, and are derived based on fitting suitable simple equations to a set of measurements. One of the first and the most well-known empirical models is the Delany-Bazley (D-B) model [START_REF] Delany | Acoustical properties of fibrous absorbent materials[END_REF], where a material's characteristic impedance and wave number are stated as a function of only one parameter, the flow resistivity. Some limitations of the D-B model are that it is only suitable for fibrous porous materials (such as various glass fibres), and that it is valid only for a restricted frequency range. In the low frequencies, the model returns non-causal solutions [START_REF] Miki | Acoustical properties of porous materials -Modifications of Delany-Bazley models[END_REF] that violate the Kramers-Krönig relations [START_REF] Waters | Causality-imposed (Kramers-Kronig) relationships between attenuation and dispersion[END_REF]. Since their introduction, several improvements to the D-B model have been made, which address the causality requirement but also widen the frequency range and types of applicable materials [START_REF] Miki | Acoustical properties of porous materials -Modifications of Delany-Bazley models[END_REF][START_REF] Qunli | Empirical relations between acoustical properties and flow resistivity of porous plastic open-cell foam[END_REF][START_REF] Komatsu | Improvement of the Delany-Bazley and Miki models for fibrous sound-absorbing materials[END_REF][START_REF] Oliva | Sound absorption of porous materials -Accuracy of prediction methods[END_REF][START_REF] Kirby | On the modification of Delany and Bazley fomulae[END_REF][START_REF] Jones | Simplification of the Delany-Bazley approach for modelling the acoustic properties of a poroelastic foam[END_REF].

Empirical models are simple, but clearly not well suited for use in material characterisation due to the fact that they contain very few physical material parameters, and do not model the materials with enough accuracy to do inversion. For this reason, we will now consider semi-phenomenological models that are at least partly based on the physics of sound propagation in porous media. Again, there are many models available, such as the Wilson model [START_REF] Wilson | Relaxation-matched modeling of propagation through porous media, including fractal pore structure[END_REF] where the thermal and viscous effects are modelled as relaxation processes, or the pore size distribution model [START_REF] Horoshenkov | Asymptotic limits of some models for sound propagation in porous media and the assignment of the pore characteristic lengths[END_REF][START_REF] Horoshenkov | A three-parameter analytical model for the acoustical properties of porous media[END_REF] which assumes that the pore sizes are log-normally distributed. In this thesis, we will use the Johnson-Champoux-Allard-Lafarge (JCAL) model, which is currently perhaps the most commonly used model and has been found to represent porous media accurately in many studies.

The Johnson-Champoux-Allard-Lafarge model

The JCAL model [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF][START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF][START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF] is a popular semi-phenomenological model that describes the equivalent density and bulk modulus of rigid frame porous media based on six physical parameters, defined in the following, that are independent of the saturating fluid. The JCAL model is directly used as the forward model in publication I, since our measurement procedure provides the materials' densities and bulk moduli. The model is based on considering the asymptotic behaviour of the density and bulk modulus at low and high frequencies, and joining the asymptotes together with simple expressions. In the JCAL model, viscous losses are treated in the expression of the density and thermal losses are treated in the expression of bulk modulus. Let us first define the physical parameters, while keeping in mind that they are defined generally for any porous media, and are not dependent on the JCAL, or any other model, that uses them.

Open porosity φ.

Porosity is one of the most important parameters describing a porous medium. It is defined as the volume of the fluid V f divided by the total volume of the material Static viscous permeability k 0 . Alongside porosity, viscous permeability k(ω) is the other main parameter describing porous media, and the most important parameter for describing sound absorbing materials. Viscous permeability describes the ease at which fluid can flow through a porous material, and is defined by Darcy's law:

V t : φ = V f V t . ( 2
k(ω) = - ηφ v ∇p , (2.5)
where ∇p is the pressure gradient, v is the average fluid velocity in a porous medium, and η is the dynamic viscosity of the saturating fluid. The static limit k 0 = lim ω→0 k(ω) depends only on the structure of the porous frame. Flow resistivity σ is the reciprocal of permeability σ = η/k 0 . The values of permeability can span several orders of magnitude, from practically impermeable 10 -20 m 2 to highly permeable 10 -8 m 2 [START_REF] Schön | Physical properties of rocks: Fundamentals and principles of petrophysics[END_REF].

Geometrical tortuosity α ∞ . When the pore network is not perfectly straight and parallel to the direction of wave propagation, the waves travel along longer paths. This increase in the path length results in an apparent increase of the density of the material compared to the density of the saturating fluid ρ f . At high frequencies the viscous skin depth δ = (2η/ρ f ω) 1 2 becomes tiny compared to the pore size, and the fluid behaves approximately as an ideal, non-viscous, fluid [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF]. In this limit, we can define the geometrical tortuosity α ∞ , which is a real-valued parameter that relates the effective density to the density of the fluid as ρ e = α ∞ ρ f . Theoretically, tortuosity is also connected to k(ω) [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF], and in fact k 0 and α ∞ can be seen as parameters describing the effective density of the material at low and high frequencies, respectively.

Geometrical tortuosity can be defined with the help of fluid velocities. Let us denote the microscopic fluid velocity of nonviscous fluid flow at r by v m (r). Tortuosity can then be defined as

α ∞ = V v 2 m (r)dV V v m (r)dV 2 , (2.6)
where the integrals are carried out over the homogenized representative volume V.

We can also define tortuosity with a simpler concept, related to the path lengths.

Let L be the length of a straight path through a porous material, and L the actual tortuous path (see Figure 2.2). Then have the greatest influence on Λ , and the constricted regions of the pore space have the greatest influence on Λ.

α ∞ = L L 2 .
(2.7)

By definition, α ∞ ≥ 1, and typical values of tortuosity range between 1 and 3.

Characteristic viscous length Λ. In many applications, we are interested in the ratio of pore volume to surface area. In the JCAL model, the characteristic viscous and thermal lengths are pore shape factors and provide an indication of this ratio.

The characteristic viscous length Λ describes viscous effects at middle and high frequencies, and is defined as [69]

Λ = 2 V |v i (r)| 2 dV A |v i (r w )| 2 dA , ( 2.8)
where the integrals are carried over the pore space, v i (r) is the fluid velocity in the pore, and v i (r w ) is the fluid velocity on the pore wall. The definition of Λ is such that constrictions in the pore network (see Figure 2.2), where the velocities are greater, contribute to the integral the most. It can therefore be thought of as an indication of the size of the constrictions. Typically, the viscous characteristic length ranges from 10 µm to 500 µm.

Characteristic thermal length Λ . Thermal effects at high frequencies are described by the characteristic thermal length Λ . It is defined in a similar way to Λ, except that it does not consider fluid velocity in the pores, but only the ratio of average pore volume to average pore surface area [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF]:

Λ = 2 V dV A dA .
(2.9)

As the exchange area between the solid and the fluid is greater in the open regions (see Figure 2.2), these areas contribute to Eq. (2.9) the most. Therefore, thermal characteristic length can be seen as an indication on the size of the open pore areas. By definition, Λ ≥ Λ, and typically for fibrous sound absorbing materials Λ = 2Λ [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF].

Static thermal permeability k 0 . Analogously to the dynamic viscous permeability, the dynamic thermal permeability k (ω) relates the pressure gradient to the macroscopic excess temperature τ as [START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF] k

(ω) = - κφ τ iωp , (2.10)
where κ denotes thermal conductivity. The static thermal permeability

k 0 = lim ω→0 k (ω),
describes the thermal effects at low frequencies.

Expressions for ρeq and Keq

Now we can state the equations for the density and bulk modulus in the JCAL model. The equivalent density ρeq of the material is defined

ρeq = ρ f φ α(ω) , (2.11)
and the dynamic tortuosity is written as

α(ω) = α ∞ + iν ω φ k 0 1 - iω ν 2α ∞ k 0 φΛ 2 , (2.12)
where ν = η/ρ f is the kinematic viscosity.

The expression for the equivalent bulk modulus Keq is given by

Keq = γP 0 φ γ - γ -1 α (ω) -1 , (2.13)
where γ is the specific heat ratio, P 0 is the static fluid pressure, and α (ω) is called the thermal tortuosity, as a homologue to α(ω). The quantity α (ω) is written as

α (ω) = 1 + iν ω φ k 0 1 - iω ν 2k 0 φΛ 2 , (2.14)
where ν = ν/Pr, and Pr is the Prandtl number. The Prandtl number is an approximately constant property of fluids, describing the ratio of momentum diffusivity to thermal diffusivity.

POROELASTIC MEDIA

In general, the frame of a porous medium is not motionless and vibrates in response to an acoustic excitation. This happens especially when the fluid-solid coupling is strong, for example when the saturating fluid is dense, such as water or oil, when the material the frame is made of is easily deformed, or when the acoustic frequencies are low. Obviously, the solid will vibrate from direct mechanical excitation as well. Accurate prediction of sound propagation in these cases can only be achieved if we consider a model where the fluid and the frame move simultaneously. The Biot theory [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range[END_REF][START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range[END_REF] describes this situation, and is commonly considered to be the most accurate model for poroelastic media.

The Biot model

Biot's poroelastic wave equations can be thought of as coupling the linear elastodynamic equation for wave propagation in solids, and the Helmholtz equation for wave propagation in fluids. Biot derived his equations from the principles of continuum mechanics, assuming the existence of potentials and the principle of least action [START_REF] Bourbié | Acoustics of porous media[END_REF]. Later it has been shown that the same equations can be derived by a homogenisation approach starting from the microstructure [START_REF] Burridge | Poroelasticity equations derived from microstructure[END_REF].

The original Biot model [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range[END_REF][START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range[END_REF], published in 1956, represents a porous medium with six parameters, which are the displacements of the solid phase u and the fluid phase U. Some assumptions of the theory are isotropy and homogeneity of the medium, and that the long-wavelength condition is fulfilled. In 1962, Biot published an alternative formulation of the theory [START_REF] Biot | Mechanics of deformation and acoustic propagation in porous media[END_REF], which describes the medium with the solid displacement and a fluid/solid relative displacement w = φ(Uu). This formulation is more flexible because it allows the treatment of inhomogeneous materials and simplifies some boundary conditions, and we therefore use it in this thesis. Other representations of the Biot equations have also been proposed [START_REF] Atalla | A mixed displacement-pressure formulation for poroelastic materials[END_REF][START_REF] Dazel | An alternative Biot's displacement formulation for porous materials[END_REF].

We will next describe the basic equations and parameters in the Biot theory. Further information, and detailed derivation of the theory, can be found for example in [START_REF] Coussy | Mechanics and Physics of Porous Solids[END_REF][START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF][START_REF] Bourbié | Acoustics of porous media[END_REF]. In Biot's alternative formulation [START_REF] Biot | Mechanics of deformation and acoustic propagation in porous media[END_REF], we can write the equations of motion in the frequency domain as

ω 2 ρ f w + ω 2 ρu = -∇ • σ, ω 2 ρ f u + ω 2 ρeq w = ∇ • p, (2.15)
where σ is the total stress tensor, ρ = (1φ)ρ s + φρ f is the bulk density of the medium, ρ s is the density of the solid, and the rest of the parameters are defined in the same way as in the rigid frame case. In particular, we use the Johnson model of dynamical tortuosity (2.11)-(2.12) in the expression of equivalent density ρeq . The stress-strain relations in the alternative formulation are of the form

σ = 2N + (λ c ∇ • u + α B M∇ • w)I, p = -M(∇ • w + α B ∇ • u).
(2.16)

Here = 1 2 ∇u + (∇u) T is the strain tensor, and I denotes an identity tensor. The quantities N, α B , λ c , and M are elastic coefficients. The parameter N is the classic shear modulus and α B is the Biot-Willis coefficient [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF]. The Biot-Willis coefficient controls the rigidity of the frame and is bounded by φ ≤ α B ≤ 1, where the lower limit corresponds to a rigid frame and the upper limit to a completely limp frame. In addition, λ c = λ + α 2 B M, where λ and M are other elastic parameters. It can be shown that when the parameters in the model are that of a perfect fluid (φ = 1, α ∞ = 1, N = 0, k 0 → ∞) the Biot equations reduce to the linear lossless acoustic wave equation, and when the parameters are that of a pure solid (φ = α B = 0, α ∞ → ∞), we get the elastic wave equation [START_REF] Bourbié | Acoustics of porous media[END_REF].

Rather than use the elastic parameters defined by Biot, we would like to characterise materials with more commonly used parameters. We therefore relate the α B , M, and λ c to porosity, shear modulus, and the bulk moduli K f , K s (the bulk modulus of the elastic solid from which the frame is made of), and K b (the bulk modulus of the porous frame in a vacuum), with the following relations

α B = 1 - K b K s , (2.17)
M = α B + K s K f -1 φ -1 K s , (2.18
)

λ c = (1 + φ)α B + K b K f -1 φ M - 2 3 N, (2.19)
see [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF][START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF].

Waves in a poroelastic medium

A fluid medium can only support the compressional wave whereas an isotropic elastic medium can support both compressional and shear wave propagation. We can therefore expect that these three waves will also propagate in a poroelastic medium, a consolidation of a solid and a fluid phase. Indeed, in a homogeneous isotropic media the Biot equations admit three different kinds of solutions which are the three types of waves that can propagate in porous media, called the fast and slow longitudinal waves and the shear wave. However, the waves do not propagate in a single phase in isolation, but all waves propagate both in the solid and the fluid phase and can be distinguished by their velocity and mode of propagation. The faster of the two compressional waves corresponds to vibration in which overall and fluid displacements are nearly in phase, and the slower one to vibration where they are out of phase. For the shear wave the movements are in phase since the fluid does not exert shear force and influences the propagation only through inertial effects [START_REF] Bourbié | Acoustics of porous media[END_REF]. Figure 2.3 depicts the generation of the three Biot waves by mode conversion at the poroelastic boundary. A compressional wave is initially propagating in a fluid medium, meets a poroelastic slab at an angle and transmits a part of its energy into the poroelastic material (reflections are not shown) where the propagation continues as Biot waves. The waves refract at different angles depending on their velocity, and can eventually be recorded after they are converted back to a compressional wave in the fluid. Note that at normal incidence, a plane wave excitation does not generate a shear wave because there is no sideways excitation at the interface.

Losses in water-saturated materials

Let us now consider attenuation of the waves in a poroelastic material. Like in the rigid frame situation, viscous effects in the fluid are an important loss factor, and we can use the dynamical tortuosity [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF] model (2.11)-(2.12) to account for them. Possible thermal losses in the saturating fluid can also be accounted for with the Champoux-Allard-Lafarge model [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF][START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF] (2.13)- (2.14). However, thermal effects can often be neglected when the medium is saturated with a liquid due to liquids' compressibility being much lower than that of gases.

Figure 2.3:

A schematic showing the three Biot waves in a porous medium (after [START_REF] Plona | Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies[END_REF]). Reflected waves are not shown.

Another source of attenuation in porous media are the viscoelastic losses in the solid frame. A viscoelastic material responds to an applied stress with a small delay, which means that the behaviour is time-dependent and the material can be said to have a memory [START_REF] Carcione | Wave fields in real media: Wave propagation in anisotropic, anelastic[END_REF]. In the frequency domain, we can model linear viscoelastic wave propagation with a complex elastic modulus, where the real part determines the phase velocity and the imaginary part determines the damping. When damping in the solid is small, we can use real and imaginary parts that are constant with frequency as the simplest attenuation model, but this is in general not valid [START_REF] Turgut | An investigation of causality for Biot models by using Kramers-Krönig relations[END_REF]. The requirement of causality, which in the frequency domain is expressed with the Kramers-Krönig relations [START_REF] Waters | Causality-imposed (Kramers-Kronig) relationships between attenuation and dispersion[END_REF][START_REF] Bourbié | Acoustics of porous media[END_REF], necessitates that the real and imaginary parts are frequency dependent.

Let us consider a general complex elastic modulus M(ω) = M R (ω) + iM I (ω), which consists of frequency-dependent real and imaginary parts. We can now define a dimensionless quality factor Q, which describes the dissipation in the medium:

Q(ω) = M R (ω) M I (ω) . ( 2

.20)

A highly dissipative medium has a small Q, whereas for a purely elastic medium Q → ∞ and the elastic modulus is real and constant. Note that the quality factor is always defined as (2.20), but different viscoelasticity models give different values of Q.

Models for describing the behaviour of Q with respect to frequency are often derived based on simple mechanical models, springs and dampers. The simplest ones are the Maxwell and Kelvin-Voigt models, in which a spring and a damper are connected in series and in parallel, respectively [START_REF] Carcione | Wave fields in real media: Wave propagation in anisotropic, anelastic[END_REF]. These models are not very realistic, though, because the Maxwell model predicts infinite damping at zero frequency and no damping when ω → ∞, and the Kelvin-Voigt model the opposite.

A more realistic model is the Zener or standard linear solid model, which combines a spring and a Kelvin-Voigt model. Attenuation in the Zener model is again frequency dependent, but now has a peak at a central frequency ω m , and approaches zero as ω → 0 and ω → ∞. It is also possible to construct an attenuation profile that has multiple peaks by connecting several Zener elements in parallel. This model is called the generalised Zener model.

The Zener model and especially its generalised version give a lot of freedom to represent an attenuation profile of a material. However, the downside is that each Zener element adds three free parameters that we need to set (or estimate during inversion). Often the frequency dependence of the quality factor is not known, but some experiments have shown that Q stays nearly constant over a wide range of frequencies [START_REF] Bourbié | Acoustics of porous media[END_REF]. We can then make use of the constant Q model, derived by Kjartansson in 1979 [START_REF] Kjartansson | ConstantQ-wave propagation and attenuation[END_REF], which is a causal model where the quality factor is constant over the whole frequency range. The advantages of using this model are that it only includes two parameters, and is mathematically very simple. A complex modulus in the constant Q model is expressed as

M(ω) = M 0 iω ω 0 2 π arctan(1/Q) , (2.21)
where ω 0 is a reference frequency, M 0 is the modulus of M(ω) at ω 0 , and Q is the quality factor.

MODELS TO COMPUTE THE ACOUSTIC RESPONSE OF A MUL-TILAYER SYSTEM

Let us now consider the construction of a forward model for the characterisation problem in publications II and III. This means constructing a model that predicts an acoustic response of a layered system that we can measure (in our case the reflection and transmission coefficients), when the equations governing wave propagation and the physical parameters of the system are known. In this thesis, motivated by computational efficiency, we use plane waves to model the wave propagation. Sound propagation in complex-shaped structures has to be usually solved with element-based methods such as the Finite Element [START_REF] Ihlenburg | Finite element analysis of acoustic scattering[END_REF], Ultra Weak Variational Formulation [START_REF] Cessenat | Application of an Ultra Weak Variational Formulation of Elliptic PDEs to the Two-Dimensional Helmholtz Problem[END_REF], or Discontinuous Galerkin [START_REF] Wilcox | A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media[END_REF] methods. These methods require a certain number of elements per acoustic wavelength and with high frequencies can become intolerably complex computationally, especially with inverse problems that we wish to solve by sampling methods where we have to repeatedly solve the forward problem. In this case the computational efficiency of the forward model can determine whether or not the problem is solvable at all in practice. A considerable saving of computational effort can be made if we have a layered system. This allows us to use a number of matrix methods to analytically represent the wave propagation [START_REF] Lowe | Matrix techniques for modeling ultrasonic waves in multilayered media[END_REF], and solve the motion equations in a fraction of the time it takes with element-based methods.

Waves in layered media

Let us consider a layered system as in Figure 2.4, where a plane wave initially propagating in medium Ω [0] is incident on the system consisting of n layers. The layers are assumed to be flat, stacked along the x 2 -dimension, and to extend to infinity with Ω [i] and the locations of the interfaces by I i .

in the other directions. The layers can consist of a fluid (or an equivalent fluid), an elastic solid, or a saturated poroelastic material. The termination condition can either be a rigid backing, in which case there is no transmission, or another layer extending to infinity in the x 2 -direction. Assuming isotropy, the coordinate axis can always be rotated so that there is no propagation in the x 3 direction, and the problem reduces to two dimensions [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF]. Let us represent an arbitrary field in the frequency domain as χ(x 1 , x 2 ; ω). Because of the plane wave assumption, we can perform a Fourier transform along x 1 and write the field as χ(x; ω) = χ(k 1 , x 2 ; ω)e ik 1 x 1 , where k 1 = -k [0] sin(ϕ inc ), and k [0] = ω/c f . We will omit the ω and k 1 dependence for clarity, and write χ(k 1 , x 2 ; ω) ≡ χ(x 2 ).

In each layer, equations for the displacements, stresses, and pressure can be stated as the superposition of the wave fields in the layer. These physical fields are solutions to the wave equation (acoustic, elastic, or poroelastic depending on the material) in an infinite medium, and can be expressed as a function of the complex amplitude A and wave number k type 2 of each wave, where type indicates the mode of the wave. The idea of multilayer matrix methods is then to combine these solutions of an infinite medium together by introducing boundary conditions at interfaces between the layers. In each layer we can have a backward and a forward propagating wave, so in a fluid layer we can have two, in an elastic layer four, and in a poroelastic layer six wave fields. We can thus write an expression of, for example, the x 2 solid displacement û in a poroelastic material as

û2 (x 2 ) =A 1 e -ik fast 2 x 2 + A 2 e ik fast 2 x 2 + A 3 e -ik shear 2 x 2 + A 4 e ik shear 2 x 2 + A 5 e -ik slow 2 x 2 + A 6 e ik slow 2 x 2 , (2.22)
where amplitudes 1 and 2 correspond to the fast longitudinal wave, amplitudes 3 and 4 correspond to the shear wave, and amplitudes 5 and 6 correspond to the slow longitudinal wave. We can state the other displacements and stresses as a function of the wave amplitudes a = [A 1 , . . . , A 6 ] T as well, using the stress-strain relations (2.16). This derivation can be quite cumbersome, however, and in practice we can use the following simpler method based on state vector formalism.

Let us combine the displacements, stresses, and pressure into a state vector ŝ(x 2 ). The state vector components can be selected arbitrarily as long as they are continuous at the interfaces, and it can be shown that for a given layer the number of the state vector components is equal to the number of the forward and backward propagating waves in the layer [START_REF] Dazel | A stable method to model the acoustic response of multilayered structures[END_REF]. For example for fluids it is convenient to select pressure and velocity normal to the interface since they are continuous over boundaries. For a poroelastic material we can choose ŝ(x 2 ) = [ σ12 , σ22 , p, û1 , û2 , v2 ], i.e. the normal components of the total stress tensor, σ12 and σ22 , the pressure p, the solid displacements û1 and û2 , and the normal component of the total velocity v2 = -iω( ŵ2 + û2 ).

Using the state vector, we can write the equations of motion for any type of layer (see for example [START_REF] Gautier | Propagation of acoustic waves in a one-dimensional macroscopically inhomogeneous poroelastic material[END_REF]) as a system of ordinary differential equations as dŝ(x 2 )

dx 2 + Aŝ = 0, (2.23)
where A is the state matrix, which depends on the components selected for the state vector, and is a function of the material parameters, properties of the incident wave, and frequency. We can use this approach to model inhomogeneous materials too, in which case A will also be a function of x 2 . It can be shown (see [START_REF] Dazel | A stable method to model the acoustic response of multilayered structures[END_REF]) that we can now write a matrix equation where the state vector is given as a function of the wave amplitudes a as

ŝ(x 2 ) = [M(x 2 )]a. (2.24)
To find the field matrix [M], we calculate the eigendecomposition of A = ΦΓΦ -1 , where Φ and Γ are matrices that include the eigenvectors and the eigenvalues, re-

spectively. Then [M(x 2 )] = Φ exp{-(x 2 -x ref )Γ}, where x ref is a reference location.

The transfer matrix method

The transfer matrix method (TMM), also known as the Thomson-Haskell method [START_REF] Thomson | Transmission of elastic waves through a stratified solid medium[END_REF][START_REF] Haskell | The dispersion of surface waves on multilayered media[END_REF], is a well-known general purpose multilayer modelling tool which can be used to calculate the acoustic response of layered systems. The method assumes that the propagating waves are plane waves, but finite fields can be studied by superposition of the plane wave solutions [START_REF] Pialucha | The reflection coefficient from interface layers in NDT of adhesive joints[END_REF]. The basic idea in the TMM is to represent the propagation in each layer with a matrix that relates the state vector of the system at the back of the layer to that at the front of the layer [START_REF] Lowe | Matrix techniques for modeling ultrasonic waves in multilayered media[END_REF]. In the following, the notations are as in Figure 2.4. The goal in the TMM is to find a matrix T which relates the state vector ŝ [1] (I 1 ) of layer Ω [1] at the first interface to the state vector ŝ [1] (I 1 ).

[n-1] (I n ) of layer Ω [n-1] at the last interface, i.e. ŝ[n-1] (I n ) = T ŝ
(2.25)

The state vectors at the beginning and end of the layers can then be connected to the known incident, unknown reflected, and unknown transmitted (if no rigid termination) waves, which lets us calculate the acoustic response of the system by solving a matrix equation.

We can find the matrix T in the following way. The amplitudes of the waves inside the same layer are obviously the same near both interfaces, so based on (2.24) we can write for the layer Ω [1] :

[M [1] (I 1 )] -1 ŝ [1] (I 1 ) = a [1] = [M [1] (I 2 )] -1 ŝ [1] (I 2 ).

(2.26)

This means we can find a single layer transfer matrix:

ŝ [1] (I 2 ) = [T [1] ]ŝ [1] (I 1 ),

where [T [1] ] = [M [1] (I 2 )][M [1] (I 1 )] -1 . Now, if layer Ω [2] is made of the same type of material as layer Ω [1] , and the components of the state vectors are continuous at the boundary, we must have ŝ [2] (I 2 ) = ŝ [1] (I 2 ) and therefore ŝ [2] (I 2 ) = [T [1] ]ŝ [1] (I 1 ).

(2.28)

If the layers are of different materials we have to include interface matrices, but for simplicity we do not consider them here and direct the reader to [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF][START_REF] Brouard | A general method of modelling sound propagation in layered media[END_REF]. The same process can be repeated until all layers have been considered: [1] ]ŝ [1] (I 1 ).

ŝ[n-1] (I n ) = [T [n-1] ][T [n-2] ] • • • [T
(2.29)

One of the strengths of the TMM is that we can describe multilayered systems by a single transfer matrix, simply by multiplying together the matrices of the single layers: [1] ].

T = [T [n-1] ][T [n-2] ] • • • [T
(2.30)

Finally, to find the response of the layered system, we express it in terms of the waves in the first and the last medium:

a [n] = [M [n] (I n )] -1 T [M [0] (I 1 )]a [0] . (2.31)
In a usual case, the first and last layers are fluids (this corresponds to a layered system in air or water for example), and we can then have a longitudinal wave propagating towards and away from the system at both ends. To find the reflection and transmission coefficients R and T, respectively, we set the amplitude of the incoming wave to 1 at layer Ω [0] , and to 0 at layer Ω [n] so there is not a wave propagating towards the system from the back. Then the amplitude of the reflected wave corresponds to R and the transmitted wave to T, i.e. a [0] = [1, R] T and a [n] = [0, T] T . Now we have two equations and two unknowns so the coefficients can be solved.

An important issue with the TMM is that it is known to exhibit instability at high frequencies and/or large layer thicknesses. This is known as the "large f d problem" [START_REF] Lowe | Matrix techniques for modeling ultrasonic waves in multilayered media[END_REF]. The problem is of numerical origin (finite precision computers) and is caused mainly by the poor condition number of the transfer matrix, due to the combination of both exponentially decaying and growing components. We observed this instability in our simulations, which meant the TMM approach could not be used. Even though the cause is well-known, the proposed solutions to this differ, and are also dependent on the application. One option is the information vector method [START_REF] Dazel | A stable method to model the acoustic response of multilayered structures[END_REF]. Another solution is to use the so-called global matrix method [START_REF] Knopoff | A matrix method for elastic wave problems[END_REF], which was adopted in this thesis.

The global matrix method

The global matrix method (GMM) [START_REF] Lowe | Matrix techniques for modeling ultrasonic waves in multilayered media[END_REF][START_REF] Knopoff | A matrix method for elastic wave problems[END_REF][START_REF] Chin | Matrix methods in synthetic seismograms[END_REF][START_REF] Schmidt | Efficient global matrix approach to the computation of synthetic seismograms[END_REF] is an alternative matrix formulation to the TMM, originally proposed in [START_REF] Knopoff | A matrix method for elastic wave problems[END_REF]. The main idea in the GMM is to write the forward and backward propagationg waves inside the material explicitly as a function of their amplitudes, and choose the origin of the waves to be on the interface they originate from. We then assemble the equations for these waves in a single global matrix that represents the whole system. A disadvantage of the GMM is that the resulting matrix can be quite big, which can lead to longer calculation times. A significant advantage, however, is that this method is stable for all frequencies and thicknesses [START_REF] Chin | Matrix methods in synthetic seismograms[END_REF][START_REF] Schmidt | Efficient global matrix approach to the computation of synthetic seismograms[END_REF]. The stability comes from the freedom to choose the origin of each wave on the interface where they originate, and the exponentially growing and decaying components are not coupled from one interface to the other.

To show the principle of the GMM, let us again consider the layered system of Figure 2.4. Assuming layers of the same type, we can use the continuity of the state vector components on the first interface:

[M [0] (I 1 )]a [0] = [M [1] (I 1 )]a [1] .

(2.32)

This can be also stated as

[[M [0] (I 1 )][-M [1] (I 1 )]] a [0] a [1] = 0. (2.33)
We can repeat the process for all interfaces and assemble the equations into a global matrix:

     [M [0] (I 1 )] [-M [1] (I 1 )] 0 • • • 0 0 0 [M [1] (I 2 )] [-M [2] (I 2 )] • • • 0 0 . . . . . . . . . . . . . . . . . . 0 0 0 • • • [M [n-1] (I n )] [-M [n] (I n )]      •           a [0]
a [1] a [2] . . .

a [n-1] a [n]           = 0.
(2.34)

If adjacent layers are not of similar media (fluid-porous for example), the state vectors are not of the same size and we have to again use interface matrices [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF].

Once the global matrix (2.34) is assembled, just like with the TMM, we have to know some of the wave amplitudes in order to solve the response of the system. Typically, we set one incoming wave amplitude to one and the rest to zero, and in the case of a layered system surrounded by fluid we then have a [0] = [1, R] T and a [n] = [0, T] T . The unknowns in the problem then consist of the amplitudes within the layers (a [1] , . . . , a [n-1] ), as well as of R and T. The known parts of a [0] and a [n] , along with the parts of the global matrix corresponding to them, can now be moved to the right hand side of equation (2.34), and the unknown amplitudes can be solved by a single matrix inversion.

Bayesian methods for parameter estimation

In this chapter, we discuss the inverse problem in the characterisation of porous media, namely estimating model parameters from acoustic measurements. Parameter estimation is a topic in publications I -III, but computational methods for the inversion are developed especially in publication II. We will treat the parameter estimation problem in the Bayesian framework, and this chapter starts by giving a brief review of Bayesian statistics. For a more comprehensive treatment we refer the reader to the following textbooks on Bayesian statistics and inversion: [START_REF] Calvetti | An introduction to Bayesian scientific computing: ten lectures on subjective computing[END_REF][START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF][START_REF] Kaipio | Statistical and Computational Inverse Problems[END_REF][START_REF] Gelman | Bayesian Data Analysis[END_REF][START_REF] Stuart | Inverse problems: A Bayesian perspective[END_REF]]. In the latter part of this chapter, we concentrate on methods of computational inference, mainly using sampling methods.

PARAMETER ESTIMATION IN THE BAYESIAN FRAMEWORK

In the Bayesian approach, all unknowns are modelled as random variables, and we express the information about them as probability densities. The solution of the inverse problem is the posterior probability density (ppd), which combines the information that we learn from measurement data and what we know from other sources (the prior). The ppd can then be used to answer all questions that we can pose in terms of probabilities. These include "What are the parameter values that maximise the posterior probability?", "What is the mean value of the posterior?", and "What are the intervals that include most of the posterior probability?".

The posterior is constructed using the famous Bayes' formula:

π(θ|y) = π(y|θ)π(θ) π(y) ∝ π(y|θ)π(θ), ( 3.1)
which is an expression of the conditional probability of the unknown parameters θ given the measurement data y. In order to find the posterior we have to consider the probability distribution of y in terms of θ, which is done by constructing a model for the likelihood π(y|θ). Likelihood tells us how the measurements would be distributed if we knew the unknown θ. In addition, we need to include the possible prior information on θ, denoted by the density π(θ), in which we express what we know about the unknowns prior to taking the measurement data into account.

Lastly, in formulation (3.1) the term π(y) is fixed, and only acts as a normalising constant. We can therefore ignore it and consider the unnormalised posterior that is proportional to the product of the likelihood and the prior.

If the inverse problem has more than two unknowns, visualising the full posterior is difficult. Typically we then summarise the information contained in the posterior with point and interval estimates. Two perhaps the most common statistical point estimates are the maximum a posteriori (MAP) and the conditional mean (CM) estimates. The MAP estimate is the point which maximises the posterior probability

θ MAP = arg max θ π(θ|y). (3.2)
Note that this maximiser may not exist or be unique [START_REF] Kaipio | Statistical and Computational Inverse Problems[END_REF]. Computing the MAP estimate requires the solution of an optimisation problem, and the procedure of computing the MAP estimates is similar to the one of the least squares.

The CM is defined as the expected value of the posterior, that is, the center-ofmass

θ CM = R n θ θπ(θ|y)dθ, ( 3.3)
where θ ∈ R n θ and n θ denotes the number of the unknowns. The conditional mean is a conditionally biased but unconditionally unbiased estimate [START_REF] Melsa | Decision and estimation theory[END_REF], and can be also shown to minimise the trace of the posterior covariance [START_REF] Kaipio | Statistical and Computational Inverse Problems[END_REF]. Computing the CM estimate is an integration problem. Usually the integral (3.3) is high-dimensional which renders quadrature integration approaches infeasible. In this case we typically approximate the integral using Markov chain Monte Carlo techniques, described later in this chapter.

Examples of the spread estimates include the conditional covariance

cov(θ|y) = R n θ (θ -θ CM )(θ -θ CM ) T π(θ|y)dθ, ( 3.4)
if it exists, and the Bayesian (marginal

) credible intervals CI k (p) = [a I , b I ] ⊂ R, defined as b I a I π(θ k |y)dθ k = p/100, (3.5)
where the marginal posterior density of θ k is defined as

π(θ k |y) = R n θ-1 π(θ|y)dθ 1 • • • dθ k-1 dθ k+1 • • • dθ n θ , (3.6)
and the interval CI k (p) contains p % of the marginal posterior mass. Definition of the p % credible interval is not unique. One way to define the interval is so that the tails of the density function contain an equal amount of posterior mass Credible intervals can be directly interpreted as statements on the probabilities of the parameter values [START_REF] Gelman | Bayesian Data Analysis[END_REF], i.e. they represent the uncertainty in the estimates.

It is also possible to define a credible set, where instead of the endpoints of the marginal densities we are looking for an equiprobability hypersurface in R n θ which encloses p % of the posterior mass.

Likelihood

Let us now describe the construction of the likelihood model. The likelihood is a function that describes the distribution of the measurements when we know what θ is, and what the measurements conditions are [START_REF] Calvetti | An introduction to Bayesian scientific computing: ten lectures on subjective computing[END_REF]. Likelihood thus contains the forward model, and information on measurement noise and some of the model uncertainties.

We will start by considering the measurement noise. Most often, as is done in this thesis as well, measurement noise is modelled as additive and mutually independent of the unknowns. Let y ∈ R n y , where n y denotes the length of the measurement vector. With the additive noise assumption, an observation model for the measurements can be written

y = h(θ) + e, (3.9)
where h : R n θ → R n y is the forward model that maps the unknown parameters to measurable data, and e ∈ R n y denotes the additive measurement noise component.

Because θ and e are mutually independent, the expression for the distribution of the measurements when θ is fixed can be stated as

π(y|θ) = π e (y -h(θ)), (3.10)
where π e (•) denotes the distribution of the noise. The noise is often assumed to be normally distributed, which means that the likelihood can be written as

π(y|θ) = 1 (2π) n y det Γ e exp - 1 2 (y -h(θ) -e * ) T Γ -1 e (y -h(θ) -e * ) , (3.11)
where Γ e is the noise covariance matrix, and e * the noise mean. If the noise covariance is known, the term in front of the exponential is constant and only acts as a normalising factor. In this thesis we assume that e * = 0 and that all covariances are positive definite. Then, by denoting the square root of the inverse of the noise covariance with L T e L e = Γ -1 e , we can write the likelihood as

π(y|θ) ∝ exp - 1 2 L e (y -h(θ)) 2 . (3.12)
If the noise covariance is unknown, we can model it as another random variable that we estimate, and write it inside the exponential term

π(y|θ, Γ e ) ∝ exp - 1 2 L e (y -h(θ)) 2 + log(det Γ e ) . ( 3

.13)

In publications I-III we model the noise as independent and identically distributed (i.i.d.), so that the general form of the noise covariance matrix is

Γ e = σ 2 e I, (3.14)
where σ 2 e is the noise variance and I denotes a n y × n y identity matrix.

Prior models

The prior density π(θ) is used to represent information about the unknowns that we have prior to taking the measurement data into account. In principle, constructing the prior is simple since the main requirement is that it gives a high probability to parameter values we expect to see, and a low probability to values we do not expect to see. The prior density should also constrain the unknowns to feasible values. Since the posterior is formed as a product of the likelihood and the prior, the prior has the biggest effect along the posterior dimensions that the data is not very informative. The exact analysis can be carried out using the eigenvalue analysis of posterior covariance. In parameter estimation this means that the prior is highly influential for parameters of which the data does not contain much information about.

One difficulty in designing the prior may arise when the prior knowledge is qualitative in nature, because it can be hard to translate into a quantitative representation.

The simplest form of the prior is a uniform density

π(θ) = 1, θ ∈ D 0, otherwise, (3.15)
where D is the domain (of feasible values) to which we want to restrict θ. We do not normalise the density because we are only interested in the relative probabilities between parameter values, and the constant normalising factor disappears in the proportionality (3.1). Priors of the form (3.15) are often used to specify positivity constraints, as was done in I, and they can also be used together with other types of priors, as was done in II and III.

Perhaps the most commonly used form of prior is the multivariate normal distribution, where the probability density is written as

π(θ) ∝ exp - 1 2 L θ (θ -θ * ) 2 . (3.16)
Here θ * denotes the expected value of the prior, and L T θ L θ = Γ -1 θ is a matrix square root of the inverse of the prior covariance matrix. We use multivariate normal priors (truncated with in indicator function of the form (3.15)) in publications II and III.

MARKOV CHAIN MONTE CARLO

Computation of the conditional mean and interval estimates requires integrating over the posterior density. The integrals are nearly always approximated numerically since an analytical representation of the posterior is rarely obtainable. In low dimensions we can use deterministic integration methods such as Gauss quadratures, but as the number of unknowns grows, the so-called "curse of dimensionality" renders such approaches infeasible. The alternative is to approximate the high-dimensional integrals by statistical methods, which are based on obtaining representative samples from the probability density of interest, and using the samples to approximate the expectations and other summaries. A common approach is Monte Carlo integration [START_REF] Robert | Monte Carlo statistical methods[END_REF], where the samples are drawn randomly around the parameter space. The basic Monte Carlo method can be made more efficient with methods such as importance sampling, where the samples are still drawn at random but more often in regions of high probability.

Markov chain Monte Carlo (MCMC) is a class of particularly useful sampling methods, where samples from the parameter space are obtained by simulating random walks in it [START_REF] Brooks | Handbook of Markov chain Monte Carlo[END_REF]. The random walk must form a Markov chain, which by definition means that the next step in the chain can depend only on the chain's previous position. The basic idea behind a MCMC simulation is that we start at a point (or in general, a state) X 1 in the parameter space, do a random change to it that satisfies certain conditions, and repeat enough times. In doing so we generate a large sample of points that are distributed according to a given target distribution π(•), such as the posterior. MCMC differs from ordinary Monte Carlo methods in that the samples are no longer independent. However, despite the dependence of the samples, running a MCMC simulation may be significantly more efficient, and often the only feasible solution, for integrating complicated probability densities.

A strength of the sampling methods is that once enough samples have been obtained, we can easily approximate unknown quantities of π(•), such as the expected value of a function g:

E π g = X g(θ)π(θ)dx ≈ 1 n n ∑ i=1 g(θ i ) , (3.17)
where X is the support of π(θ), and n is the sample size. For example the CM estimate is computed by taking the mean value of the obtained posterior samples. In addition, computing the posterior marginal densities and credible intervals from the samples is trivial. The reason we can make the approximation in (3.17) is the ergodic theorem, which implies that if Θ = (θ 1 , θ 2 , θ 3 , . . . ) is a Harris ergodic Markov chain with π(•) as its limiting distribution, then with probability 1,

ḡn := 1 n n ∑ i=1 g(θ i ) → E π g as n → ∞, (3.18)
if E π |g| < ∞ [START_REF] Flegal | Markov chain Monte Carlo: Can we trust the third significant figure?[END_REF]. For the precise definitions, see [START_REF] Robert | Monte Carlo statistical methods[END_REF][START_REF] Brooks | Handbook of Markov chain Monte Carlo[END_REF][START_REF] Gamerman | Markov chain Monte Carlo: stochastic simulation for Bayesian inference[END_REF]. To construct such a Markov chain we use a transition kernel K, which is a conditional probability density such that θ k+1 ∼ K(θ k , θ k+1 ). A simple example is the random walk process, θ k+1 = θ k + k , where the random term k does not depend on the previous samples. Obviously, when n < ∞ the approximation ḡn is not exact. We will discuss the error variance later in this chapter.

Metropolis-Hastings update

The transition kernel, or the update mechanism, with which the chain is generated is a central factor in the behaviour and efficiency of the MCMC sampler. The Metropolis-Hastings (MH) accept-reject update [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] forms the basis of most of the update mechanisms. Let us again denote the target probability density with π(•), and the current location of the chain in the parameter space by θ k . Then the steps for a single MH update are:

1. Propose a candidate θ for the next element in the chain by drawing it from the proposal density q(θ k , θ ).

Calculate the acceptance ratio

a(θ k , θ ) = min 1, π(θ )q(θ , θ k ) π(θ k )q(θ k , θ ) . (3.19)
3. With probability a(θ k , θ ) set θ k+1 = θ (accept proposal), and otherwise set θ k+1 = θ k (reject proposal).

Note that the chain is propagated even if the proposed candidate is rejected. In addition, when sampling from posterior densities, the normalising factor (π(y) in (3.1)) is cancelled out. The MH update is, in principle, very simple to implement. However, in practice there are enormous differences in the efficiency depending on the choice of the proposal density.

As stated above, a special case of the MH algorithm is the random walk Metropolis (RWM) algorithm, where the proposal is of the form θ = θ k + k , where k is stochastically independent of θ k . One of the most common method of constructing a RWM proposal is to use a Gaussian distribution, so the proposal has the form and k ∼ N (0, Σ), where Σ is the proposal covariance matrix. A useful property of the RWM proposals is that they are symmetric, i.e. q(θ k , θ ) = f ( θθ k ). With symmetric proposal densities the acceptance ratio (3.19) simplifies to

a(θ k , θ ) = min 1, π(θ ) π(θ k ) . (3.20)
This is called the Metropolis acceptance ratio and it is the original MCMC algorithm that was proposed in 1953 [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]. In the Metropolis algorithm, if π(θ ) is greater than or equal to π(θ k ) the proposal is always accepted, whereas if π(θ ) is smaller than π(θ k ) the proposal is accepted with the probability a(θ k , θ ).

With Gaussian proposal densities, a big question is how to choose the the proposal covariance Σ. As already noted, the proposal has a huge effect on the efficiency of the algorithm, and different ways to choose Σ has been an area of active research.

Efficiency measures

Although in the limit n → ∞, any correctly implemented MCMC sampler will generate a set of samples that are distributed according to the limiting distribution, an efficient sampler produces usable estimates quickly, whereas an inefficient one may not converge in the lifetime of the universe. We need therefore ways to assess the efficiency of the samplers.

Autocorrelation An efficient MCMC sampler will mix well, i.e. it moves quickly all around the relevant parts of the parameter space, and will achieve convergence quickly (as opposed to pseudo-convergence where the chain seems to have converged but is actually stuck in a local maximum). Unfortunately, no diagnostic can reliably detect the convergence of a Markov chain. What we can do, however, is to check the autocorrelation of the chain, which gives an indication on how many iterations we need to produce one independent sample, and thus says something about the efficiency of the sampler. The ideal Markov chain would have an integrated autocorrelation time of one, τ = 1, equivalent to i.i.d. sampling.

Samples in a Markov chain are roughly independent when we consider a subsampled chain, where we take every τ:th sample. To measure the efficiency of the sampler, we are mainly interested in the number of independent samples, or the effective sample size (ESS), of the chain. The ESS can be estimated by considering the autocorrelation time for one variable at a time, but this ignores the correlations between parameters and gives an inaccurate picture of the sample quality. Vats et al. [START_REF] Vats | Multivariate output analysis for Markov chain Monte Carlo[END_REF] have proposed a multivariate version of the ESS, which considers the covariance structure of the sample. They define it as

mESS = n |Λ| |Σ ∞ | 1/n θ , (3.21)
where n is the sample size, n θ is the dimension of the parameter space, Λ is the sample covariance matrix, Σ ∞ is the asymptotic covariance matrix in the Markov chain central-limit theorem (see eq. (3.24)), and | • | represents the matrix determinant. An estimate Σ n for Σ ∞ can be found by a multivariate batch means estimator:

Σ n = b n a n -1 a n -1 ∑ k=0 ( θk -θ)( θk -θ) T . (3.22)
To calculate Σ n we need to cut the chain into a n batches that have the length b n (i.e. n = a n b n ), θk is the mean value of the batch k:

θk = 1 b n kb n ∑ i=(k-1)b n +1 θ i , (3.23)
and θ is the mean value of all of the samples (i.e. the CM estimate). There are several ways to decide the batch size, but in general it should grow with the chain length, for example as b n = n 1/2 or b n = n 1/3 . Also higher correlation in the parameters calls for a larger batch size, suggesting that b n ∝ τ.

Stopping criteria and diagnosing convergence

We need to have some criteria to stop the MCMC simulation. The simplest one is to run the chain for N predetermined steps and then stop, but this is far from ideal.

To make meaningful inferences we need to be relatively sure that the chain has first of all converged and secondly had enough time to sample the target distribution.

Autocorrelation times together with looking at the time series of the individual chains provide a tool to assessing whether we can stop the simulation or not.

A more objective stopping criteria can also be constructed. For example, we can decide to stop the simulation once a predetermined accuracy of a statistic is reached. Commonly we look at the estimate of the mean, θ, and its estimation error, the Monte Carlo standard error (MCSE), defined in the following.

Burn-in

Before it is sensible to start estimating properties such as the mESS and MCSE, we need to make sure the chain has reached the steady-state. In the beginning of the MCMC run, especially if the starting location of the chain was far away from the main support of the pdf, the values do not represent the target density very well and it can be good to remove them so that they do not cause bias in the estimates. This initial period is called the burn-in period. Removing the burn-in is important especially when we are using adaptive MCMC methods that we will describe later. However, note that removing a burn-in phase from the beginning of the simulation is sometimes completely unnecessary [START_REF] Brooks | Handbook of Markov chain Monte Carlo[END_REF]. This is the case when we do not use adaptation, and do sufficiently long runs so that the effect caused by the initial locations is negligible. During the first 15,000 iterations the chain is drifting towards the region of the main probability mass, but after that it stays more or less stationary. Figure 3.1 shows an example of a MCMC run, where the chain was not started in the region where the most of the probability mass is, and therefore in the beginning we can see the chain drifting towards this area. After about 15,000 iterations, the drifting seems to stop and we have reached the steady state. We could then remove these 15,000 iterations from the beginning as the burn-in.

Monte Carlo standard error

The difference ḡn -E π g is called the Monte Carlo error, and although in practice we can never know its true value, we can approximate it by estimating the variance of its asymptotic distribution. Let us examine the unidimensional case, i.e. θ is a single unknown. The central limit theorem (CLT) states that [START_REF] Flegal | Markov chain Monte Carlo: Can we trust the third significant figure?[END_REF]:

√ n ( ḡn -E π g) → N 0, σ 2 g , (3.24)
as n → ∞, where σ 2 g := var π {g(θ 1 )} + 2 ∑ ∞ i=2 cov π {g(θ 1 ), g(θ i )}, and we assume that θ 1 , θ 2 , θ 3 , . . . are identically distributed and θ 1 is distributed according to the target distribution π(•). The CLT can be stated in the multivariate setting too, where we replace σ 2 g with its multivariate counterpart Σ g . The CLT shows that that the accuracy of a Monte Carlo simulation increases as √ n. To estimate the Monte Carlo error we need an estimator σ2 g for the asymptotic variance σ 2 g . Non-overlapping batch means are one of several ways to do this. If we denote the mean of the k:th batch [g(θ (k-1)b n +1 ), ..., g(θ kb n )] with Ȳk (see Eq. (3.23)), the batch means estimate for σ 2 g is

σ2 g = b n a n -1 a n ∑ k=1 ( Ȳk -ḡn ) 2 , (3.25)
where again a n is the number of batches and b n is the batch length. Now we can estimate the Monte Carlo standard error of ḡn [START_REF] Flegal | Markov chain Monte Carlo: Can we trust the third significant figure?[END_REF]:

MCSE = σg √ n , (3.26)
and we can also calculate the (1α) • 100 % confidence interval for the mean. We know that the error of the mean has a t-distribution with n -1 degrees of freedom (the CLT tells us that asymptotically, the error is normally distributed). So we are looking for values q 1 and q 2 that contain the true value with probability (1α).

From probability theory we know that [100]

q 1 < ḡn -E π g σg / √ n < q 2 iff. ḡn -q 2 σg √ n < E π g < ḡn -q 1 σg √ n , (3.27)
where

P q 1 < ḡn -Eg σg / √ n < q 2 = 1 -α . (3.28)
For example the half width of the (1α) • 100 % confidence interval for Eg is then given by

CI 1-α = t a n -1 1 - α 2 σg √ n , (3.29)
where t a n -1 (•) denotes the inverse cumulative distribution function of the t-distribution with a n -1 degrees of freedom.

Stopping criteria

Let us now examine a few criterion for stopping the simulation.

Fixed-length

The standard way to run MCMC is to iterate for N steps, check visually if the run has "converged" and either stop or run for longer. This is not very accurate, though, and better ways have been proposed.

Fixed-width

In the fixed-width criterion the desired level of accuracy is predefined and the simulation is stopped once it is reached. One example of such a stopping rule is given in [START_REF] Flegal | Markov chain Monte Carlo: Can we trust the third significant figure?[END_REF], where the stopping criterion is

t a n -1 σg √ n + εI(n < n * ) < ε , (3.30)
where ε is the maximum error level, I is an indicator function, and the role of n * is to prevent a premature stopping due to a bad estimate of σg . In practice, this is calculated for each parameter separately and the simulation is terminated when the uncertainties in all of the parameters are under the error limit.

Fixed-width relative to standard deviation Another way to stop the simulation is to base the maximum error level ε on a percentage of the standard deviation of the parameters, instead of a predefined level of accuracy. This allows for some flexibility in the stopping criterion since the method allows for more absolute error in parameters that have inherently more variability while keeping the relative error the same. This is the method we chose to use in publications II and III. The stopping criterion is

t a n -1 σg √ n + ε SD I(n < n * ) < ε SD , (3.31)
where ε SD represents now the allowed percentage error in the standard deviation σ θ of the actual parameter, for example ε SD = 0.10 • σ θ .

Multivariate fixed-width It can be shown [START_REF] Vats | Multivariate output analysis for Markov chain Monte Carlo[END_REF] that stopping based on the effective sample size is equivalent to the fixed-width relative to standard deviation rule. Calculating the multivariate ESS would take the correlations between parameters into account, making the stopping rule even more robust.

MCMC SAMPLERS

In this section, we will describe some commonly used updating mechanisms for the Markov chain. We call the update mechanisms 'samplers' since they determine how the MCMC code samples the target density.

A robust sampler is one that can deal with probability densities that are challenging to sample, such as multimodal or strongly non-normal ones. We need a robust sampler so we can be reasonably confident that we are indeed sampling from the full distribution instead of just a local maximum. Unfortunately, the robustness of a sampler also depends on the application, and hence universal recommendations that work in all cases are not possible.

Overview

Since the "second-generation MCMC revolution" starting approximately in the 1990's [START_REF] Brooks | Handbook of Markov chain Monte Carlo[END_REF], a huge number of MCMC samplers have been proposed. These include adaptive Metropolis [START_REF] Haario | An adaptive Metropolis algorithm[END_REF][START_REF] Andrieu | A tutorial on adaptive MCMC[END_REF][START_REF] Haario | Adaptive proposal distribution for random walk Metropolis algorithm[END_REF][START_REF] Atchadé | On adaptive Markov chain Monte Carlo algorithms[END_REF][START_REF] Andrieu | On the ergodicity properties of some adaptive MCMC algorithms[END_REF][START_REF] Roberts | Coupling and Ergodicity of Adaptive Markov Chain Monte Carlo Algorithms[END_REF][START_REF] Roberts | Examples of Adaptive MCMC[END_REF][START_REF] Vihola | Robust adaptive Metropolis algorithm with coerced acceptance rate[END_REF], Metropolis-adjusted Langevin [START_REF] Roberts | Exponential Convergence of Langevin Distributions and Their Discrete Approximations[END_REF][START_REF] Roberts | Optimal scaling of discrete approximations to Langevin diffusions[END_REF][START_REF] Atchadé | An Adaptive Version for the Metropolis Adjusted Langevin Algorithm with a Truncated Drift[END_REF], Hamiltonian Monte Carlo [START_REF] Duane | Hybrid Monte Carlo[END_REF][START_REF] Neal | Bayesian Learning for Neural Networks[END_REF][START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Hoffman | The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo[END_REF], Gibbs [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF][START_REF] Gelfand | Sampling-Based Approaches to Calculating Marginal Densities[END_REF], reversible jump MCMC [START_REF] Dettmer | Trans-dimensional geoacoustic inversion[END_REF][START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF][START_REF] Green | Trans-dimensional Markov chain Monte Carlo[END_REF][START_REF] Brooks | Efficient construction of reversible jump Markov chain Monte Carlo proposal distributions[END_REF][START_REF] Bodin | Seismic tomography with the reversible jump algorithm[END_REF][START_REF] Bodin | Transdimensional tomography with unknown data noise[END_REF], delayed rejection [START_REF] Mira | On Metropolis-Hastings algorithms with delayed rejection[END_REF][START_REF] Haario | DRAM: efficient adaptive MCMC[END_REF], randomize-then-optimize [START_REF] Bardsley | Randomize-Then-Optimize: A Method for Sampling from Posterior Distributions in Nonlinear Inverse Problems[END_REF][START_REF] Bardsley | Randomize-Then-Optimize for Sampling and Uncertainty Quantification in Electrical Impedance Tomography[END_REF], parallel tempering [START_REF] Earl | Parallel tempering: Theory, applications, and new perspectives[END_REF][START_REF] Vousden | Dynamic temperature selection for parallel tempering in Markov chain Monte Carlo simulations[END_REF][START_REF] Swendsen | Replica Monte Carlo simulation of spinglasses[END_REF][START_REF] Geyer | Markov chain Monte Carlo maximum likelihood[END_REF][START_REF] Miasojedow | An Adaptive Parallel Tempering Algorithm[END_REF][START_REF] Araki | Adaptive Markov chain Monte Carlo for auxiliary variable method and its application to parallel tempering[END_REF][START_REF] Sambridge | A Parallel Tempering algorithm for probabilistic sampling and multimodal optimization[END_REF], affine invariant ensemble sampling [START_REF] Christen | A general purpose sampling algorithm for continuous distributions (the t-walk)[END_REF][START_REF] Goodman | Ensemble samplers with affine invariance[END_REF][START_REF] Foreman-Mackey | emcee: The MCMC Hammer[END_REF], elliptical slice sampling [START_REF] Murray | Elliptical slice sampling[END_REF], stochastic Newton [START_REF] Martin | A Stochastic Newton MCMC Method for Large-Scale Statistical Inverse Problems with Application to Seismic Inversion[END_REF][START_REF] Petra | A Computational Framework for Infinite-Dimensional Bayesian Inverse Problems, Part II: Stochastic Newton MCMC with Application to Ice Sheet Flow Inverse Problems[END_REF], and likelihood-informed subspace [START_REF] Cui | Likelihoodinformed dimension reduction for nonlinear inverse problems[END_REF][START_REF] Constantine | Accelerating Markov Chain Monte Carlo with Active Subspaces[END_REF] algorithms. This list is by no means exhaustive, and in the following we will only describe the samplers that were used in I-III. ∼ U ([-σ, σ]), where Σ denotes the proposal covariance matrix, and σ a vector of proposal widths. However, other types of proposals have also been successfully used [START_REF] Yang | Searching for efficient Markov chain Monte Carlo proposal kernels[END_REF]. In the following we will use the term "proposal (co)variance" to mean generally the width of the proposal density. The RWM algorithm is very simple to implement, but can be difficult to tune properly, the reason being that no matter what kind of proposal we are using, we have to fix its variance before the simulation.

Random walk Metropolis and acceptance rate

If the variance of the proposal is set too large, the proposals θ are far away from the current location θ k , and likely the value of the density at π(θ ) is much lower than at π(θ k ). Then the acceptance ratio (3.20) is small, the proposal will be rejected often, and the chain will not move. This wastes computation resources and leads to long autocorrelation times. On the other hand, if the proposal variance is set too small, virtually all proposals will be accepted but the chain still moves very slowly around the parameter space. Again we have long autocorrelation times. Clearly there must exist better values for the proposal variance between these two extremes, where the proposals are sometimes accepted and sometimes rejected, and the chain moves around the parameter space more efficiently. However, finding the efficiencywise optimal, or even a good, proposal covariance by trial and error when little is known about the target distribution, is problematic to say the least.

This behaviour of the chain can be monitored with the acceptance rate, which is defined simply as the long-term average of the acceptance ratios (3.20). Again, a general optimal value for the acceptance rate does not exist, but as a rule of thumb as long as the acceptance rate is between 10 and 60 % the extremes of bad sampling are avoided. However, it has been shown [START_REF] Gelman | Efficient Metropolis jumping rules[END_REF][START_REF] Roberts | Weak convergence and optimal scaling of random walk Metropolis algorithms[END_REF] that efficiency-wise the optimal acceptance rate for uncorrelated Gaussian target distributions is 0.44 in the one-dimensional case, declining to 0.23 when n θ → ∞. In addition, the highest efficiency is achieved when the proposal covariance mimics the covariance of the target distribution, scaled by approximately 2.4 2 /n θ .

Adaptive algorithms

Even if the size and spatial orientation of the proposal covariance could be tuned by trial and error in low-dimensional cases, in high dimensions it is virtually impossible. Fortunately, the proposal covariance can be adapted automatically as the simulation goes on and we learn more about the target distribution.

One of the first adaptive algorithms was the adaptive Metropolis (AM) [START_REF] Haario | An adaptive Metropolis algorithm[END_REF][START_REF] Haario | Adaptive proposal distribution for random walk Metropolis algorithm[END_REF], where it was proposed to adapt the size and orientation of the proposal distribution on the fly as the simulation progresses. The basic idea is to approximate the target density with a multivariate Gaussian and use the approximation as the proposal distribution:

θ = θ k + k , k ∼ N (0, Σ k ) , (3.32)
where

Σ k = Σ 0 , k ≤ k 0 s d cov (θ 1 , ..., θ k ) + s d εI, k > k 0 . (3.33)
We need to select a starting proposal covariance Σ 0 , but usually it is enough to choose something small so that proposals are accepted and a first estimate of the target distribution covariance can be calculated after k 0 iterations. The scaling parameter s d is typically set to s d = 2.4 2 /n θ , which is a good choice especially if the target density is close to Gaussian [START_REF] Gelman | Efficient Metropolis jumping rules[END_REF]. The parameter ε is introduced to prevent the covariance Σ k from becoming a singular matrix, and is selected to be very small compared to the size of the support of the target density. The proposal covariance Σ k can be updated iteratively, so that computing it induces practically no extra computational cost.

It is clear that with adaptive methods the resulting chain is no longer strictly Markovian, since θ n now depends on all the previous states. However, it is shown in [START_REF] Haario | An adaptive Metropolis algorithm[END_REF] that under certain conditions the chain is still ergodic and the adaptive algorithm can be used. Later it was realised that the AM algorithm can be viewed as belonging to a more general set of stochastic control algorithms [START_REF] Andrieu | Controlled MCMC for optimal sampling[END_REF] and convergence of more general adaptive algorithms was also proved [START_REF] Atchadé | On adaptive Markov chain Monte Carlo algorithms[END_REF][START_REF] Andrieu | On the ergodicity properties of some adaptive MCMC algorithms[END_REF]. The main requirement for the adaptive algorithms to be ergodic is that the adaptation is "vanishing", i.e. that for every successive iteration we adapt less and less.

The AM algorithm conditions hold also for a case where we do not consider all the past states of the chain, but only an increasing part of the near history. Calculating the covariance from (θ n/2 , θ n/2+1 , ..., θ n ) instead of all the samples can be better in some situations. This can be the case for example when the starting location is bad and the first samples would orientate the covariance unfavourably.

Parallel tempering

The AM algorithm works very well in unimodal and convex cases, even if there are strong correlations between the parameters. However, it can be slow to converge and travel around the parameter space in more complicated situations. Multimodal distributions are notoriously hard to sample efficiently, because the individual peaks are usually separated by large areas of almost zero probability. Even if the target density is not multi-modal as such, it can include narrow peaks of high probability where the chain is likely to get stuck for long periods of time, especially if using adaptive methods.

Parallel tempering is one way to make sampling such distributions easier [START_REF] Earl | Parallel tempering: Theory, applications, and new perspectives[END_REF][START_REF] Swendsen | Replica Monte Carlo simulation of spinglasses[END_REF][START_REF] Geyer | Markov chain Monte Carlo maximum likelihood[END_REF]. It is based on raising the target density to a power, i.e. tempering, as

π(θ) 1/T
where T can be interpreted as the temperature. At high temperatures the possible peaks of the density are flattened out, making it easier for the chain to travel between them. For a simple example with a one-dimensional density, see Fig. 3.2. At T = 1 the target density consists of two peaks separated by a wide zero probability area. This is a case where it would be difficult to construct an efficient sampler, because to sample within the peaks the proposal needs to be narrow but to jump between the peaks it needs to be wide. When the temperature is increased, the peaks start to flatten out. At T = 100 the modes have mixed into one, making the density easier to sample.

Naturally, tempering changes the target density into something else that we are not primarily interested in. The idea in parallel tempering is to sample in parallel n T Markov chains at temperatures 1 = T 1 < T 2 < • • • < T n T , so we always keep sampling the density we are interested in, but also densities in higher temperatures. Then we construct a mechanism that couples the different temperature chains together, so they can exchange information on their location in the parameter space. Chains in densities with higher temperatures can in effect tell the chains sampling the colder densities where they should move. In practice, this exchange of information is done by swapping the locations of two chains.

Between-chain Metropolis update

Swapping the locations must preserve the distributions of the respective tempered chains, and this is achieved by once again using the Metropolis rule (3.20). This gives a general name for parallel tempering and other similar auxiliary variable methods, Metropolis-coupled Markov chain Monte Carlo, abbreviated sometimes as MCMCMC or MC 3 .

The Metropolis acceptance ratio can be derived as follows in the between-temperature swapping case. First, let us define the inverse temperature β i ≡ 1/T i , and consider the joint density of two chains, π(θ i , θ j ), that are at temperatures β i and β j :

π(θ i , θ j ) = π(θ i ) β i π(θ j ) β j , (3.34)
where θ i is the current location of the Markov chain that samples from the distribution π(θ i ) β i (respectively θ j for π(θ j ) β j ). Proposing to swap the locations is equal to proposing new values (θ i , θ j ) = (θ j , θ i ). Then we can write the Metropolis ratio:

A j,i = π(θ i , θ j ) π(θ i , θ j ) = π(θ j , θ i ) π(θ i , θ j ) = π(θ j ) β i π(θ i ) β j π(θ i ) β i π(θ j ) β j , ( 3.35)
and the acceptance probability of the swap is a = min(1, A j,i ). Often for numerical reasons we want to deal with the logarithms of the proposal densities. In that case the Metropolis ratio simplifies to

A j,i = log π(θ j ) π(θ i ) β i -β j . ( 3.36)
In Bayesian inference, we often do not want to temper the whole posterior density but only the likelihood:

π T (θ|y) ∝ π(y|θ) 1 T π(θ).
(3.37)

Then if we set T n T = ∞, at the highest temperature we are sampling only the prior, which is often easy to sample. In addition, values outside the prior support we have defined as infeasible, so restricting the hottest chain to search in the prior volume is a natural limit for the space. The Metropolis ratio for the case (3.37) reduces to

A j,i = log π(y|θ j ) π(y|θ i ) β i -β j , (3.38)
because the prior is cancelled out. In principle, swaps can be proposed between any two chains, but in practice it is enough to propose to swap adjacent chains since the swap acceptance ratio drops quickly with a growing temperature difference.

Temperature ladder

An essential factor affecting the efficiency of parallel tempering algorithms is the number and distribution of the temperatures T 1 < T 2 < • • • < T n T (the temperature ladder). Analogously to the MH acceptance rate issue, if the temperatures are too close to each other the distributions are nearly the same and thus almost all swaps are accepted. But there is no advantage to sampling the same distribution, so we want to spread the temperatures out. On the other hand, if the temperature difference is too big, none of the swaps are accepted and we are simulating multiple auxiliary chains for nothing. It turns out that a similar acceptance rate as the withinchain rate (23 % in the asymptotic limit) is optimal for the between-chain swaps as well [START_REF] Kone | Selection of temperature intervals for paralleltempering simulations[END_REF]. Then our goal should be to first make sure the swapping rate between adjacent chains is equal for all pairs, and close to 23 %. A related topic is the maximum temperature we should use, if we do not use the tempered likelihood form (3.37) where it can be set to infinity. The higher T n T is, the higher is the number of temperatures needed which increases the computational cost. However, for the sampling to be efficient T n T must be high enough so that the peaks have flattened out.

It is not possible to choose a perfectly efficient temperature ladder before starting the simulation because we do not know the type of the target distribution. A good starting point for the temperature ladder is to distribute the temperatures geometrically [START_REF] Earl | Parallel tempering: Theory, applications, and new perspectives[END_REF][START_REF] Kofke | On the acceptance probability of replica-exchange Monte Carlo trials[END_REF]. This provides a nearly constant acceptance rate for swaps in many systems. However, big changes may be still needed before the temperature ladder works optimally, and doing it by trial and error is very time-consuming. One option is to adapt the temperature values automatically during the simulation, where T 1 = 1 and T n T are kept constant and the rest are adapted so that the acceptance rate for between-chain swaps is equal between all temperatures. This maximises the communication efficiency between the chains. For examples on the temperature adaptation, see [START_REF] Vousden | Dynamic temperature selection for parallel tempering in Markov chain Monte Carlo simulations[END_REF][START_REF] Miasojedow | An Adaptive Parallel Tempering Algorithm[END_REF].

Review of thesis results

In this Chapter, we review the main results in publications I-III. Methods that we used to solve the forward and inverse problems in these publications were discussed in Chapters 2 and 3. We will therefore concentrate here mainly on describing the measurement setups and findings based on the parameter estimates.

PUBLICATION I: CHARACTERISATION OF RIGID FRAME POROUS MATERIALS FROM IMPEDANCE TUBE MEASUREMENTS

In publication I, we studied the characterisation of air-saturated porous media using the Bayesian approach, and compared the results to a least squares method.

The main contribution of the paper was to provide a robust and easy-to-use characterisation method for rigid frame materials. In particular, we computed parameter uncertainty estimates and examined correlations between the parameters. The measurements were done in a standard transmission impedance tube (see for example [START_REF] Song | A transfer-matrix approach for estimating the characteristic impedance and wave numbers of limp and rigid porous materials[END_REF]).

Measurement setup

We selected five materials, three polyurethane foams and two wools, that are commonly used, for example in sound and heat insulation applications in buildings and cars, see Fig. 4.1. Naturally, the proposed characterisation method extends beyond these materials and can be used with any air-saturated media that we can measure in an impedance tube and model as an equivalent fluid. These include glass beads, sandy soils and gravel, or even living plants [START_REF] Horoshenkov | A review of acoustical methods for porous material characterisation[END_REF]. The impedance tube where the samples were tested in is shown in Fig. 4.2. The sound source is on the left, the sample in the middle, and there are two microphones on both sides of the sample. The tube is closed off by a semi-anechoic ending. By sending a logarithmic sine sweep, and using the scattering matrix formalism, we measured the waves reflected and transmitted through the tested sample at frequencies 300-5800 Hz. The reflection and transmission coefficients were then converted into the equivalent density ρmeas eq and the equivalent bulk modulus Kmeas eq of the material.

Inversion

As mentioned earlier, we represented the porous materials using the JCAL model, which gives us directly the density ρmodel eq and bulk modulus Kmodel eq . We write the observation model as

y = h JCAL (θ) + e, (4.1)
where θ = [φ, α ∞ , Λ, Λ , k 0 , k 0 ], and y = [ ρ meas eq (ω), Kmeas eq (ω)] T . The cost functional for the least squares solution was formed as

f (ω; θ) = W 1 ρ meas eq (ω) -ρ model eq (ω; φ, α ∞ , Λ, k 0 ) 2 + W 2 Kmeas eq (ω) -Kmodel eq (ω; φ, Λ , k 0 ) 2 , (4.2)
where the weights W 1 = 1/ρ 2 RMS and W 2 = 1/K 2 RMS are added to compensate for the different amounts of influence the measurements of the density and bulk modulus have on the estimate of φ. Note that the viscous and thermal effects are decoupled when we form the minimization problem this way, and only porosity is found in both expressions. The LS solution is then obtained as

θ = arg min θ∈D f (ω; θ), (4.3)
where D is the domain to which we want to restrict the parameter values. In this paper, the only restriction is to bound the values of the unknowns to physically reasonable ones, i.e. φ ∈ [0, 1], α ∞ ∈ [1, 10], and so on. The minimisation was carried out using the Nelder-Mead algorithm [START_REF] Lagarias | Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions[END_REF]. We found that with some materials the algorithm often converged to different minima depending on the initial values. To increase the chances of finding the global minimum we used a multistart approach, where the minimisation was started several times with the initial values drawn randomly around the parameter space.

Let us now describe the construction of the prior and likelihood densities so that we can compute the Bayesian estimates. To achieve results comparable to those given by the LS solution, we use a uniform prior that provides the same physical constraints as in the LS case. We assume that the measurement is corrupted by white noise, with the covariance

Γ e = Γ 1 0 0 Γ 2 , ( 4.4)
where Γ 1 = σ 2 e,ρ ρ 2 RMS I and Γ 2 = σ 2 e,K K 2 RMS I. The quantities σ e,ρ and σ e,K denote the standard deviations of the measurement error, estimated from repeated measurements, and I represents an identity matrix. The posterior density can then be stated as

π(θ|y) ∝ exp - 1 2 L e (y -h JCAL (θ)) 2 π unif (θ), (4.5)
where L e L T e = Γ -1 e . We used the adaptive Metropolis algorithm to compute statistics of the posterior.

Results

The fit of the JCAL model (with the MAP estimate as the input) to measurements of the density and bulk modulus of two example materials is shown in Fig. 4.3. The largest model discrepancy is observed at the solid frame resonance frequency (∼ 1 kHz for the foam and ∼ 2 kHz for the wool), where the coupling between the air and the frame is higher and the frame vibrates in response to the acoustic excitation. The effect of this model discrepancy on the inversion is relatively small, however, because the frequency range where the resonance effects appear is narrow and the model errors exhibit symmetry around the reconstructed curves. We observed similar model fits to all measured samples (see publication I). dimensional marginal distributions and the MAP and LS point estimates. The two point estimates coincide very accurately, as expected, and are always found within the 50 % regions of distribution. We can observe a strong positive correlation between porosity φ and the thermal characteristic length Λ . However, both parameters are still resolved well as shown by the narrow support of the posterior distribution.

The posterior densities of all the foam materials showed only one mode. They were thus easy to sample and we did not encounter identifiability issues. The posterior densities of the wool samples were still mostly unimodal, but especially the thermal parameters (Λ and k 0 ) were not as well resolved as with the foams. One reason for this is that the wool samples showed more anisotropy, and consequently

PUBLICATION II: COMPUTATIONAL METHODS FOR ESTIMAT-ING THE MATERIAL PARAMETERS OF POROELASTIC MEDIA

In the second publication we examined the estimation of the parameters of a poroelastic plate from simulated ultrasound measurements. The plate was assumed to be submerged in water so that acoustic waves can excite the frame of the material.

The motivation for this study was to develop statistical inversion methods that work with the strongly non-linear poroelastic models, design a simple measurement that could be carried out in practice, and assess what kind of information can be extracted from the this kind of measurements. A numerical study is needed so that we can compare the estimates and the true parameter values to for example evaluate possible estimation biases.

The simulation model

We considered a measurement setup up shown in Fig. 4.5, where Ω [1] denotes the poroelastic layer, and Ω [0] and Ω [2] are semi-infinite halfspaces of water. The measurement data we wanted to obtain were the acoustic reflection and transmission coefficients R exp and T exp , respectively. To simulate the measurements while avoiding the inverse crime [START_REF] Kaipio | Statistical and Computational Inverse Problems[END_REF], we used a finite-element model constructed in COMSOL Multiphysics R software v. 5.3. We truncated the computational domain by the following boundary conditions: a plane wave radiation condition on the left side (Γ 1 ) to create the incident plane wave, an absorbing impedance boundary condition on the right (Γ 2 ), and a periodic Floquet condition on the top and bottom to simulate a geometry extending to infinity in the x 1 direction. We considered frequencies between 100 and 400 kHz, with 5 kHz steps, and added Gaussian noise to the final measurements. The simulation was repeated for angles between 0 • and 45 • , with 5 degree steps.

Inversion

Our goal was to estimate all parameters included in the Biot model that are related to the poroelastic plate (i.e. we assumed that the parameters of the saturating fluid are known). The real and imaginary parts of the elastic parameters were estimated separately, and we also estimated the angle of the incoming wave. The vector of all unknowns thus reads θ = [Λ, α ∞ , log 10 (k 0 ), φ, Re K b , Im K b , Re K s , Im K s , Re N, Im N, ρ s , L, ϕ inc ].

We expressed permeability on a logarithmic scale because of the huge range of values it can take. In addition, we estimated the measurement noise level(s) simultaneously with the material parameters. We formed the measurement vector by stacking the reflection and transmission coefficients together as y(ϕ inc ) = [R exp (ϕ inc ), T exp (ϕ inc )] T in the same vector, and we considered inversion from one angle at a time. The measurement noise was again assumed to be additive, and the noise covariance read

Γ e = σ 2 e R I 0 0 σ 2 e T I , (4.6)
where σ 2 e R and σ 2 e T are the noise variances in the reflection and transmission measurements, respectively. We can now denote the full vector of unknowns by θ = [θ, σ e R , σ e T ].

The prior was constructed using a multivariate Gaussian density, π( θ) ∼ N ( θ * , Γ θ), whose tails were truncated using an indicator function B( θ):

B( θk ) = 1, if θmin,k ≤ θk ≤ θmax,k 0, otherwise. (4.7)
The indicator function was used to enforce hard boundaries such as positivity or negativity constraints, as well as the physical constraint K b ≤ K s (1φ) [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF]. The final step in the inverse problem was to select a suitable MCMC sampler to calculate posterior statistics. We could not make the AM algorithm, used in publication I, work satisfactorily in the present problem, because the sampler had a tendency to get stuck in a local maximum for long periods of time, which could lead to false diagnosis of convergence if we did now know the true parameter values. A solution we found to work very well was to use parallel tempering, especially of the form (3.37) where we could set the highest temperature to infinity. Especially when the PT algorithm was combined with an adapted proposal density and adapted temperature ladder, the sampler was easily able to jump from any starting point within the prior support to sampling close to the true parameter values.

An example of the simulated noisy measurement is shown in Fig. 4.6. The figure also shows the (noiseless) model prediction that corresponds to the CM estimate, and the posterior predictive distribution. We can observe that the measurement points are all found within the ±3 standard deviations of the posterior predictive distribution, which suggests that the noise levels were estimated correctly.

Results

The most interesting measurement to consider is the one at normal incidence, since it is simpler to carry out than measurements at oblique incidences. We found that many parameters can be identified from the normal incidence data, in the sense that the posterior variance is much smaller than the prior variance. These parameters include porosity φ, tortuoisity α ∞ , viscous characteristic length Λ, and density of the solid ρ s .

The marginal posterior distributions of some of the parameters are shown in Fig. 4.7, which also shows the marginal prior density, point estimates, and the 95 % credible intervals. An interesting observation is that while the marginal distribution of permeability k 0 has a peak at the location of the true value, the distribution is highly skewed towards higher values of permeability. It seems that there is a clear lower bound for k 0 , but once past a certain value, changes in permeability do not show in the data any more. Apart from k 0 , the marginal posterior densities seem to be centred around the true parameter value, and the CM estimate tends to be very close to the truth.

The joint marginals posterior distributions reveal more information on the parameters and their identifiability. rameters that exhibited the most correlation. The strongest correlation is between the real parts of the bulk and the shear moduli, which are almost exactly linearly negatively correlated, and both parameters admit values from zero to way past the true values. However, the true values of the parameters are on the long and narrow support of the posterior which shows that the parameter estimates and their uncertainty bounds still contain the correct value. Another parameter that exhibits strong correlation with the elastic moduli is α ∞ .

The strong correlation between K b and N is a result of how the elastic moduli affect the wave speeds. The speed of the fast longitudinal wave is controlled by a linear combination of K b and N (in addition to other factors such as density), whereas the only elastic modulus that controls the speed of the shear wave is N. Therefore we would need both types of waves present to resolve both moduli, which we do not have at normal incidence.

At oblique angles of measurement, the shear was is generated, and consequently we are able to estimate K b and N with much better accuracy. This can be clearly seen in Fig. 4.9, which shows some of the parameter estimates, along with the estimate of the 95 % credible interval, as a function of the measurement angle. Unsurprisingly, the estimates of K b and N get better as the angle increases, but also the estimate of α ∞ improves with a growing angle of incidence. This may be due to the fact that α ∞ is correlated with both elastic moduli at normal incidence, as Fig. 4.8 shows, and can be more accurately estimated as the estimates of the elastic moduli are more accurate. However, the accuracy of the estimates for parameters such as porosity and solid density does not change significantly with the incidence angle.

PUBLICATION III: ESTIMATING THE MATERIAL PARAMETERS OF POROELASTIC OBJECTS FROM ULTRASOUND MEASURE-MENTS IN WATER

In publication III, we studied the applicability of the inversion method introduced in II to real measurements of a porous ceramic plate, a glass bonded silica, type QF-20 (Filtros Ltd., East Rochester, NY). The measurements were carried out in a water tank using focused ultrasound transducers, and in total 216 points on the plate were measured on a 80 × 230 mm regular grid, with 10 mm spacing between the measurement points. There were two main reasons for doing such a comprehensive set of measurements. First, estimating the material parameters independently at two or more locations very close to each other provides a sensibility check for the results. Large sudden jumps in parameter values could indicate a problem with the inversion and/or measurement method. Secondly, the resulting numerous parameter (and uncertainty) estimates can be plotted as a map that shows the potential variability of the parameters as a function of the measurement location. This is one way to quantify the inhomogeneity of the object.

The measurement setup

The frequency range of the measurements depends on the validity of the model and practical constraints for producing ultrasound in water. The effects of the Biot waves can only be seen below frequencies where scattering effects start to dominate, which we found for the studied material submerged in water to start at around 400 kHz. To produce low-frequency ultrasound, we used Panametrics V389 (Olympus Scientific Solutions Americas Inc., MA) point-focused transducers whose central frequency was 500 kHz. The usable frequency range of these transducers goes down to 200 kHz. The transducers were immersed in the water tank and connected to a computer-controlled measurement system. The porous plate was placed in the focal zone of the transducers, and the -6 dB beam diameter of a pulse/echo signal at 400 kHz at this distance was 10 mm. The emitting transducer was excited by a spike signal, which produced a short ultrasonic pulse in water. The received signals were then recorded on a computer, and 40 pulses were averaged for each measurement to reduce the effect of random noise.

We then computed the reflection and transmission coefficients of the plate, R exp and T exp , respectively, using the spectral ratio technique. First, the reflected and transmitted pulses are recorded using the two configurations shown in Fig. 4.10, and then a reference signal for both configurations is also recorded. For reflection, a single transducer is used in the pulse/echo mode and the reference signal is obtained by pointing the transducer upwards to the water-air interface. In the transmission configuration, two transducers are used and the reference signal is obtained in water. Finally, we take the Fourier transform to represent the measurements in the frequency domain, and use for example the Wiener deconvolution filter [START_REF] Chen | On effective spectrum-based ultrasonic deconvolution techniques for hidden flaw characterization[END_REF] to obtain R exp and T exp .

Inversion

The inversion method used here is similar to the one introduced in publication II, but a few modifications were made. In II, we modelled the visco-elastic losses in the solid with an imaginary part that was constant with respect to frequency. However, such an approach can be shown to be slightly non-causal [START_REF] Turgut | An investigation of causality for Biot models by using Kramers-Krönig relations[END_REF], which is why we implemented the constant Q model [START_REF] Kjartansson | ConstantQ-wave propagation and attenuation[END_REF] of Kjartansson to represent the losses with a causal model. To differentiate the elastic parameters from publication II, we denote their modulus here with a subscript of zero, like K b,0 , and the associated quality factor by Q K b (respectively for K s and N).

In addition to the plate's physical properties, we considered the measurement noise levels, thickness of the plate, and a distance mismatch term ε R , as nuisance parameters that we are not primarily interested in but nevertheless improve the inversion results. The term ε R is defined as the difference in the length of the signal path between the pulse/echo reflection measurement and the reference measurement, and has an effect on the phase of R exp .

The prior was constructed using truncated normal distributions, whose mean values were selected based on information that was available from the manufacturer of the plate and other studies. The prior variance was set to be large so that a wide range of parameter values were possible. In addition to the constraints employed in II using the indicator function, we included a constraint 2/3N 0 ≤ K b,0 which corresponds to requiring that the Poisson's ratio of the solid frame is positive.

Parameter inference at each measurement location was carried out independently of the other points.

Results

There is a large number of measurements and inversion results to examine, which is why we will summarise the results by showing only the point estimates and the width of their credible intervals. However, before presenting this summary let us concentrate on a randomly selected measurement location, and plot the model fit and marginal posterior distributions. This gives a more detailed view on the inversion, and lets us compare the findings in publication II.

Figure 4.11 shows the measurement of R exp and T exp from one of the measured locations, along with the model prediction that is computed using the CM estimate. The figure also shows the posterior predictive distribution using standard deviations. Perhaps the only clearly visible systematic discrepancy between the model and measurements can be seen at low frequencies of the transmission coefficient. Elsewhere the measurement points are found within plus or minus two standard deviations, and this was a trend observed in the other measurement locations as well.

Continuing with the same measurement location, the marginal posterior distributions of α ∞ , log 10 (k 0 ), φ, K b,0 , N 0 , and ρ s are shown in Fig. 4.12. As was the case for the estimates of the numerical study, tortuosity, porosity, and solid density are well identified in the sense that their marginal posterior distributions are narrow compared to the prior distributions. Permeability also exhibits a trend similar to that in publication II, where it has a definite lower bound, but can admit values up to the highest allowed value of 10 -8 m 2 . However, as opposed to the finding of publication II, both K b,0 and N 0 seem to be well identified even from the normal incidence data. One clear reason for this is the condition 2/3N 0 ≤ K b,0 that was added to the prior. Indeed, the estimated values of K b,0 ≈ 13 GPa and N 0 ≈ 18 GPa are close to this relation, and the estimated Poisson's ratio of the porous frame is very close to zero. In contrast, at this location the estimate for the bulk modulus of the solid is K s,0 ≈ 22 GPa, and the Poisson's ratio is about 0.18.

The effect of the prior on K b,0 and N 0 can also be seen in the joint marginal densities, shown in Fig. 4.13, where the line 2/3N 0 = K b,0 is clearly visible. Other correlations can be seen between the solid density and the bulk and elastic moduli, as well as between tortuosity and the solid bulk modulus. However, the observed correlations are lower than in the numerical case. This may be explained by the fact that with the real measurements we have a lot more model error, and therefore the model together with the data behave in a different way.

Let us now summarise the results in Fig. 4.14, which shows the CM and 95 % CI estimates, for selected parameters, of all 216 measurements. The estimates are plotted according to where the measurement data they are based on was taken at, and hence the images show how the parameter estimates and their uncertainties vary spatially. We can see immediately that overall the parameters show variability over the scale of the plate, but the small-scale variability from one measurement location to the next is relatively small. Only at a very small point there is a sudden jump in the estimated values (visible especially in the estimates of Λ and α ∞ ).

Because the inversion at each point is carried out independently of the other locations, the smooth varying of the parameters must be found in the measurement data. The smoothness in the data may be explained by two factors. The first possibility is that the resolution of the measurements is significantly lower than the spacing of the measured locations, i.e. the effective area seen by the transducers is much larger than the 10 mm pulse/echo beam diameter. If this were the case, however, we would see a low resolution in all parameters whereas now the estimates of φ and ρ s show structures much smaller than the ones seen in α ∞ and K b,0 . We can then conclude that the observed smooth variability in the parameters is a property of the object.

Finally, Fig. 4.15 shows the spatial variability of the estimated nuisance parameters. None of the nuisance parameters can be seen to vary in the same way as the Biot parameters, which further shows that the variability in the pore parameters comes from the measured plate. Another notable observation is that the nuisance parameters can be identified accurately, as shown by the relatively narrow widths of the 95 % credible intervals. For example, the distance from the pulse/echo transducer to the plate is estimated to vary by up to 3 mm (the distance between the transducer and the water-air interface is constant in the reference measurement, so changes in ε R show directly the changes in distance between the transducer and the plate), but the associated uncertainty interval for the distance estimate is less than 0.06 mm. The estimates of ε R also show that the plate is not completely planar. In addition, the thickness of the plate can be seen to vary by almost a millimetre, with an uncertainty of less than 0.15 mm. These nuisance parameters are clearly identified independently, and accounting for them during the inversion thus improves the estimates of the Biot parameters. The CM and 95 % CI estimates of selected parameters at every measured location (the pixelated images), and a schematic that shows the plate and the locations inside it that were measured (denoted by the crosses). The parameter estimates were calculated separately for each measurement. Therefore, the colour of the pixels denotes the value of the of the estimated quantity, and the location of the pixels corresponds to the location where the measurement was made.

Summary and conclusions

In this thesis, we developed Bayesian characterisation methods for poroviscoelastic and rigid frame porous media. The motivation was to improve the reliability and repeatability of the characterisation, by considering measurement and model errors, incorporating prior information, and finally approximating the true posterior probability density of the parameters using sampling methods. In particular, our goal was to provide feasible error estimates for each unknown parameter. We considered materials that were saturated either with air or water, and tested the proposed methods with real measurement data.

In publication I, we proposed a Bayesian method for the characterisation of rigid frame materials, and compared the results to a least squares (LS) estimate, and other characterisation methods. We used a standard transmission impedance tube to measure the materials' equivalent density and bulk modulus, and the Johnson-Champoux-Allard-Lafarge (JCAL) model [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF][START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF][START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF] to connect the measurement data to six physical parameters. The parameter inference in the Bayesian framework was then carried out using an adaptive Markov chain Monte Carlo (MCMC) technique, and a point estimate was reported in the form of the maximum a posteriori (MAP). The results showed that the measurement data carried information on all JCAL parameters, and that the MAP estimates were within a few percent of those found using the reference methods. In addition, the LS estimates were found to be nearly identical to the MAP estimates, but only if the minimisation of the cost functional was done using a multistart approach. Otherwise the LS solution often returned a local minimum. The statistical approach also showed that most parameters are identified independently, but with sometimes certain parameters may exhibit strong correlation with one another. However, the correlations were found to be different for each material, which is why it is informative to summarise the posterior using MCMC each time a new material is measured.

Publication II was a numerical study in which we developed inversion methods for poroviscoelastic media, and designed an experimental setup where the characterisation could be also done in practice. To make the inversion using sampling methods possible, we needed to have a simple measurement geometry which could be modelled with fast numerical methods. For this reason, the measurement setup consisted of a flat slab of porous material submerged in a water tank. Such a setting can be considered as a layered system, and if we assume plane wave propagation, several fast methods are available to solve the equations. The porous medium was modelled using the Biot equations, with the Johnson et al. model [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF] to account for viscous losses. The forward model was then constructed using a global matrix method, which is both a fast and stable way to solve wave propagation problems in layered media. We simulated the measurement data, which consisted of the acoustic reflection and transmission coefficients of the porous plate, with a finite-element model, and considered several angles of incidence.

A key contribution was to construct an adaptive parallel tempering MCMC algorithm, and to demonstrate that it is a viable method to compute the parameter estimates in the considered setting. The results showed that many interesting pa-rameters, including porosity, can be estimated from just the measurements at normal incidence, and incorporating measurements from oblique angles would not improve the estimates much. These parameters include porosity and permeability. However, if accurate information on the elastic parameters is desired, measurements at oblique angles of incidence are recommended so that the effect of the shear waves can be seen.

In publication III, we implemented the experimental design from publication II in a water tank measurement system in the Laboratoire d'Acoustique de l'Université du Mans (LAUM), carried out extensive measurements of a porous ceramic plate, and applied the inversion method introduced in II to the measurements. The goal was to evaluate the method with real measurements and to estimate the possible spatial changes in the parameters. The actual measurement system consisted of focused ultrasound transducers connected to a computer-controlled positioning system, which was used to measure the reflection and transmission coefficients at several locations. However, due to the limitations of the system, we were able to make the measurements only at normal incidence. The inversion was carried out independently at each measurement location. In addition to the physical parameters of the plate, we estimated some nuisance parameters related to the measurement setup, such as the measurement noise level.

The measured reflection and transmission coefficients were found to vary significantly from location to location, suggesting that the plate is inhomogeneous. Nevertheless, the fit of the Biot model to the measurement data was equally good for nearly every measurement, and the measurements were contained within the envelope of the posterior predictive distribution, which indicates that the employed model is a feasible one. However, we found that the residuals between the model and the measurements were not distributed exactly as white noise. This implies that there are features of the measurement system and the measured object that are not accounted for in our forward model, such as possible depth-wise inhomogeneities of the plate, errors resulting from the plane wave approximation, and attenuation due to scattering effects. Modelling errors can obviously introduce bias into the parameter estimates, but quantifying this was not pursued in this thesis. We also found that the observed model residuals were significantly smaller than the variability in the measurements between different locations of the ceramic plate.

The inversion results for each unknown were plotted as two-dimensional graphs which showed the conditional mean (CM) estimate and the length of the associated 95 % credible interval at each measurement location. The parameter estimates were different depending on the location of the measurement, but at the same time the changes in the values were gradual and smooth. In addition, the parameters did not all vary in the same way, and in particular the pore structure parameters did not vary in the same way as the nuisance parameters, such as the plate thickness, which suggests that the observed variation is real and not an artefact of the model errors. In conclusion, the results indicate that the employed computational model can be used to feasibly estimate the parameters of viscoelastic porous media.

Although the results of this thesis have shown that the proposed methods are a promising step towards the robust characterisation of porous materials, a lot of future work remains. For instance, throughout the thesis we assumed that we have only plane wave propagation, which can be easily justified for an impedance tube, but is only an approximation in free-field measurements. Modelling the finite, nonplanar, field of the ultrasound transducers would be an obvious improvement for the inversion, and would let us quantify possible biases the plane-wave approxima-tion leads to in the parameter estimates. One way to model the true ultrasound field is to use element-based methods such as the finite element or the boundary element methods. This would also enable the characterisation of objects that are of irregular shape. However, modelling complex objects and acoustic fields using elementbased approaches will substantially increase the computational burden, which may quickly render the proposed MCMC method infeasible. There is clearly a need for fast computational methods that can still account for a finite-sized acoustic field.

One possible approach to model more realistically the fields emitted by ultrasound transducers, without having to resort to full-wave solutions, is to use the superposition of harmonic plane waves. The idea is that the finite incident field is decomposed into infinite plane wave components using the Fourier transform, and then a plane wave solution in the layered media is computed for each component. Finally, to find the reflection and transmission coefficients, the responses of each plane wave component are summed together [START_REF] Pialucha | The reflection coefficient from interface layers in NDT of adhesive joints[END_REF][START_REF] Vecchio | Prediction of ultrasonic field propagation through layered media using the extended angular spectrum method[END_REF]. The computational requirements of the plane wave decomposition approach is naturally directly proportional to the number of plane waves considered.

Another possibility in the Bayesian framework is to construct a statistical model for the errors that are induced when we use approximate forward models, and to marginalise over them before the inversion. This is called the Bayesian approximation error (BAE) approach [START_REF] Kaipio | Approximate marginalization over modeling errors and uncertainties in inverse problems[END_REF][START_REF] Kaipio | Statistical inverse problems: discretization, model reduction and inverse crimes[END_REF]. An advantage of using the BAE is that once the statistical model for the approximation errors (the discrepancies of the predictions between an approximative and a more accurate model) have been computed, we can use the fast approximative model in the inversion, while partly recovering from the errors. Applied to characterisation of porous media, we could consider a finite element model as the accurate one, and marginalise over the approximation errors when using a plane wave model.

Solving the inverse problem using MCMC methods is, however, always going to be relatively slow. If a significantly faster, or even real-time, characterisation method is needed, we have to compute the MAP estimate which can often be obtained with only a few iterations, and use the Laplace approximation to get an estimate of the posterior covariance.

During the thesis work, we implemented and tested numerous MCMC algorithms for the problem of poroelastic material characterisation, and found that adaptive Metropolis and parallel tempering worked efficiently. However, there are certainly other good algorithms that could work equally well or better, and new MCMC samplers are being developed all the time. While parallel tempering is one of the best ways to deal with multimodal densities, other methods such as Hamiltonian Monte Carlo [START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF] may be more efficient in different settings, especially if derivatives of the posterior are easily obtained. Additionally, nowadays there are many freely available MCMC software packages which can facilitate easy out-of-the-box sampling [START_REF] Foreman-Mackey | emcee: The MCMC Hammer[END_REF][START_REF] Carpenter | Stan: A Probabilistic Programming Language[END_REF][START_REF] Lunn | WinBUGS -A Bayesian modelling framework: Concepts, structure, and extensibility[END_REF][START_REF] Salvatier | Probabilistic programming in Python using PyMC3[END_REF].

The measurements of the porous ceramic layer carried out during this thesis were done only at normal incidence. This was mostly a limitation of the measurement system in the water tank, which did not allow for accurate measurements of the transducer angle that would be needed for the oblique incidence measurements. Recently, LAUM has acquired a new computer-controlled measurement system which would allow these measurements. An obvious extension would therefore be to test the proposed inversion method with measurements at oblique angles, which would give more information especially on the bulk and shear moduli that were found to be correlated at normal incidence.

The inversion method should also be tested with control materials whose properties are known. The problem of obtaining such porous materials is that porous materials are inherently heterogeneous, and reliable control materials manufactured using traditional methods are therefore hard to come by. On the other hand, because the pores need to be very small for the setup and frequencies we have considered (on the order of micrometers), 3D printing such structures is only now starting to become possible.

Another generalisation for the characterisation of porous materials would be to account for the possible macroscopic heterogeneity and anisotropy of the samples. Macroscopically heterogeneous porous materials offer, for example, more possibilities for the design of sound absorbing materials [START_REF] Geslain | An application of the Peano series expansion to predict sound propagation in materials with continuous pore stratification[END_REF][START_REF] Mahasaranon | The effect of continuous pore stratification on the acoustic absorption in open cell foams[END_REF], and anisotropic materials are found in many porous structures such as fibrous sound and heat insulation materials and sediment structures, as well as acoustic metamaterials [START_REF] Christensen | Anisotropic Metamaterials for Full Control of Acoustic Waves[END_REF]. Gradually changing material profiles in one dimension have been estimated in the rigid frame case using the wave splitting and invariant imbedding techniques, with a LS minimisation approach [START_REF] Ryck | Reconstruction of material properties profiles in one-dimensional macroscopically inhomogeneous rigid frame porous media in the frequency domain[END_REF]. Another possibility to model macroscopically inhomogeneous layered materials is to use the reversible jump MCMC methods [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF], where the number of the layers can be another unknown. Such methods provide estimates with uncertainties for the parameters within the layers, but also include the uncertainty in the model parametrisation [START_REF] Dosso | Efficient transdimensional Bayesian inversion for geoacoustic profile estimation[END_REF][START_REF] Gehrmann | Trans-dimensional Bayesian inversion of controlled-source electromagnetic data in the German North Sea[END_REF]. An approach for the characterisation of three-dimensional anisotropic materials has been proposed in [START_REF] Terroir | General method to retrieve all effective acoustic properties of fullyanisotropic fluid materials in three dimensional space[END_REF].

Finally

  , I would like to thank my partner Lily. You have travelled across the vii world with me and have shown endless love and support during these years, for which I cannot thank you enough. Kuopio, January 17th, 2020 Matti Niskanen viii

Figure 2 . 2 :

 22 Figure 2.2: Left: Illustration of a tortuous path. Right: The open pore space regions

Figure 2 . 4 :

 24 Figure 2.4: Representation of a multilayered system. Different media are denoted

  k |y)dθ k = (100p)/200. (3.7) Another common way to define the interval is to minimise the distance between the endpoints [a I , b I ] = arg min a I ,b I b Ia I b I a I π(θ k |y)dθ k = p/100 . (3.8)

Figure 3 . 1 :

 31 Figure 3.1: The first 45,000 iterations of one parameter of an example MCMC run.

  First let us comment on the RWM sampler, which most of the more advanced samplers are based on. The proposal density in RWM is usually Gaussian, ∼ N (0, Σ), or uniform,

Figure 3 . 2 :

 32 Figure 3.2: A mixture of two Gaussian distributions with variance 1 and means ±10, normalised after tempering. At T = 100 the peaks are widened so much that we have a unimodal distribution.

Figure 4 . 1 :

 41 Figure 4.1: Samples measured in publication I: (1) soft polyurethane foam, (2) acoustic foam, (3) melamine foam, (4) felt, (5) glass wool.

Figure 4 . 2 :

 42 Figure 4.2: The transmission impedance tube used in the measurements.

Figure 4 . 3 :

 43 Figure 4.3: Measurements of density and bulk modulus, shown as real and imaginary parts, of two of the samples (black continuous lines), and the model prediction corresponding to the MAP estimate plotted on top (dotted lines).

Figure 4 .

 4 [START_REF] Liu | Porous Materials[END_REF] shows posterior summaries for the acoustic foam sample. Each subplot shows the joint marginal distributions of the estimated parameters, the one

Figure 4 . 4 :

 44 Figure 4.4: Marginal posterior densities of the acoustic foam sample. The red curves on the bottom and left side of each plot are the 1D marginal densities, and the contour lines indicate the 50 and 95 % regions of distribution. The dots represent a sample of the MCMC points, and r is the Pearson correlation coefficient. The black cross denotes the MAP-estimate, and the red circle the LS-solution.

Figure 4 . 5 :

 45 Figure 4.5: A schematic of the computational model.

Figure 4 . 6 :

 46 Figure 4.6: An example of the simulated noisy measurement (at normal incidence) and the fit of the plane wave Biot model.

Figure 4 . 7 :

 47 Figure 4.7: The marginal prior and posterior distributions of tortuosity, permeability, porosity, and the real part of the bulk modulus, estimated from measurements at normal incidence. Also shown is are the true values, the CM estimates, and the 95 % credible intervals.

Figure 4 . 8 :

 48 Figure 4.8: The joint marginal distributions and the true values of selected parameters, estimated from measurements at normal incidence.

  Figure 4.8 shows the joint distributions of pa-

Figure 4 . 9 :

 49 Figure 4.9: Parameter estimates and 95 % credible intervals as a function of the measurement angle ϕ inc .

Figure 4 . 10 :

 410 Figure 4.10: A schematic view of the measurement system and the poroelastic plate, a) reflection measurement, b) transmission measurement. The shaded zones indicate the -6 dB size of a pulse/echo signal on the surface of the plate.

Figure 4 . 11 :

 411 Figure 4.11: An example of typical measured reflection and transmission coefficients, and the fit of the Biot model. The shaded areas denote the one, two, and three standard deviations of the posterior predictive distribution.

Figure 4 . 12 :

 412 Figure 4.12: The marginal prior and posterior distributions of selected parameters. Estimated from an example measurement in the water tank. Also shown is are the CM estimates and the 95 % credible intervals.

Figure 4 . 13 :

 413 Figure 4.13: The joint posterior distributions and the CM estimate of selected parameters, estimated from the water tank measurements.

Figure 4 .

 4 Figure 4.14: The CM and 95 % CI estimates of selected parameters at every mea-

Figure 4 . 15 :

 415 Figure 4.15: The CM and 95 % CI estimates of the nuisance parameters.

TABLE OF CONTENTS

 OF Rigid frame porous media................................................................. 2.1.1 The Johnson-Champoux-Allard-Lafarge model ........................ 2.2 Poroelastic media............................................................................. 2.2.1 The Biot model ..................................................................... 2.2.2 Losses in water-saturated materials ........................................ 2.3 Models to compute the acoustic response of a multilayer system ....... Computational methods for estimating the material parameters of poroelastic media .......................................................

	1 Introduction
	2 Wave propagation in porous media
	2.1

2.3.1 Waves in layered media.......................................................... 2.3.2 The transfer matrix method ................................................... 2.3.3 The global matrix method ..................................................... 3 Bayesian methods for parameter estimation 3.1 Parameter estimation in the Bayesian framework .............................. 3.1.1 Likelihood ............................................................................. 3.1.2 Prior models.......................................................................... 3.2 Markov chain Monte Carlo ............................................................... 3.2.1 Metropolis-Hastings update ................................................... 3.2.2 Efficiency measures................................................................ 3.2.3 Stopping criteria and diagnosing convergence ......................... 3.3 MCMC Samplers.............................................................................. 3.3.1 Overview ............................................................................... 3.3.2 Adaptive algorithms............................................................... 3.3.3 Parallel tempering.................................................................. 4 Review of thesis results 4.1 Publication I: Characterisation of rigid frame porous materials from impedance tube measurements ......................................................... 4.1.1 Measurement setup ............................................................... 4.1.2 Inversion ............................................................................... 4.1.3 Results .................................................................................. 4.2 Publication II: 4.2.1 The simulation model ............................................................ 4.2.2 Inversion ............................................................................... 4.2.3 Results .................................................................................. 4.3 Publication III: Estimating the material parameters of poroelastic objects from ultrasound measurements in water ................................ 4.3.1 The measurement setup......................................................... 4.3.2 Inversion ............................................................................... 4.3.3 Results .................................................................................. xi 1 Introduction