Directeur De Thèse Professeur

Christian Soize

Professeur, rapporteur Jury M Keith Attenborough

Directeur M Philippe Blanc-Benon

Professeur, rapporteur M Michel Bruneau

M François

M Jérôme Defrance

M Philippe Docteur

Docteur Jean

M Christian Hdr

Professeur Soize

Keywords: energy dissipation, 2, 3, 2 High-angle formulation

Dans cette recherche, un modèle de propagation d'ondes de choc sur grandes distances sur un environnement urbain est construit et validé. L'approche consiste à utiliser l'Equation Parabolique Nonlinéaire (NPE) comme base. Ce modèle est ensuite étendu afin de prendre en compte d'autres effets relatifs à la propagation du son en milieu extérieur (surfaces non planes, couches poreuses, etc.). La NPE est résolue en utilisant la méthode des différences finies et donne des résultats en accord avec d'autres méthodes numériques. Ce modèle déterministe est ensuite utilisé comme base pour la construction d'un modèle stochastique de propagation sur environnements urbains. La Théorie de l'Information et le Principe du Maximum d'Entropie permettent la construction d'un modèle probabiliste d'incertitudes intégrant la variabilité du système dans la NPE. Des résultats de référence sont obtenus grâce à une méthode exacte et permettent ainsi de valider les développements théoriques et l'approche utilisée.

A.1 Overpressure values and corresponding damage to structures from [START_REF] Lannoy | Analyse des explosions air-hydrocarbure en milieu libre : Etudes déterministes et probabilistes du scénario d'accident[END_REF] and [START_REF] Lees | Loss Prevention in the Process Industries[END_REF] Jérôme, merci pour ton énergie et les très nombreuses mais néanmoins intéressantes rencontres que j'ai pu faire grâce à toi. Philippe, merci pour ton aide précieuse et les découvertes musicales. Je pense que je dois pouvoir effectuer une deuxième thèse sur les groupes rock des années 70/80 ! Christian, merci pour ta rigueur scientifique, ô combien précieuse ! Chaque discussion fût une voie de sortie dans le tunnel, une remise sur les rails. Une thèse de doctorat est avant tout une formation au métier de chercheur. Je pense avoir eu au travers de vous trois un aperçu complet des multiples facettes de ce métier, et je vous en remercie chaudement.

Et enfin, last but not least, je tiens à exprimer ma profonde gratitude envers ceux qui, souvent de manière indirecte et involontaire, ont permis que ce travail naisse, vive, et aboutisse.

Effectuer une thèse de doctorat est une tâche longue, parfois éreintante et décourageante, jamais simple, toujours passionnante. Ce travail ne serait pas ce qu'il est si mes proches, famille, petite amie et amis, n'avaient pas été là pour me changer les idées, me ramener sur terre et me soutenir de multiples façons. Ce document est également un peu le vôtre. Merci. This research was performed within the MSME laboratory at Université Paris-Est (UPE) Marne-la-Vallée in strong partnership with the Acoustics and Lighting Department (environmental and urban acoustics group) of CSTB. What gets us into trouble is not what we don't know. It's what we know for sure that just ain't so.

Grenoble

Mark Twain (1835 -1910)

General introduction and context of research

The motivation for this work, entitled Nonlinear acoustic wave propagation in complex media -Application to propagation over urban environments, finds its origin in events such as the explosion of the AZF factory in Toulouse (France) in 2001. This accident was responsible for the death of 30 people, thousands of injured and cost around e3 billion.

Today there are still some facts about this event, such as a double echo, that scientists are unable to explain with the tools they possess. The work presented in this document tries to provide an answer to such problems, through the development of a numerical method for the study of blast wave propagation in complex media and over large distances, with the undercurrent objective of risk assessment of such industrial sites.

Industrial explosion risk and consequences

Industrial explosions since 2001

Since 2001, eight major industrial explosions have been reported. These accidents caused more than 80 casualties, tenth of thousands of injured and cost billions of euros. Tab. (1) lists these events and gives their description. With eight major accidents in less than a decade (approximately 1 accident per year), the risk and the consequences of industrial explosions are far from being nonexistent.

Two major events can be denoted from Tab. (1): the explosion of the AZote Fertilisant (AZF) factory in Toulouse (France) in 2001 and the explosion of an oil depot in Buncefield (England) in 2005. These two accidents killed 30 people, injured thousands of others and cost around e3 billion. ◮ Thermal effects from the combustion of inflammable materials. Depending on the material and the detonation conditions, the combustion time can vary from some milliseconds to several hours. Human consequences are burns; it can inflame neighboring structures.

◮ Overpressure effects which result from the propagation of a high pressure wave caused by an explosion. Effects on humans range from eardrum rupture to death by lung hemorrhage. Buildings can be heavily damaged and even collapse.

◮ Toxic effects resulting from a leak of a toxic substance. The most exposed are employees and emergency services.

In this work we only consider the effects from the overpressure wave. Possible damage to buildings and injuries caused by an overpressure wave are briefly presented below.

Damage to constructions caused by blast waves

Common units to characterize an overpressure amplitude are millibars and psi. In this document the overpressure amplitude is given in kPa, so that it can easily be linked to acoustical quantities. As a reference for comparison, a charge of 1 kg of Trinitrotoluene (TNT) produces a positive peak overpressure of approximately 10 kPa at 10 meters from the explosion.

Tab.

(2) lists the overpressure levels and the corresponding observed damage on structures.

The figures are taken from [START_REF] Merrifield | Simplified calculations of blast induced injuries and damage[END_REF] (figures from different sources can be found in Appendix A). For very low overpressures, typically less than 1 kPa, very few damage to buildings will occur. From 1 kPa to 5 kPa, minor damage will be observed: cracks in walls, windows shattered, tiles displaced, etc. These overpressure strengths are not directly harmful for humans (see next section), however, indirect harms through broken glass or projected objects can be important. From 5 kPa to 15 kPa, serious damage will be observed on buildings and the environment: inhabitable houses, partial collapse of walls and roofs, trees blown down, etc.

Injuries caused by blast waves

Direct injuries to people are presented in Tab. (3); this table does not include indirect harms caused, for example, by fallen tiles or broken glass. The threshold of fatality for a vulnerable population is very low: 7 kPa is less than what is produced by 1 kg of TNT at 10 meters.

Most common direct injuries are eardrum ruptures: more than 500 cases were reported after the AZF explosion. The threshold of fatality for a healthy population is around 15 kPa.

Psychological consequences, that do not appear in Tab.

(3), are also important. The weeks following the accident in Toulouse, more than 700 psychiatric consultations were recorded, 3 200 psychoactive drug prescriptions were given and more than 5 000 people were reported to suffer from acute stress [see [START_REF] Sauret | Le traumatisme psychologique secondaire à l'explosion de l'usine AZF à Toulouse (Psychological trauma subordinate to the explosion of the AZF factory)[END_REF]. During the five days following the AZF explosion, the number of heart attacks was three times more than the year before [see [START_REF] Rivière | Conséquences sanitaires de l'explosion survenue à l'usine AZF, le 21 septembre 2001 -Rapport final sur les conséquences sanitaires dans la population toulousaine[END_REF].

Two examples: the AZF factory and the Buncefield oil depot explosions Two of the most important industrial explosions occurred during the last decade in Europe.

Both were impressive by their amplitude and illustrate the consequences on humans and their Table 3: Commonly accepted figures for direct harm to people from blast overpressure. This table is taken from [START_REF] Merrifield | Simplified calculations of blast induced injuries and damage[END_REF].

Overpressure (kPa) Direct harm to people 7 Threshold of fatality (1-5 %) for a vulnerable population The crater is still clearly visible at the top of the picture. The explosion was heard 80 km away and was reported as occurring in multiple places; police at first believed that at least five bombs had simultaneously gone off. There is still controversy over the exact number of explosions. The results of the official inquiry were that a warehouse of ammonium nitrate had exploded following improper handling of this dangerous material, including mixing with chemical impurities. Specifically, it is believed by mandated official experts that * "a mislabeled 500 kg bin of sodium dichloroisocyanate mistakenly thought to be ammonium nitrate was dumped in the off spec ammonium nitrate warehouse.

Here under sufficiently hot and humid conditions it could have reacted with the ammonium nitrate to form nitrogen trichloride which is an exceedingly unstable compound. The decomposition of the nitrogen trichloride could have provided the heat and pressure required to detonate the ammonium nitrate which when used as an industrial explosive typically requires detonators."

Human, environmental and economic impact

The factory was close to the city: one of the most inhabited areas, Le Mirail, is just one kilometer away from the factory site. Several schools, one university campus, one hospital and a psychiatric hospital had to be evacuated. The disaster caused 30 deaths, 2 500 seriously wounded and 8 000 light casualties. Two thirds of the city's windows were shattered, causing Buncefield oil depot, Buncefield (England), 2005

The accident

The first and largest explosion occurred at 06:01 UTC near container 912. From all accounts, it seems to have been an unconfined vapor cloud explosion of unusually high strength. The explosions were heard 150 km away; there were reports it was audible in France and the Netherlands. The British Geological Survey monitored the event, which measured 2.4 on the Richter scale. It was reported that people were woken in south London, and as far west as Wokingham, where numerous people felt the shock wave after the initial explosion.

Subsequent explosions occurred at 06:27 and 06:28. Witnesses observed flames hundreds of feet high from many miles away, with the smoke cloud visible from space.

Human, environmental and economic impact

Damage from the blasts, ranging from broken windows and blown-in or warped front doors to an entire wall being removed from a warehouse, occurred more than 800 m away. Buildings in neighboring St Albans also suffered. For example, Townsend School had serious blast damage, and a window was blown out of St Albans Abbey (both 8 km from the site). Several nearby office blocks were hit so badly that almost every window, front and back, was blown in as the explosion ripped through them. Reports also indicated that cars in nearby streets caught fire. The roof of at least one house was blown off. Fig.

(2) illustrates the explosion impact on buildings and the environment: more than 400 meters from the explosions, buildings are seriously damaged. There were 43 reported injuries; two of them being seriously harmed.

Since the explosion occurred early in the morning, very few people were on site, which explain the low number of casualties. This explosion cost almost a billion pounds (see Tab. ( 4)), the main part being compensation claims.

Conclusions

Eight major explosions have occurred during the last decade. About a hundred people were killed, thousands of others were harmed and the overall cost of these events exceeded e3 billion. This work aims at developing a numerical simulation tool for the study of propagation of blast waves from industrial explosions. This would allow the physics of the phenomenon to be better understood and ultimately, this tool could be used to assess the risks inherent to a particular installation under certain circumstances (explosion strength, meteorological conditions, etc).

The physical problem of blast wave propagation falls into the general category of acoustics and more specifically of nonlinear outdoor sound propagation. Chapter 1 introduces the main features of high-amplitude wave propagation outdoors and details the work presented in this document.

1

Main features of nonlinear sound propagation outdoors

This chapter presents the features of high-amplitude sound propagation outdoors that are of interest for the study of blast wave propagation in complex media. Section 1.1 details the features of sound propagation outdoors that have a significant effect on sound pressure levels (geometrical spreading, atmospheric absorption, meteorological effects, ground effects) while Section 1.2 briefly summarizes the high-amplitude effects and their consequences on wave propagation. Section 1.3 reviews the studies specifically dedicated to blast wave propagation and finally, conclusions are given and the content of this document is described in Section 1.4.

Chapter content

Features of sound propagation outdoors

This section is a review of most areas of sound propagation outdoors that are of interest for long-range acoustic applications. Because the field is so diffuse only the dominant mechanisms are presented here. More detailed information can be found in various papers [see e.g. [START_REF] Embleton | Tutorial on sound propagation outdoors[END_REF][START_REF] Piercy | Review of noise propagation in the atmosphere[END_REF][START_REF] Ingård | A review of the influence of meteorological effects conditions on sound propagation[END_REF] or in more complete publications like the books by [START_REF] Salomons | Computational atmospheric acoustics[END_REF] and [START_REF] Attenborough | Predicting outdoor sound[END_REF].

The first two areas covered in this chapter are geometrical spreading and atmospheric absorption (Sections 1.1.1 and 1.1.2, respectively). Together these are the dominant mechanisms determining the sound levels outdoors. Ground effects are treated in Section 1.1.3 and effects of meteorological conditions on sound propagation in Section 1.1.4.

Geometrical spreading

Waves spread in three dimensions when the sound source is small compared to the distances being considered. The resulting attenuation depends on the propagation distance and is frequency independent. For a spherical sound source in a homogeneous medium in free field the acoustic power is uniformly spread on a spherical wave front. The wave front area being proportional to the square of the sphere radius, the acoustic intensity decays by 1/r 2 and the acoustic pressure by 1/r, r being the distance from the source to the receiver. At twice the distance from the source, the wave front area is four times as large, and the Sound Pressure Level (SPL) decreases by about 6 dB. Each time the distance is doubled, the SPL decreases by the same quantity.

Sound waves spread cylindrically from a line of sources which are all similar but radiate independently. The area of the cylindrical wave front is proportional to the distance, the SPL thus decreases by about 3 dB per doubling of distance, at half the rate of spherical spreading.

Atmospheric absorption

Atmospheric absorption is a feature of wave propagation that is always present in outdoor sound propagation. Although it may be neglected in some applications, this phenomena has to be carefully taken care of in long-range applications. The global absorption phenomenon is mainly due to three physical effects:

◮ classical absorption caused by the transport processes of classical physics (shear viscosity, bulk viscosity, thermal conductivity, mass diffusion and thermal diffusion),

◮ rotational relaxation of molecules in air,

Features due to the ground

In this section the different roles of the ground in sound propagation are detailed. First, the case of sound reflections on a flat and rigid ground and its resulting interference pattern are detailed. The second section deals with the case of non-rigid ground surfaces; different impedance models are briefly presented.

Path-length differences (PLD)

When both the source and the receiver are above the ground a phase change occurs due to the different lengths of direct and reflected waves paths. This phase change occurs in addition to the change due to the ground characteristics.

Ground impedance

The reflection coefficient R p for a plane wave on a locally reacting plane surface is

R p = sin(ψ) -Z 0 /Z sin(ψ) + Z 0 /Z , (1.1)
where ψ is the grazing angle, Z 0 = ρ 0 c 0 is the characteristic impedance of air and Z is the acoustic impedance of the surface. In the frame of a locally reacting surface, an hypothesis which has been proved to be valid for most of the surfaces encountered in outdoor sound propagation applications [see for example the papers and the book [START_REF] Attenborough | Acoustical impedance models for outdoor ground surfaces[END_REF][START_REF] Attenborough | Ground parameter information for propagation modelling[END_REF][START_REF] Salomons | Computational atmospheric acoustics[END_REF], the acoustic characteristics of the ground may be represented by its acoustical impedance Z = R + jX. This impedance value may depend on the frequency but not [START_REF] Zwikker | Sound absorbing materials[END_REF] impedance model. The ZK model originates from linearized acoustic equations for a plane wave traveling in a rigidly-framed porous material. In one dimension they can be generalized as follows:

-

∂ x v = Ω 0 ρ 0 c 2 0 ∂ t p ′ , (1.2a) Φ Ω 0 ρ 0 ∂ t v + σ 0 v = -∂ x p ′ . (1.2b)
In the above Eqs. (1.2), ρ 0 is the ambient medium density, c 0 is the ambient sound speed, v is the particle velocity, p ′ is the acoustic pressure disturbance, Φ the material tortuosity, Ω 0 the ground porosity, σ 0 the flow resistivity of the material. The time variable is t while the space variable is x, and ∂ i means partial derivation with respect to variable i.

The tortuosity is defined as the ratio of the arc length of the pore structure to the distance between its end points. The porosity is the ratio of the air volume contained in the porous material to its total volume.

The derivation presented in the paper by [START_REF] Salomons | Computational atmospheric acoustics[END_REF] yields an expression for the characteristic impedance Z c of the porous material, defined by

Z c = ρ 0 c 0 Φ Ω 2 0 + i σ 0 Ω 0 ρ 0 ω , (1.3)
in which c 0 is the ambient sound speed. This model has been shown to give good agreement with more complex models over a wide range of frequencies [see e.g. [START_REF] Wilson | Time-domain modeling of the acoustic impedance of porous surfaces[END_REF][START_REF] Wilson | Time-domain calculations of sound interactions with outdoor ground surfaces[END_REF][START_REF] Wilson | Time-domain equations for sound propagation in rigid-frame porous media[END_REF]. [START_REF] Delany | Acoustical properties of fibrous absorbent materials[END_REF] impedance model. It is an empirical model that allows absorbing fibrous materials to be modeled. For a material with flow resistivity σ 0 (here in N.m -4 .s), the expression of the characteristic impedance Z c is

Delany

Z c = ρ 0 c 0 1 + 9.8 10 3 f σ 0 -0.75 + j11.9 10 3 f σ 0 -0.73 , (1.4)
where f is the wave frequency. Although this model is simple (the impedance only depends on the frequency and the flow resistivity) it has been shown that the results obtained are in good agreement with measurements. The DB impedance model is commonly used in outdoor sound propagation applications.

The validity of both models mainly depends on the frequency and the material considered.

Thanks to the use of four ground parameters, the ZK model can generally be used in a wider range of applications than the DB impedance model, which uses a very simple parametrization.

Note that for low frequencies the ground impedance absolute value becomes very high and thus the reflection coefficient is close to unity (almost acoustically rigid ground).

Effects of meteorological conditions

In this section meteorological conditions usually encountered in outdoor sound propagation and their effects on wave propagation are described.

Relevant meteorological phenomena

The complete meteorological description of an environment is a complex task. A lot of quantities such as temperature, wind velocity, air density, air pressure, intervene in this description. These quantities are obviously inter-independent and most of the time vary in space and time. For a more exhaustive description of the different meteorological phenomena we refer the reader to specialized literature such as the book by [START_REF] Munn | Descriptive micrometeorology[END_REF].

Nevertheless, from an acoustical point of view one can denote two dominant parameters that will modify sound propagation paths: wind speed gradients and temperature gradients.

These two quantities locally affect the effective sound speed and thus change the way sound propagates.

Temperature gradients. During daytime the variation of temperature with height for a large flat area may be represented by the expression Wind velocity gradients. Wind velocity is a three-dimensional vector quantity difficult to represent due to its spatial and temporal instability. The variation of the average wind speed V with height z in the vicinity of the ground for a flat area is approximately as shown in Fig. (1.5). This wind speed profile may be represented with a logarithmic law; for altitudes greater than z 0 the average wind speed can be expressed by

T = T 0 -K t log (z/z 0 ) , (
V = K v log (z/z 0 ) . (1.6)
The parameter z 0 is determined by the roughness of the surface, often approximately a characteristic height of obstacles. The constant K v is related to the surface roughness and the average wind velocity above this layer. If we neglect the vertical component of the average wind velocity (which is often very small in comparison to the average values for the horizontal component), we can write the average sound speed with the following formula:

c(x, z) = c 0 + V x , (1.7)
in which V x is the average wind speed in the horizontal direction and c (x, z) is the spatiallydependent average sound speed. [START_REF] Piercy | Review of noise propagation in the atmosphere[END_REF].

Combined effects of wind velocity gradients and temperature gradients

Combined effects of temperature gradients and wind velocity gradients can yield very complex situations where sound ray paths can adopt various trajectories. Below are presented two cases taken from the article by [START_REF] Ingård | A review of the influence of meteorological effects conditions on sound propagation[END_REF] which illustrates complex situations that can be found in outdoor sound propagation applications.

Case 1 (see Fig.

(1.8(a))): there is a change of wind gradient at the layer 3 km above the ground. Above this layer the wind gradient can no longer compensate for the temperature gradient, and the rays will all be bent upwards in this region. The rays in the vertical plane in the wind direction are shown to reach out to a distance of 38.8 km from the source. Beyond this distance there will be a shadow zone.

Case 2 (see Fig.

Sound propagation in urban environments

With the increasing need of transport during the last decades propagation of sound in urban areas has become a major topic of outdoor sound propagation. This section presents the different methods used to study sound propagation in urban areas on different scales:

propagation around a single building and propagation in a district.

Single building scale

Most of the works on sound propagation around a single building concentrates on determining the Sound Pressure Level (SPL) in the shielded zone behind the building. It mainly consists in using ray tracing methods and the Uniform Theory of Diffraction (UTD) by [START_REF] Kouyoumjian | A uniform theory of diffraction for an edge in a perfectly conducting surface[END_REF] and its extension to multiple diffraction. Most recent works on this topic include the publications by [START_REF] Li | Absorbent parallel noise barriers in urban environments[END_REF] and [START_REF] Kirkpatrick | Sound propagation in the vicinity of an isolated building: an experimental investigation[END_REF]. The number of rays becoming incredibly high for multiple diffractions these methods are restricted to small number of diffracting edges and are hence not adapted to long-range sound propagation simulations. Moreover there is no diffraction theory available for high-amplitude wave propagation.

Street or district scale

Renterghem [START_REF] Renterghem | Numerical evaluation of sound propagation over green roofs[END_REF]; [START_REF] Renterghem | Parameter study of sound propagation between city canyons with a coupled FDTD-PE model[END_REF]; [START_REF] Heimann | Three-dimensional linearised Euler model simulations of sound propagation in idealised urban situations with win effects[END_REF] use Finite Difference Time Domain (FDTD) methods to study sound propagation over urban areas.

The geometries of buildings (overall dimensions, balconies, windows, etc.) and properties (absorbing surfaces, diffusely reflecting facades, etc.) are explicitly accounted for in the numerical methods. This yields accurate but time consuming simulations. [START_REF] Picaut | Numerical modeling of urban sound fields by a diffusion process[END_REF] uses a different approach to the FDTD method. An energy method (based on a diffusion equation, a mathematical extension to Sabine's concept of diffuse sound fields) is used to obtain the sound fields in a network of rectangular streets. This method has given good agreement with other solutions for propagation in a single rectangular street but the application of the method to a district seems difficult. Indeed, diffusion coefficients that depend on the facade roughness have to be experimentally determined. This can be a long and difficult process for large scale applications where there are several kinds of building facades.

Experimental studies of sound propagation in districts were performed by [START_REF] Picaut | Experimental study of sound propagation in a street[END_REF] and [START_REF] Thorsson | Noise levels on the shielded side in cities using a flat city model[END_REF].

Acoustically induced vibrations

Acousto-seismic coupling. Depending on the ground structure it is possible as a result of the poro-elasticity of ground for the above ground acoustic blast wave to couple to seismic waves, particularly at a frequency where the dispersive Rayleigh wave speed is the same as that of the speed of sound above the ground [see [START_REF] Sabatier | The interaction of airborne sound with the porous ground: The theoretical formulation[END_REF][START_REF] Madshus | Air-ground interaction in long range propagation of low frequency sound and vibration -field tests and model verification[END_REF].

Coupling to buildings vibrations. The first resonance frequencies associated with flexural vibration modes of buildings are generally in the range [0, 10] Hz, which is also the range of high energy of blast waves. It is then be possible for the blast wave to couple into building vibrations directly. This phenomenon is investigated, for example, in the paper by [START_REF] Sutherland | Sonic Booms And Building Vibration Revisited[END_REF].

Due to the high amplitudes involved, the attenuation associated with acoustic to seismic coupling at the ground or at building surfaces may be important and ideally should be taken into account. However, the study of ground and buildings vibrations being out of the scope of this work, it will further be assumed that the acousto-seimic coupling is weak.

Nonlinear effects in sound propagation

This section deals with the different phenomena occurring during high-amplitude wave propagation.

Wave steepening

The classical, linear theory of acoustics assumes that a wave speed only depends on its propagation medium. In the absence of any perturbation of the fluid medium, the sound speed is thus assumed to be constant with respect to space and time. While this assumption is valid for most of outdoor sound propagation applications, special care must be taken when high-amplitude waves are under concern.

Indeed, for plane waves in a homogeneous medium, we can derive from Euler equations

* v ≈ p ′ ρ 0 c 0 , (1.8a) c ≈ c 0 + (∂ p c) 0 p ′ , (1.8b) which combine to give c ≈ c 0 + βv , (1.9)
in which c 0 is the adiabatic speed of sound and c is the actual speed of sound. In Eq. (1.9), β is the coefficient of hydrodynamic nonlinearities defined by

β = 1 + 1 2 B A ,
(1.10) * Eqs. (1.8) are given as equations 11-1.5 in [START_REF] Pierce | Acoutics: An introduction to its physical principles and applications[END_REF] with

A = (ρ∂ ρ p) 0 , (1.11a) B = ρ 2 ∂ 2 ρ p 0 . (1.11b)
For a perfect gas, β is defined by

β = γ + 1 2 , (1.12)
γ being the ratio of specific heats. Eq. (1.9) means that the sound speed is directly proportional to the wave particle velocity v. This contribution is neglected for low-amplitude waves but may be important for high-amplitude waves. For example if a wave has a peak pressure of 5 kPa the sound speed increases by 15 m.s -1 at the positive pressure peak location * † . For a sine pulse the sound speed difference between its maximum and minimum pressures is hence 30 m.s -1 , about 10 % of the ambient sound speed. The positive pressure peak will eventually catch up the negative pressure peak, creating a shock. This phenomenon is known as wave steepening and is illustrated in Fig.

(1.9). The time waveform is shown at six different reduced distances σ, defined by

σ = x x , (1.13)
in which x is the shock formation distance. The shock formation distance mainly depends on the initial waveform. For one-dimensional wave and sinusoidal waveforms in a lossless medium, it is given by

x = ρ 0 c 2 0 βkP 0 , (1.14)
in which P 0 is the initial wave amplitude and k is the wavenumber. As one can see in Eq. (1.14) the shock formation distance x is frequency-dependent. At the distance x =x (or σ= 1), and for a lossless medium, the slope between negative and positive pressures is infinite.

Wave steepening is responsible for three effects of high-amplitude wave propagation:

◮ generation of harmonics, ◮ shock coalescence, ◮ and anomalous energy dissipation. * Assuming plane wave propagation, an adiabatic equation of state and a first-order approximation of the relation between pressure and density. † For air the coefficient of hydrodynamic nonlinearities β is approximately 1.2 under normal atmospheric conditions.

Figure 1.9: Illustration of wave steepening. As the wave propagates, a shock forms. For propagation in a lossless medium, the slope between negative and positive pressures is infinite when σ = 1. This figure is taken from [START_REF] Wochner | Numerical simulation of multi-dimensional acoustic propagation in air including the effects of molecular relaxation[END_REF].

Generation of harmonics

A high-amplitude sine waveform will decay into a shock wave (an N-shaped wave). The consequence is that harmonics of the fundamental frequency, which initially did not exist, will rise. In a lossless medium, the energy being conserved * , the energy contained in the fundamental frequency decreases to compensate for higher harmonics. The further the observation points is, the higher the harmonics amplitudes and and the lower the energy of the fundamental is. Hence, it is possible that at some distance one observes some frequencies that just did not exist in the source spectrum.

For simple waves one can easily derive analytical expressions for the wave Fourier components amplitudes, but for realistic signals with complex frequency content, the shock formation distance being dependent on the initial Fourier component amplitude and frequency, one has to use more advanced numerical methods. [START_REF] Webster | Experimental investigation of outdoor propagation of finite-amplitude noise[END_REF] used an array of loudspeakers to demonstrate and measure the harmonics generation phenomenon. As the observation points moves away from the source, the energy in higher frequencies increases, while the energy contained in the initial frequency decreases, so that total energy is conserved. This figure is taken from [START_REF] Wochner | Numerical simulation of multi-dimensional acoustic propagation in air including the effects of molecular relaxation[END_REF].

Shock coalescence

In the frame of the weak shock theory (see e.g. section 11.3 in [START_REF] Pierce | Acoutics: An introduction to its physical principles and applications[END_REF]), once a shock has formed due to wave steepening, its speed v sh can be approximated by .15) in which p ′ + and p ′ -are the pressure values in front of and behind the shock, respectively. Eq. (1.15) means that at some distance two shocks with different values of p ′ + + p ′ -move at different speeds. Hence, one shock may catch up the second: they may merge into a single shock. This phenomenon is called shock coalescence. When shocks coalesce, the number of zero-crossings typically decreases, increasing the characteristic time scale of the waveform.

v sh = c 0 + 1 2 β p ′ + + p ′ - ρ 0 c 0 , ( 1
One may thus observe transfer of energy from high to lower frequencies. shock coalescence: a high-amplitude Gaussian modulated sine pulse burst is shown at three different distances. The initial waveform is given at distance x 1 . At distance x 2 , the larger shocks which have a higher velocity catch up the weaker shocks, and finally, at distance x 3 , zero-crossings disappear and shocks coalesce into a single one. [START_REF] Gallagher | Experiments on the Nonlinear Characteristics of Noise Propagation from Low and Moderate Reynolds Number Supersonic Jets[END_REF] measured scale model jets and observed wave steepening and a shift of energy to lower frequencies which may have resulted from shock coalescence.

Anomalous energy dissipation

In the absence of shocks nonlinear effects do not change the energy associated with a wave. As seen in Section 1.2.1, energy will be rearranged over the frequency spectrum due to harmonics generation (lower to higher frequencies). For plane waves, the energy per unit area transverse to propagation direction is written as

E (t) = 1 ρ 0 c 2 0 +∞ -∞ p ′2 (x, t) dx .
(1.16)

If the propagation medium is lossless then the potential energy E (t) is independent of time.

However, if there is a shock, and since there is a discontinuity, then the integral in Eq. (1.16) has to be split in two parts. After a few algebraic manipulations [see [START_REF] Rudnik | On the attenuation of a repeated sawtooth shock wave[END_REF] we can obtain

d t E (t) = - βP 3 (t) 3ρ 2 0 c 3 0 , (1.17a) P (t) = P 0 1 + tβP 0 L 0 ρ 0 c 0 1/2 , (1.17b)
where P 0 and L 0 are the initial peak amplitude and the initial length of the positive and negative phases, respectively. The quantity P (t) being positive, the energy time derivative

d t E (t) is negative.
Hence, there is a dissipation process due to the presence of the shock, even in a lossless medium.

1.2.3 Propagation of high-amplitude waves in porous media: the Darcy law and the Forchheimer correction

Classical (linear) behavior of rigidly-framed porous materials has been presented in Section 1.1.3: the material introduces an attenuation and a delay to the reflected wave, modifying the Sound Pressure Level (SPL) at the receiver. In this section the nonlinear behavior of porous materials is briefly presented.

For steady flows it is well accepted that the pressure gradient in a porous medium with flow resistivity σ 0 is in most cases proportional to the particle velocity, such that

-∂ x p ′ = σ 0 v . (1.18)
Eq. (1.18), called Darcy's law, can be viewed as a definition of the static flow resistivity σ 0 .

While this relation is valid for low flow velocities, measurements have shown that the Darcy law is no longer valid for high flow rates. [START_REF] Forchheimer | Wasserbewegung durch Boden (Water movement through soil)[END_REF] proposed to correct the external force term with a term proportional to the square of the particle velocity, such that

-∂ x p ′ = σ 0 (1 + ξv) v , (1.19)
where ξ is the Forchheimer nonlinearity parameter. Eq. (1.19) implies that the global flow resistivity increases with increasing particle velocity. [START_REF] Nelson | Propagation of finite-amplitude sound in air-filled porous materials[END_REF] slightly modified this expression to adapt it to wave propagation:

-∂ x p ′ = σ 0 (1 + ξ |v|) v .
(1.20)

Eq. (1.20), compared to Eq. (1.19), accounts for the fact that the flow resistivity increases with the particle velocity absolute value. From Eq. (1.20) we can write a nonlinear flow resistivity σ defined by

σ = σ 0 (1 + ξ |v|) .
(1.21) [START_REF] Beavers | Non-Darcy flow through fibrous porous media[END_REF]; [START_REF] Beavers | Fluid flow through a class of highly deformable porous media, Part I: Experiments with air[END_REF]; [START_REF] Joseph | Nonlinear equation governing flow in a saturated porous medium[END_REF] have demonstrated the validity of the Darcy-Forchheimer law for various applications. Behavior of porous materials under high-amplitude acoustic excitation has been studied using the Darcy-Forchheimer law with some success by [START_REF] Kuntz | Attenuation of intense sinusoidal waves in air-saturated, bulk porous materials[END_REF] [and later by [START_REF] Wilson | Forchheimer-type nonlinearities for high-intensity wave propagation of pure tones in air-saturated porous media[END_REF]McIntosh & Lambert, 1990;Lambert & McIntosh, 1990;[START_REF] Auregan | Measurement of the nonlinear behaviour of acoustical rigid porous material[END_REF][START_REF] Umnova | High amplitude pulse propagation and reflection from a rigid porous layer[END_REF][START_REF] Umnova | Behavior of rigid-porous layers at high levels of continuous acoustic excitation: Theory and experiment[END_REF]. The main conclusion about high-amplitude wave reflection at porous boundaries is that the apparent porous layer rigidity increases with wave amplitude. However, note that [START_REF] Umnova | Behavior of rigid-porous layers at high levels of continuous acoustic excitation: Theory and experiment[END_REF] demonstrated that, depending on the layer thickness and material properties, the reflection coefficient associated with rigid-porous samples could decrease with increasing signal amplitude.

Blast wave propagation

This section gathers together bibliographical data on the acoustics of detonations. A synthesis of theoretical and experimental studies can be found in [START_REF] Bobin | Méthodes d'étude de la propagation d'une onde de souffle dans un réseau de galeries : synthèse des donnés bibliographiques (Study methods of blast wave propagation in a network of tunnels: synthesis of bibliographical data[END_REF]; [START_REF] Parmentier | Synthèse des résultats expérimentaux relatifs aux détonations d'explosifs sphériques[END_REF].

Similarity law

The hypothesis of spherical propagation in a homogeneous atmosphere is made. Under this hypothesis the similarity law is verified [see [START_REF] Kinney | Explosive shocks in air[END_REF][START_REF] Swisdak | Explosions effects and properties. Part 1. Explosion effects in air[END_REF]. For two explosions with energies W 1 and W 2 , we have

K 2 = K 1 W 2 W 1 1/3 , (1.22)
where K 1 and K 2 can be ◮ either distances from the explosion where one has the same overpressure peaks max p (t),

◮ either times of shock arrival t 0 , ◮ either positive phase durations t + -t 0 , ◮ or the overpressure impulses defined by I =

t + t 0 p (t) dt.
If ambient pressure P 0 and temperature T 0 are different from an explosion to another one, the similarity law Eq. (1.22) must be corrected according to

K 2 = K 1 W 2 W 1 P 0 1 P 0 2 1/3 T 0 1 T 0 2 1/2
.

(1.23)

The application of Eq. (1.23) allows us (as long as the atmosphere is homogeneous and the explosions energies are known) to restrain the study to the case of a single explosion. Usually, charts and curves are given for explosions occurring in standard atmospheric conditions at sea level.

TNT equivalence

For explosive materials of different compositions the source energy can be measured in terms of kilograms TNT equivalent. This means that 1 kg of a certain explosive material is equivalent to 1 kg of TNT times an equivalence factor. These factors are most of the time determined through experiments [see [START_REF] Swisdak | Explosions effects and properties. Part 1. Explosion effects in air[END_REF][START_REF] Hyde | User's guide for microcomputer programs CONWEP and FUNPRO[END_REF]. Tab. (1.1) lists some equivalence factors found in the literature for some common explosives for the positive pressure peak. Note that the equivalence factor may depend on the mass of material and that they may be different for different wave properties (e.g. maximum pressure/positive phase duration).

Table 1.1: TNT equivalence factors for some common explosives. Figures are taken from [START_REF] Parmentier | Synthèse des résultats expérimentaux relatifs aux détonations d'explosifs sphériques[END_REF].

Material

Equivalence factor for pressure TNT 1

Dynamite (Nitroglycerin) 0.9 C4 1.37

Analytical solutions

Theoretical studies performed during the last five decades aimed at overcoming the lack of experimental data. The first attempt to develop a theory for blast wave propagation was done by von Neumann [1961]. In the point source theory all the energy is concentrated at a single point where the pressure is infinite. Gas dynamics equations allow analytical solutions to be derived for the pressure field at any distance from the point source. This theory assumes that the air is a perfect gas for which specific heats are constants and that the propagation medium is homogeneous.

A variant of this method was developed by [START_REF] Taylor | The formation of a blast wave by a very intense explosion[END_REF]. Taylor supposes that the energy is concentrated in a sphere with radius R 0 and thus solves the infinite pressure issue and makes the method more adapted to numerical calculations. These analytical solutions have been used for comparisons with nuclear tests in New Mexico; good agreement has been found between measurements and theory. [START_REF] Laumbach | A point explosion in a cold exponential atmosphere[END_REF] solved the gas dynamics equations in a cold atmosphere with exponentially varying volumetric mass. The agreement with the results of [START_REF] Taylor | The formation of a blast wave by a very intense explosion[END_REF] is very good for a homogeneous atmosphere. [START_REF] Kinney | Explosive shocks in air[END_REF] derived analytical solutions for blast wave propagation from spherical charges. The model gives the pressure history p (t) at a given reduced distance Z in a homogeneous medium, for an equivalent charge W of TNT. To the author's knowledge, the Kinney-Graham (KG) model is the most accurate and complete one. The Kinney-Graham (KG) model is presented in Appendix B.

These theoretical studies are essential to the understanding of the phenomenon of blast wave propagation. However, these solutions are most of the time obtained in the frame of very restrictive hypotheses (homogeneous medium, point sources, etc.) and cannot be used for simulations of high-amplitude wave propagation in complex media.

Experimental studies

Institut Saint Louis (ISL) has performed many experiments for spherical explosives in free field [see for example Froböse, 1968;[START_REF] Froböse | Mesures dans les ondes de choc engendrées par la detonation d'une sphère en explosifs plastique[END_REF][START_REF] Froböse | Production et enregistrements d'ondes de choc engendrées par l'éclatement d'une étincelle ou par détonation en vue de l'étude de leurs effets sur l'audition du cobaye[END_REF]. These experiments have been used for comparisons with analytical solutions, or to study nonlinear phenomena but measurements were performed at short ranges from the sources.

Van der Eerden & Védy [2005]; [START_REF] Attenborough | Overview of the theoretical development and experimental validation of blast sound-absorbing surfaces[END_REF][START_REF] Attenborough | Blast sound absorbing surfaces[END_REF] numerous acoustic and seismic data to be gathered at large distances (more than 400 km from the explosion). The Misty Picture experiment is still often used as a benchmark for blast wave propagation codes [see [START_REF] Gainville | Ray tracing for long range atmospheric propagation of infrasound[END_REF][START_REF] Piserchia | Numerical modeling of infrasound propagation at very long distance[END_REF], especially for applications to propagation in inhomogeneous media (strong refraction effects, turbulence).

Numerical methods

This section presents the most commonly used models to simulate high-amplitude wave propagation in inhomogeneous environments. Specifically, general time-domain methods and models based on a paraxial approximation are detailed.

General time-domain algorithms

Navier-Stokes equations and Euler equations for inviscid fluid are used for aeroacoustic simulations in Computational Fluid Dynamics (CFD). With the increase of the available compu-tational power and the advances in numerical methods Euler equations have become popular for the study of wave propagation over moderate distances.

The Euler equations are composed of a continuity equation, N equations of conservation of momentum, where N is the number of spatial directions under consideration, and one or more equations for energy (or entropy) conservation.

Various constitutive equation sets have been derived those last years. Most recent models

[see for example [START_REF] Wochner | Numerical simulation of multi-dimensional acoustic propagation in air including the effects of molecular relaxation[END_REF][START_REF] Sparrow | A numerical method for general finite amplitude wave propagation in two dimensions and its application to spark pulses[END_REF] are able to simulate absorption from shear and Bulk viscosity, thermal conductivity, and molecular relaxation processes, and account for realistic atmospheric conditions.

For a two-dimensional domain with coordinates (x, z), the equation set for propagation in a lossless medium is

∂ t ρ T + ∂ x (ρ T u) + ∂ z (ρ T w) = 0 , (1.24a) ∂ t (ρ T u) + ∂ x ρ T u 2 + ∂ z (ρ T uw) = -∂ x p T , (1.24b) ∂ t (ρ T w) + ∂ x (ρ T uw) + ∂ z ρ T w 2 = -∂ z p T , (1.24c)
∂ t (ρ T e 0 ) + ∂ x (ρ T ue 0 ) + ∂ z (ρ T we 0 ) = -∂ x (p T u) -∂ z (p T w) , (1.24d)
where e 0 is the energy per unit mass. The energy equation Eq. (1.25) and the ideal gas law Eq. (1.26) close the equation system, such that

ρ T e 0 = ρ T c v T + ρ T |V| 2 2 , (1.25)
and

p T = ρ T RT , (1.26)
in which T is the gas temperature, c v is the specific heat capacity at constant volume and R is the gas constant.

For the sake of brevity detailed notations are not given here. This model requires advanced numerical methods to be solved, one of the most commonly used being the finite-difference method. [START_REF] Wochner | Numerical simulation of multi-dimensional acoustic propagation in air including the effects of molecular relaxation[END_REF] used a third-order accurate Runge-Kutta (RK) scheme for time discretization and a Weighted Essentially Non Oscillatory (WENO) scheme for space. Details on these finite-difference schemes can be found in Appendix C. The RK scheme requires the calculation of two fields at each time step, and for each of them, the evaluation of every derivative appearing in Eqs. (1.24) has to be performed at each grid position. Although the numerical model presented above yields very accurate solutions for nonlinear sound propagation, its computational complexity makes it difficult to use for long-range propagation applications.

Paraxial approximations: the KZK equation and the NPE model

With hypotheses of weak nonlinearities and a paraxial approximation, the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation [see [START_REF] Zabolotskaya | Quasi-plane waves in the non-linear acoustics of confined beams[END_REF][START_REF] Kuznetsov | Equations of nonlinear acoustics[END_REF] can be derived from Euler equations * . The KZK equation writes

∂ z p = βp ρ 0 c 3 0 ∂ τ p + δ 2c 3 0 ∂ 2 τ p + c 0 2 τ -∞ ∂ x p 2 dτ , (1.27)
where τ = t-z/c 0 and δ accounts for absorption from air. In its original formulation the KZK equation does not include the effects from spatially-varying sound speed. However, [START_REF] Blanc-Benon | Propagation of finite-amplitude sound through turbulence: Modeling with geometrical acoustics and the parabolic approximation[END_REF] modified the equation so that it can account for sound speed variations and turbulence [see also [START_REF] Ganjehi | Ondes de choc acoustiques en milieu hétérogène, des ultrasons au bang sonique (Acoustic shock waves in heterogeneous medium, from ultrasounds to sonic bang[END_REF].

The Nonlinear Parabolic Equation (NPE) is derived using similar assumptions. The NPE has first been developed by [START_REF] Mcdonald | Time domain formulation for pulse propagation including nonlinear behaviour at a caustic[END_REF] and has been successfully used for underwater acoustics simulations [see [START_REF] Castor | Long-range propagation of finite-amplitude acoustic waves in an ocean acoustic waveguide[END_REF] and blast wave propagation in air

[see van der Eerden & Védy, 2005;[START_REF] Attenborough | Overview of the theoretical development and experimental validation of blast sound-absorbing surfaces[END_REF][START_REF] Attenborough | Blast sound absorbing surfaces[END_REF][START_REF] Leissing | Nonlinear outdoor sound propagation -A numerical implementation and study using the nonlinear progressive wave equation[END_REF]. The NPE is based on the resolution of a nonlinear wave equation over a moving window that surrounds the wavefront. While reducing domain size, and thus computational cost, the moving window principle prevents backward propagation. For the derivation of the original NPE model we refer the reader to articles by [START_REF] Mcdonald | A tutorial on the nonlinear progressive wave equation (NPE) -Part 1[END_REF] or [START_REF] Caine | A tutorial on the nonlinear progressive wave equation (NPE) -Part 2. Derivation of the three dimensional cartesian version without use of perturbation expansions[END_REF].

The NPE model for a two-dimensional domain with Cartesian coordinates (x, z) is based on the following equation:

D t R = - c 0 R 2r -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx .
(1.28)

The ambient sound speed is c 0 while c 1 is the sound speed perturbation in the window, i.e. c 1 = c (x, z) -c 0 , where c (x, z) is the spatially-dependent sound speed. The overdensity variable R= ρ ′ /ρ 0 , in which ρ ′ the acoustic density perturbation and ρ 0 is the ambient medium density, is dimensionless. The moving-window operator D t is defined by

D t = ∂ t + c 0 ∂ x .
(1.29)

The assumptions used to derive this model from the Euler equations are ◮ weak nonlinearities, ◮ weak sound speed perturbations,

◮ and propagation along a main direction.

Eq. (1.28) can be used to propagate weak shocks over moderate distances within a domain with spatially-varying sound speed.

Various modifications and additions to this original model have been made during the past two decades: spherical and cylindrical coordinate system versions [see Too & Ginsberg, 1992a] and high-angle formulation [see [START_REF] Mcdonald | High-angle formulation for the nonlinear progressive wave equation model[END_REF] have been developed, and Too & Lee [1995a] extended the NPE with an additional term to account for thermoviscous effects. Propagation in multiple media has also been successfully studied using this model by Ambrosiano et al. The main feature of paraxial approximations is their capability to perform long-range simulations with low numerical costs. Indeed, the use of single variable one-way wave equation makes them suitable tools for large scale applications. Limitations of these paraxial methods derive from the underlying assumptions: the hypothesis of weak nonlinearities prevents the use of these methods in the near field of the source and moreover, the moving-window technique restricts the initial acoustic field to transient signals.

Blast waves temporal and spectral characteristics

In the case of an explosion the amplitude of the wave in the near field of the source is obviously highly dependent on the explosive charge. At moderate distances from the source (a few hundred meters), the positive amplitude peak is on the order of some kPa while the negative amplitude peak is lower, on the order of a few hundred Pa. Note that due to nonlinear dissipation processes the wave amplitude does not decay at a geometrical rate, even in a homogeneous medium (see Section 1.2.2). For a charge of 30 T of TNT * at a distance of 200 m (or 1000 m) from the source the positive and negative peak pressures are 18.9 kPa (or 2.6 kPa), and -5.7 kPa (or -236 Pa) (see Tab. (1.2) and Fig.

(1.12)).

Time scales of acoustic waves from explosions are on the order of a few milliseconds. The positive phase duration is usually shorter than the negative one and as the wave propagates further from the source, the difference becomes more and more important. For instance, for an acoustic wave due to a charge of 30 T of TNT, the positive and negative phase durations are 2.9 ms and 34 ms at a distance of 200 m, and 4.1 ms and 91 ms at a distance of 1000 m (see Tab. (1.3) and Fig.

(1.12)).

Generally speaking, the frequency where the amplitude is the highest is on the order of a few dozen Hertz. As an example, a charge of 30 T of TNT at a distance of 200 m produces a wave whose spectrum has a peak amplitude at a frequency of approximately 30 Hz. The waveform being highly impulsive, the bandwidth where the spectrum has a high energy is often broad.

A -3 dB bandwidth yields a frequency range [0, 70] Hz (see Fig.

(1.12)). As the negative phase duration becomes more and more important as the wave propagates the energy in the high-frequency range decreases with distance. At a distance of 1000 m the most energetic frequency is approximately 5 Hz.

Conclusions and outline of the document

In this section, conclusions on blast wave propagation in complex environments are first given in Section 1.4.1. With the help of these remarks on the features that are of interest for this research the proposed work is described in Section 1.4.2 and finally, the outline of this document is given in Section 1.4.3.

Conclusions

Features of nonlinear sound propagation in complex environments Section 1.1 has summarized features of sound propagation that are of special importance when propagation of acoustic waves over large distances is under interest. Within the propagation domain one must take into account atmospheric conditions, since refraction effects can largely modify sound rays paths: thanks to an inversion layer the sound from the Buncefield explosion was heard in France, more than 150 km away from the depot.

Considering the frequencies studied, dissipation effects will be on the order of a few dB (see In Section 1.2 nonlinear effects on sound propagation have been reviewed. The combined effect of wave steepening and shock coalescence induces a shift of energy from the original frequency band to the high frequency band due to steepening, and from the original frequency band to lower frequency band due to shock coalescence. Dissipation at the shock introduces some additional absorption that has to be taken into account.

The combined effects of meteorological conditions, of ground topography and impedance, and of high-amplitude effects can yield extremely complex situations that can only be handled through the use of advanced models and numerical methods.

On the size of the propagation domain and computational times

The main objective of this work is to develop a model for high-amplitude sound propagation in complex media, with application to the study of blast waves from explosions. This involves simulations in large domains whose typical distances are 10 km wide and 2 km high, with wave frequencies up to one hundred Hertz. With a small-signal sound speed of 340 m/s, this yields a domain size of approximately 3 000 by 600 wavelengths. Although no numerical method has yet been chosen, one can choose as a first guess a spatial resolution of 10 points/λ. This yields grids with 3.10 4 × 6.10 3 nodes, giving 1.8 billion grid points for two-dimensional simulations.

Hence, the numerical complexity of the available numerical methods is a fundamental choice criteria.

It has been chosen to restrict simulations to two-dimensional (or three-dimensional with azimuthal symmetry) calculations. Realistic (three-dimensional) cases can be approximated by repeated 2.5D problems in different compass directions. This allows the problem size to be considerably reduced compared to full problems. Note that this hypothesis can become problematic in case of strong crosswinds. Even with two-dimensional simulations, the problem size is still huge and the effective problem size has to be reduced. This can be done in two ways:

◮ either by reducing the number of points in the calculation grid, ◮ or by reducing the model complexity.

Reducing the number of grid points. Two methods can be used to lower the number of calculation points in the grid:

◮ Adaptive Mesh Refinement (AMR)

◮ and moving-window methods.

Both rely on the simple principle that it is not necessary to mesh the domain where there is no acoustic perturbations.

AMR techniques use a grid hierarchy: first the top-level grid is created with a coarse mesh.

For each of its cell an algorithm determines if the solution is accurate enough. If not the cell is refined with a finer sub-grid, and so on until the desired level of accuracy is reached. The main difficulty with AMR methods is to deal with boundary conditions between sub-grids of different levels.

Another way to decrease the grid size is to use a moving-window technique, in which a window moving with the wave surrounds the waveform. The implementation of such moving-windows is straightforward: the grid remains Cartesian and no special boundary conditions are needed.

Note that the NPE uses a frame-following coordinate formulation and hence naturally uses the moving-window principle.

Since we deal in this work with impulsive and finite-length signals, using any grid size reduction method should be extremely efficient.

Reducing model complexity. Another way of lowering the computational effort is to reduce the model complexity. Indeed, a model such as the KZK equation or the NPE model is based on a single-variable one-way wave equation. There is only one unknown and at most first-order time derivatives. Oppositely, the two-dimensional Euler equations are composed for the simplest models of five equations with five unknowns for the simplest models, and the calculation of the different fields requires the evaluation of a dozen spatial derivatives.

Although Euler equations solutions provide very accurate solutions, the associated numerical effort is often prohibitive in the context of long-range sound propagation applications.

Proposed work and methodology

Proposed work. A feature of long-range sound propagation outdoors that has been kept apart until now is propagation over urban environments. Specifically, we are interested in the effect of the surface irregularities (buildings) on the acoustic field above the urban layer.

This problem could be studied with deterministic numerical models where buildings geometries are explicitly given. For example, ray tracing methods or models based on the Euler equations allow the environmental context to be taken into account and would technically be suited for this application. However, for long-range propagation applications, these methods suffer from their numerical complexity. Moreover, the suitability of these models for the application under interest can be questioned: considering the high complexity of the real system, the model approximations (e.g. using the uniform theory of diffraction for ray-tracing methods) and the uncertainties on the model parameters (e.g. buildings geometries), the computational model could be improved introducing a probabilistic model.

In this work, a different approach is proposed. It consists in using a simplified model for the urban and atmospheric layers. Propagation is modeled with a Nonlinear Parabolic Equation (NPE). Since the calculation domain is limited to a small area around the signal, computational cost is generally reduced compared to Euler equations methods. For the propagation in the urban layer, a NPE for porous ground layers (which has to be developed) will be used.

Therefore the simplified model is composed of two domains:

◮ the atmospheric layer, where propagation is modeled by a NPE for air,

◮ and the urban city layer, where propagation is modeled by a NPE for porous ground surfaces.

Equations to couple the two domains complete the simplified model of sound propagation over urban cities.

Using this model as a predictive model for the real system will only allow one to get a raw approximation of the propagation phenomenon. Hence a probabilistic approach of uncertainties will be developed and used to enhance the model capabilities.

The work presented in this document can be separated in two main tasks. Concerning the porous layer, the model uses a minimal parametrization in which the layer is assumed to be equivalent to a continuous fluid medium. A wave causes a vibration of air particles contained in the ground pores, while the ground frame does not vibrate. Next, a first-order boundary interface condition to link the two Nonlinear Parabolic Equation (NPE) models (air and porous layers) is developed, which finalizes the construction of the deterministic NPE model for sound propagation in complex media.

Methodology for the development of a stochastic NPE model. The objective is to develop a stochastic model for sound propagation over urban cities using a NPE model originally designed for sound propagation over porous ground layers, in which the urban city is taken into account through independent random porous layer parameters, denoted Γ, Λ, Θ.

The probability distributions of these random variables depend on a parameter w. The model output is the pressure at the receiver which is noted P r (ω). A second propagation model, in which the urban city is explicitly accounted for, is used to provide reference solutions. Given a parameter u which characterizes the urban city geometry and several probability models, a urban city realization is generated, and then used in the reference model to obtain the pressure P exp r (ω) at the receiver. To construct the stochastic sound propagation model, for a given parameter u, one has:

Step 1: to construct the probability models of the urban city geometrical parameters. The probability models are determined with the help of Information Theory [see [START_REF] Shannon | A mathematical theory of communication[END_REF] and the Maximum Entropy Principle [see [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. Once the probability models of the geometrical parameters are determined, different city realizations corresponding to a given parameter u can be generated.

Step 2: to use the city realizations generated in Step 1 and to perform simulations with the reference model in order to obtain statistical information on the model output

P exp r (ω).
Step 3: to construct the probability models of the stochastic NPE model random parameters Γ, Λ and Θ. These models are determined using the same method as in Step 1 and depend on parameter w.

Step 4: using the outputs from the reference model P exp r (ω) (Step 2), to identify parameter w corresponding to the given parameter u previously fixed. This identification is done by solving an inverse stochastic problem: the "distance" between P r (ω; w)

and P exp r (ω) is minimized so that the optimal parameter w opt is obtained. Once w opt is determined, the NPE stochastic model can be used to study nonlinear wave propagation over urban cities.

Main features of high-amplitude wave propagation outdoors, such as dissipation effects, refraction effects or hydrodynamic nonlinearities, are naturally present in NPE models. The constructed stochastic model, thanks to its low numerical cost, could thus be used to study nonlinear wave propagation in complex environments. The urban city parameter u is given.

Simplified wave propagation model

Outline of this document

The document is organized as follows.

Chapter , 1987;[START_REF] Mcdonald | A tutorial on the nonlinear progressive wave equation (NPE) -Part 1[END_REF] while the second uses some assumptions about the air particles motion [see [START_REF] Caine | A tutorial on the nonlinear progressive wave equation (NPE) -Part 2. Derivation of the three dimensional cartesian version without use of perturbation expansions[END_REF].

In this section a derivation procedure that yields the original NPE model is presented step-bystep. Although the derivation presented here is different from what can be found in previously cited papers, it is highly inspired from them.

Notations and variables definition

The propagation domain considered is two-dimensional with main axes x (horizontal direction) and z (vertical direction). The sound speed c (x, z) is allowed to vary with position but is stable with respect to time and total density ρ T and total pressure p T variables are noted as follows:

ρ T = ρ 0 + ρ ′ , (2.1a)
p T = p 0 + p ′ , (2.1b)
in which ρ 0 and p 0 are ambient density and ambient pressure, respectively, and ρ ′ and p ′ are acoustic perturbations of these quantities. The components of the flow velocity vector V are u and w, which are the flow velocities in the x-and z-directions, respectively. Ambient medium properties are assumed to be fixed in space and time and without acoustic perturbations the flow velocities are equal to zero. Partial derivation with respect to the variable i is noted ∂ i .

Proposed derivation

In a two-dimensional coordinate system and in a lossless media * the continuity and conservation of momentum equations are

† ∂ t ρ T = -∂ x (ρ T u) -∂ z (ρ T w) , (2.2a) ∂ t (ρ T u) = -∂ x p T + ρ T u 2 -∂ z (ρ T uw) , (2.2b) ∂ t (ρ T w) = -∂ z p T + ρ T w 2 -∂ x (ρ T uw) .
(2.2c)

The

∂ 2 t ρ T = ∂ 2 x p T + ρ T u 2 + ∂ 2 z p T + ρ T w 2 + 2∂ x ∂ z (ρ T uw) . (2.3)
We then make the assumption that the propagation is predominant in the x-direction and thus eliminate nonlinear terms in z-derivatives in Eq. ( 2.3). Specifically, terms ∂ 2 z ρ T w 2 and ∂ x ∂ z (ρ T uw) are discarded. This yields

∂ 2 t ρ T = ∂ 2 x p T + ρ T u 2 + ∂ 2 z p T .
(2.4)

To reduce Eq. ( 2.4) to a single variable equation the total pressure p T is substituted by a second-order expansion in ρ ′ from an assumed adiabatic equation of state defined by

p T = p 0 + c 2 ρ ′ + c 2 γ -1 2ρ 0 ρ ′ 2 , (2.5)
in which γ is the ratio of specific heats. It is furthermore assumed that the sound speed c (x, z) can be written

c (x, z) = c 0 + c 1 (x, z) , (2.6)
where c 0 is a constant and c 1 (x, z) is a spatially varying perturbation. Substituting Eq. (2.6) in Eq. (2.5) gives

p T = p 0 + (c 0 + c 1 ) 2 ρ ′ + (c 0 + c 1 ) 2 γ -1 2ρ 0 ρ ′ 2 . (2.7)
With this new expression for the total pressure p T Eq. (2.4) becomes

∂ 2 t ρ T = ∂ 2 x (c 0 + c 1 ) 2 ρ ′ + (c 0 + c 1 ) 2 γ -1 2ρ 0 ρ ′ 2 + ρ T u 2 + ∂ 2 z (c 0 + c 1 ) 2 ρ ′ . (2.8)
In order to further simplify Eq. (2.8) the assumption that the sound speed perturbation c 1 (x, z) is small compared to c 0 , i.e. c 1 (x, z)≪c 0 , is made, and hence the term ∂ 2 z c 2 1 ρ ′ in Eq. (2.8) can be eliminated. This yields

∂ 2 t ρ T = ∂ 2 x (c 0 + c 1 ) 2 ρ ′ + (c 0 + c 1 ) 2 γ -1 2ρ 0 ρ ′ 2 + ρ T u 2 + c 2 0 ∂ 2 z ρ ′ .
(2.9)

The last step to reduce Eq. (2.9) to a single-variable equation is to eliminate the variable u, the flow velocity in the x-direction. A first-order result for plane waves is used, i.e. the velocity u is written

u = c 0 ρ ′ ρ 0 .
(2.10) Substituting Eq. (2.10) in Eq. (2.9), and keeping only terms of order up to two in the xderivative term gives

∂ 2 t ρ T = ∂ 2 x c 2 0 + 2c 0 c 1 ρ ′ + c 2 0 γ + 1 2 ρ ′ 2 ρ 0 + c 2 0 ∂ 2 z ρ ′ .
(2.11)

An important remark can be done at this stage of the derivation. Indeed, if one assumes low amplitude waves (i.e. neglects the nonlinear term in Eq. (2.11)) and assumes that the propagation medium is homogeneous (i.e. c 1 = 0), one obtains

∂ 2 t ρ T -c 2 0 ∇ 2 ρ T = 0 , (2.12)
which is the classical, linear wave equation.

Eq. (2.11) now contains only one variable, the density perturbation ρ ′ . The main principle of NPE models is the simulation of sound propagation in a frame surrounding the wavefront and moving at speed c 0 . The propagation is hence assumed to be outgoing only and a framefollowing formulation of Eq. (2.11) is sought. A moving-frame operator D t is introduced. It is defined by (2.13) and taking the second time derivative yields

D t = ∂ t + c 0 ∂ x ,
∂ 2 t = D 2 t -2c 0 D t ∂ x + c 2 0 ∂ 2 x . (2.14)
Provided one is distant from the source and the solution can be written in the form

ρ (t, x, z) = F (x -c 0 t) G (x, z) , (2.15)
then it can be shown that if the (first-order) parabolic approximation holds [see [START_REF] Caine | A tutorial on the nonlinear progressive wave equation (NPE) -Part 2. Derivation of the three dimensional cartesian version without use of perturbation expansions[END_REF] then

D 2 t = 0 . (2.16)
Rewriting Eq. (2.11) with the moving-frame operator yields a one-way single-variable wave equation, defined by (2.17) in which β= γ + 1 2 is the coefficient of nonlinearity. Integrating Eq. (2.17) with respect to the variable x and rearranging the equation gives a first version of the Nonlinear Parabolic Equation (NPE) model, such that

-2c 0 D t ∂ x ρ ′ = ∂ 2 x 2c 0 c 1 ρ ′ + c 2 0 β ρ ′ 2 ρ 0 + c 2 0 ∂ 2 z ρ ′ ,
D t ρ ′ = -∂ x c 1 ρ ′ + c 0 β 2 ρ ′ 2 ρ 0 - c 0 2 ∂ 2 z ρ ′ dx .
(2.18)

Eq. (2.18) is then written with a dimensionless density perturbation R= ρ ′ ρ 0 ; this yields

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx . (2.19)
Note that equivalently, one could use the dimensionless quantity 2.18) gives an error consistent to the order of accuracy sought.

Q = p ′ ρ 0 c 2 0 [see Ambrosiano et al. , 1990]. Since Q = R + O ǫ 2 , replacing R by Q in Eq. (
The above equation Eq. (2.19) is the NPE model as presented in papers by [START_REF] Mcdonald | Time domain formulation for pulse propagation including nonlinear behaviour at a caustic[END_REF]; [START_REF] Mcdonald | A tutorial on the nonlinear progressive wave equation (NPE) -Part 1[END_REF]; [START_REF] Caine | A tutorial on the nonlinear progressive wave equation (NPE) -Part 2. Derivation of the three dimensional cartesian version without use of perturbation expansions[END_REF]. However note that the cylindrical spreading term c 0 R/2r, where r is the distance from the source, has been dropped for clarity.

The derivation of the original NPE model without the use of the perturbation expansion method can be found in Appendix D.

Properties of the Nonlinear Parabolic Equation model

Assumptions used for the derivation of the NPE model

To reduce the NPE Eq. (2.19) to a single-variable one-way wave equation, several hypothesis are made. Listing these simplifications will help highlighting the model limitations.

◮ Main propagation along the x-direction: this hypothesis allows one to keep only linear terms in z-derivative in Eq. (2.3), and thus eliminate the vertical flow velocity variable w from the model.

◮ Weak sound speed perturbations: assuming c 1 (x, z)≪c 0 lets further simplify the model by keeping only the first-order term in z-derivative (the term c 2 0 ∂ 2 z ρ ′ ) in Eq. (2.8).

◮ Weak nonlinearity: the pressure is replaced by a second-order expansion in ρ ′ , the flow velocity u by a first-order approximation. The NPE model is hence limited to the simulation of waves with moderate overpressure values.

◮ Paraxial approximation: the use of a moving-frame operator D t implies a one-way wave propagation (no back-scattering). The use of the NPE model is hence restricted to far-field simulations of finite-length signals.

Modularity of the NPE model

One can identify three different differential operators in Eq. (2.19):

◮ the term -∂ x (c 1 R) accounts for refraction effects due to spatially-varying meteorological conditions,

◮ the term -c 0 β 2 ∂ x R 2 accounts for nonlinear effects, ◮ and the last term - c 0 2 ∂ 2 z R
-∂ x c 0 β 2 R 2 .

Relation to the frequency-domain Parabolic Equation (PE)

Dropping the nonlinear term in Eq. (2.19) yields

D t R = -∂ x (c 1 R) - c 0 2 ∂ 2 z R dx , (2.20)
and substituting R = e i(kx-ωt) f (x, z) in Eq. (2.20), in which ω and k are the wave pulsation and wavenumber, respectively, and f (x, z) is slowly varying in the x-direction compared to the exponential term, and retaining only dominant terms (terms from the rapid oscillations of the exponential) gives

∂ x f = -ik c 1 c 0 f + i 2k ∂ 2 z f .
(2.21) [START_REF] Mcdonald | Time domain formulation for pulse propagation including nonlinear behaviour at a caustic[END_REF]. We refer the reader to a paper by [START_REF] Tappert | The parabolic approximation method[END_REF] for pioneering work on PE for linear underwater acoustics applications. A synthesis of Parabolic Equations (PE) methods can be found in the paper by [START_REF] Lee | Parabolic equation development in the recent decade[END_REF].

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx ◮ Mowing-

Relation to the inviscid Burgers equation

The inviscid Burgers equation is * (2.22) where τ = t -x/c 0 introduces a retarded time frame. The above equation is very similar to the NPE model Eq. ( 2.19) in one dimension, written as This equation, without the thermoviscous absorption term, is *

∂ x p ′ = β 2ρ 0 c 3 0 ∂ τ p ′2 ,
D t R = -∂ x c 1 R + c 0 β 2 R 2 . ( 2
∂ x p = 1 c 0 ∂ t c 1 p + β 2ρ 0 c 0 p 2 + c 0 2 t -∞ ∂ 2 z p dt . (2.24)
This equation includes hydrodynamic nonlinearities and the effects of sound speed variations with altitude. The NPE and KZK equations are similar with the roles of time and distance reversed. These models mainly differ by their numerical implementation: the NPE uses an accurate finite-difference scheme (Flux Corrected Transport (FCT) algorithm, see Chapter 4)

to evaluate the nonlinear operator [see [START_REF] Castor | Long-range propagation of finite-amplitude acoustic waves in an ocean acoustic waveguide[END_REF].

Extensions to the original Nonlinear Parabolic Equation model

The formulation of the NPE defined by Eq. (2.19) only accounts for sound speed variations and nonlinear effects. The NPE model can be extended to handle geometrical spreading and dissipative effects from thermoviscous absorption.

Geometrical spreading

For the case of azimuthal symmetry the NPE model contains an additional operator accounting for geometrical spreading. The NPE model becomes

D t R = - c 0 R 2r -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx , (2.25)
in which r = c 0 t + x is the distance from the source and x denotes the range in the moving window.

Thermoviscous effects

This section describes the procedure to include thermoviscous effects in the NPE propagation model described in Section 2.1. The derivation is based on Too & Lee [1995a]. * Eq. (2.24) is given as equation ( 25) in Blanc-Benon et al. [2002] In a two-dimensional coordinate system the continuity and conservation of momentum equations including thermoviscous effects are [see [START_REF] Wochner | Numerical simulation of finite amplitude wave propagation in air using a realistic atmospheric absorption model[END_REF]]

∂ t ρ T = -∂ x (ρ T u) -∂ z (ρ T w) , (2.26a) ∂ t (ρ T u) = -∂ x p T + ρ T u 2 -∂ z (ρ T uw) + µ B ∂ 2 x u + ∂ x ∂ z w + µ (∂ x Φ xx + ∂ z Φ xz ) , (2.26b) ∂ t (ρ T w) = -∂ z p T + ρ T w 2 -∂ x (ρ T uw) + µ B ∂ 2 z w + ∂ z ∂ x u + µ (∂ z Φ zz + ∂ x Φ zx ) , (2.26c)
in which µ and µ B are the fluid shear and Bulk viscosities, respectively and Φ ij represents the rate of shear tensor which is defined by

Φ ij = ∂ j v i + ∂ i v j - 2 3 ∇ • Vδ ij , (2.27)
where the v i are components of the particle velocity vector V and δ is the Kronecker delta.

The derivation procedure that yields a NPE model with thermoviscous effects is similar to the one presented in Section 2.1.2: Eq. (2.26a) is differentiated with respect to time, Eqs. (2.26b, 2.26c) are substituted in the newly developed expression, and only terms of first order in z-derivatives are kept. Furthermore, it can be proved that to a certain order of accuracy, and under the assumptions made, w= 0. The proof of this statement is voluntarily omitted here: it is redundant to the derivation presented in Section 3.2. These assumptions yield

∂ 2 t ρ T = ∂ 2 x p T + ρ T u 2 + ∂ 2 z p T - 4 3 µ + µ B ∂ x ∇ 2 u .
(2.28)

The following equation of state is used * for the total pressure p T :

p T = p 0 + c 2 ρ ′ + c 2 γ -1 2ρ 0 ρ ′2 -κ 1 c v - 1 c p ∇ • V , (2.29)
where κ is the thermal conductivity and c v and c p are the specific heats. Substituting Eq. (2.29) in Eq. (2.28) and following the procedure described in Section 2.1.2 yields

∂ 2 t ρ T = ∂ 2 x c 2 0 + 2c 0 c 1 ρ ′ + c 2 0 γ + 1 2 ρ ′2 ρ 0 + ∂ 2 z c 2 0 ρ ′ -c 0 δ th ∂ x ∇ 2 ρ ′ , (2.30)
where δ th is the sound diffusivity [see [START_REF] Lighthill | Waves in fluids[END_REF]Too & Lee, 1995a] and is defined by

* δ th = ρ -1 0 4 3 µ + µ B + κ 1 c v - 1 c p .
(2.31)

A NPE model that includes thermoviscous effects can be obtained in introducing the movingframe operator and rearranging Eq. (2.30). One obtains

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx + δ th 2 ∇ 2 R .
(2.32)

The second z-derivative in the Laplacian operator is then discarded; this yields

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx + δ th 2 ∂ 2 x R .
(2.33)

Eq. (2.33) is an augmented NPE model that includes dissipation effects. Simulations from this model have been compared to experimental results and excellent agreement was found with measurements by [START_REF] Moffett | Large-amplitude pulse propagation -A transient effect[END_REF][START_REF] Moffett | Large-amplitude pulse propagation -A transient effect[END_REF] and [START_REF] Lockwood | Directive harmonic generation in the radiation field of the circular piston[END_REF]. It has also been proved that the presented model was able to recreate self demodulation phenomena [see Too & Lee, 1995a].

Variants of the original Nonlinear Parabolic Equation model

The NPE model of Section 2.1 has been derived for a two-dimensional Cartesian coordinate system and uses the first-order "small-angle" parabolic approximation. The following two sections present Nonlinear Parabolic Equation (NPE) models for cylindrical and spherical coordinate systems (see Section 2.3.1) and a "high-angle" formulation of the Cartesian coordinate system version of the NPE (see Section 2.3.2).

Cylindrical and spherical formulations

Depending on the type of the coordinate system chosen different formulations of the NPE model are obtained. Too & Ginsberg [1992a,b] derived such formulations for cylindrical (Eq. (2.34)) and spherical (Eq. (2.35)) coordinate systems. These formulations are defined by

D t R = - c 0 2r -∂ r c 0 β 2 R 2 - c 0 2 ∂ 2 z R dr , (2.34)
and

D t Γ = -∂ r c 0 β 2 Γ 2 r - c 0 2 ∂ 2 θ + 1 tan θ ∂ θ Γ 2 r dr , (2.35)
in which r is the distance from the source to the calculation point, θ is the angle ∠ ( -→ x , -→ r )

and Γ = rR. Note that these formulations do not include sound speed variations. Further in this document the Cartesian coordinate system formulation of the NPE model defined by

Eq. (2.19) is used.

High-angle formulation

A high-angle formulation of the original NPE has been derived by [START_REF] Mcdonald | High-angle formulation for the nonlinear progressive wave equation model[END_REF] [see also [START_REF] Claerbout | Fundamentals of geophysical data processing[END_REF]. This new formulation allows one to get a better accuracy at high propagation angles, which is made possible by keeping terms of higher order in the expression of the Laplacian operator and by retaining the second time derivative term in Eq. (2.14).

Note that this development to higher orders only concerns the diffraction operator, not the nonlinear operator.

For azimuthally symmetric propagation problems the high-angle formulation of the NPE is defined by

D t R = - c 0 2 R r -∂ r c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R + c 0 R r 2 - 1 2c 0 D 2 t R dr = 0 .
(2.36)

According to McDonald the above formulation allows the Root Mean Square (RMS) error at high propagation angles (∼50°) to be reduced from 12% to 1.5% . This improved accuracy comes with increased calculation times, roughly a factor two according to McDonald. Calculating the solution array at time n + 1 using Eq. (2.36) implies working with an additional array at time n -1 due to the second time derivative. This considerably increases memory requirement, calculation complexity, and thus computational times.

Since in this work we are mainly interested in long-range sound propagation applications, the high-angle formulation has not been retained. A higher accuracy close to the source is a weak improvement compared to the numerical effort needed to integrate the additional terms.

Nonlinear Parabolic Equation model for propagation in multiple media

To the author knowledge, the only attempt to use the NPE model to simulate high-amplitude wave propagation in multiple media has been done by [START_REF] Ambrosiano | Nonlinear propagation in an ocean acoustic waveguide[END_REF], and later used by [START_REF] Mcdonald | Weak shock interaction with a free-slip interface at low grazing angles[END_REF]. Ambrosiano et al. used the NPE to investigate propagation of acoustic pulses in a shallow ocean waveguide. The ocean-bottom interface is described by a fluid-fluid interface which presents a density discontinuity. Boundary conditions are continuity of pressure and normal particle displacement across the interface.

Before going further with interfacial conditions one must address the problem of possible interface deformation by the wave. Indeed, high-amplitude waves may distort the interface in such a way that for a flat and horizontal interface, normal particle displacement is not vertical. Ambrosiano et al. investigated this issue; their conclusions are summarized below.

The interface deformation slope ζ x is to the first order

ζ x = -∂ x dt dt ′ ρ -1 ∂ z p ′ (2.37) = -∂ z R dx + O ǫ 7 2
.

(2.38)

The deformation of the surface normal introduces a nonlinear term in the normal velocity continuity condition across the interface, which is written as that if the bottom density was ρ 0 b = 1.5 ρ 0w and the overpressure was less than 0.04, the nonlinear term in Eq. (2.39) was less than 1% of the linear term. As the importance of this term is weak and decreases with range it was discarded in the expression of interfacial boundary conditions. Linear boundary conditions for a fluid-fluid interface were hence used to investigate the reflection and transmission of a high-amplitude waves in a shallow ocean waveguide.

ρ -1 0w ∂ z R w = ρ -1 0 b ∂ z R b -ρ -1 0w -ρ -1 0 b ∂ x R b ∂ z R b dx + O ǫ 7 2 , ( 2
In the case of wave propagation in air over a ground layer one has

ρ 0 ground ≫ ρ 0 air , (2.40)
such that Eq. (2.39) can be reduced to

ρ -1 0 air ∂ z R air = ρ -1 0 ground ∂ z R ground -ρ -1 0 air ∂ x R ground ∂ z R ground dx + O ǫ 7 2
(2.41)

Fortunately the acoustic field in the ground layer is very small, ensuring that the nonlinear term is small compared to the linear one in Eq. (2.41).

Using first-order finite-difference approximations of overdensity variables and their derivatives

Coupling the NPE model with near-field and far-field propagation methods

gives an expression for the interfacial boundary condition, which is defined by

R air j + R air j+1 = R ground j + R ground j+1 (pressure continuity) ,
(2.42a)

ρ -1 0 ground R ground j+1 -R ground j = ρ -1 0 air R air j+1 -R air j (normal velocity continuity) , (2.42b)
where the interface is taken to be midway between vertical grid points with indexes j + 1 and j. Eqs. (2.42) directly give expressions for the unknowns R w j and R b j+1 : The Euler equations implementation uses a rectangular grid and axisymmetrical coordinates.

R w j = ρ 0 b -ρ 0w ρ 0 b + ρ 0w R w j+1 + 2ρ 0w ρ 0 b + ρ 0w R b j , (2.43a) R b j+1 = ρ 0w -ρ 0 b ρ 0 b + ρ 0w R b j + 2ρ 0 b ρ 0 b + ρ 0w R w j+1 . ( 2
Outputs from this code are arrays of pressure and particle velocities on the grid at each time iteration. The NPE uses the same kind of data representation. The pressure values from the Euler equations code can thus be used to initialize the NPE calculations. When the pressure amplitude is low enough so that nonlinear effects can be neglected (typically less than 1 kPa), a linear propagation method can be used. The pressure values are then recorded along a vertical line and their Fourier transforms are used as the input for the frequency-domain PE.

Coupling several codes allows us to take advantage of each method strengths and weaknesses. For the case where the NPE is initialized by another propagation code, we require this code to be able to provide spatial waveforms at a given time. This can become tricky with frequencydomain methods since a number of inverse Fourier transforms and temporal domain to spatial domain mappings have to be performed. On the opposite coupling another time-domain method is straightforward. One has to take a spatial sampling of the acoustic field at a given time, extract the part of this sampling belonging to the moving-window (windowing), and start the NPE calculations with this initial condition.

The NPE model can be used to initialize any other propagation code that only requires overpressure values. Indeed, the NPE calculation result is the acoustic field at each time iteration. Time signals can easily be extracted from them, and eventually Fourier transformed to initialize a frequency-domain method like the PE. For the case where the NPE is used to initialize another time-domain method which requires velocity values (e.g. a Linearized Euler Equations (LEE) implementation), an approximation of the velocities can be obtained by using a first-order plane wave result (see Eq. (2.10)).

Significant applications using the Nonlinear Parabolic Equation model

In this section publications where the NPE model has been used are briefly reviewed. The model has mainly been used in two application fields: underwater acoustics and environmental acoustics.

Underwater acoustics

Setting the problem

Due to the large distances considered in this work the propagation model used must account for the effect of hilly or mountainous terrain on sound propagation. To achieve this goal a convenient method is the use of terrain-following coordinates. The consideration of topography in outdoor fluid dynamics by a transformed terrain-following coordinate system goes back to the 50's and the early mesoscale meteorological models [see for example [START_REF] Philips | A coordinate system having some special advantages for numerical forecasting[END_REF][START_REF] Mahrer | The numerical study of the airflow over mountains using the University of Virginia mesoscale model[END_REF].

The ground height is noted h (x) and its first and second derivatives with respect to x (the main propagation direction) are noted h ′ and h ′′ , respectively. Modeling with transformed coordinates often involves the transformation proposed by Gal-Chen & Sommerville [1975] in which the spatial coordinates are transformed according to While pressure, density, sound speed and horizontal flow velocity variables are unchanged, the vertical flow velocity is modified to account for ground elevation, such that

x -→ x , (3.1a) z -→ H z -h (x) H -h (x) , ( 3
w -→ w + h ′ u . (3.3)
Expressions for spatial derivatives of space and time dependent variables can be found using the chain rule. One obtains

∂ x R -→ ∂ x R -h ′ ∂ z R , (3.4a) ∂ z R -→ ∂ z R , (3.4b) ∂ x u -→ ∂ x u -h ′ ∂ z u , (3.4c) ∂ z u -→ ∂ z u , (3.4d) ∂ x w -→ ∂ x w -h ′ ∂ z w + h ′ ∂ x u -h ′ ∂ z u + h ′′ u , (3.4e) ∂ z w -→ ∂ z w + h ′ ∂ z u . (3.4f)

Model derivation

In order to derive a NPE model for propagation over non-flat surfaces we start from the nonlinear wave equation Eq. (2.11) written in dimensionless form, such that

∂ 2 t R = ∂ 2 x c 2 0 + 2c 0 c 1 R + c 2 0 βR 2 + c 2 0 ∂ 2 z R . (3.5)
Performing the replacements in Eqs. (3.4) yields a nonlinear wave equation in the transformed coordinate system. This equation writes

∂ 2 t R = D 2 x c 2 0 + 2c 0 c 1 R + c 2 0 βR 2 + c 2 0 ∂ 2 z R , (3.6)
where D 2 x is the expression of the second x-derivative in the transformed coordinate system, which is defined by

D 2 x -→ ∂ 2 x + h ′2 ∂ 2 z -h ′′ ∂ z -2h ′ ∂ x ∂ z (3.7a) -→ ∂ 2 x + L , (3.7b) where L = h ′2 ∂ 2 z -h ′′ ∂ z -2h ′ ∂ x ∂ z .
Introducing the moving-frame operator as in Section 2.1.2 yields

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx -L c 0 2 + c 1 R + c 0 β 2 R 2 dx . (3.8)
The assumptions inherent to the NPE model are now used. Specifically, small meteorological perturbations and weak nonlinearities are introduced, i.e.

c 1 ≪ c 0 2 , (3.9a) βR ≪ 1 . (3.9b)
Associated with the assumption of dominant propagation in one direction * , i.e.

∂ x R ≫ ∂ z R,
Eqs. (3.9) allows Eq. (3.8) to be reduced to

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx - c 0 2 LR dx . (3.10)
Rearranging Eq. (3.10) gives the NPE model with terrain-following coordinates defined as for a graphical representation of additional terms):

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 1 + h ′2 ∂ 2 z R + h ′′ ∂ z R dx + c 0 h ′ ∂ z R . ( 3
◮ the term 1 + h ′2 ∂ 2 z R in the integral increases transverse propagation speed by a quantity c 0 h ′2 , independently of the slope sign, ◮ the second term h ′′ ∂ z R in the integral increases or decreases the transverse propagation speed depending on the ground slope direction,

◮ and the last term c 0 h ′ ∂ z R shifts the waveform in the calculation window to correct for ground elevation.

Including thermoviscous effects in the Generalized Terrain -Nonlinear Parabolic Equation model

To include thermoviscous effects in the GT-NPE model we start from Eq. (2.33) where the dissipative term is written in the transformed coordinate system, such that

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 1 + h ′2 ∂ 2 z R + h ′′ ∂ z R dx + c 0 h ′ ∂ z R + δ th 2 ∂ 2 x + 1 + h ′2 ∂ 2 z -h ′′ ∂ z -2h ′ ∂ x ∂ z R .
(3.12)

Figure 3.1: Illustration of the coordinate transformation used to take into account the ground topography within the NPE model. Following the derivation of Section 2.1.2 in which one keeps only dominant terms allows one to get a GT-NPE model with thermoviscous effects included, defined by

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 1 + h ′2 ∂ 2 z R + h ′′ ∂ z R dx + c 0 h ′ ∂ z R ◮ Locally
D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 1 + h ′2 ∂ 2 z R + h ′′ ∂ z R dx + c 0 h ′ ∂ z R + δ th 2 ∂ 2 x R .
(3.13)

Model properties

The GT-NPE model developed can be used to simulate sound propagation over non-flat terrains. However, steepest slopes should not exceed a 45°angle, i.e. max |∂ x h| < 1. Indeed larger slopes may yield numerical instabilities unless spatial and temporal steps are adapted accordingly. Generally speaking, the stability of finite-difference schemes depends on the Courant-Friedrichs-Lewy (CFL) number defined as

CFL = ∆l ∆t c max , (3.14)
where ∆l is the smallest spatial step, i.e. ∆l = min (∆x, ∆z), in which ∆x and ∆z are the spatial steps in the x-and z-directions, respectively, ∆t is the time step and c max is the maximum wave speed. For the finite difference scheme to be stable the condition CFL < c must be satisfied, where c is some positive constant that depends on the numerical scheme used. In the transformed coordinate system the smallest spatial steps ∆s (x) depends on the terrain slope so that it decreases with increasing slope such that

min (∆s (x)) = ∆l 1 + max |∂ x h| 2 . (3.15)
For example, assuming a maximum slope max |∂ x h| = 5 implies decreasing the time step ∆t by a factor 5.1, and thus increasing the calculation time for this time iteration by about the same factor. A fortiori obstacles with sharp angles and complex shapes cannot be modeled by the terrain-following coordinates method. This would yield infinite slopes and hence infinitely small time steps.

Besides these numerical considerations, a second limitation on the slope of the ground topography arises from the approximations made during the derivation. It can be seen in Eq. (3.8)

that the integral contains a meteorological perturbation term (c 1 term) and a nonlinear effect term (R 2 term). As explained in Section 2.1.3 one interesting feature of the NPE model is its modularity: various effects can be included or excluded from the model and are segregated into separate differential operators. Keeping these additional terms under the integral would destroy this property since one would have to take into account meteorological and nonlinear effects in both directions * . Neglecting these terms implies that ◮ the propagation is mainly along the transformed x-direction, ◮ and that the ground topography derivatives are small, i.e. h ′ ≪ 1 and h ′′ ≪ 1.

Ground topography is hence restricted to gentle slopes. This ensures numerical stability and model correctness, while keeping all features of the NPE model and without any major modification concerning numerical implementation. [1987]). We use the same procedure to prove that the GT-NPE Eq. (3.11) is a proper equivalence of the Generalized Terrain -Parabolic Equation (GT-PE) proposed by [START_REF] Sack | A parabolic equation for sound propagation in two dimensions over any smooth terrain profile: the generalized terrain parabolic equation (GT-PE)[END_REF]. We drop the nonlinear term in Eq. (3.11)

and substitutes -ωt) .

R = f (x, z) e i(kx
(3.16)

Neglecting non-dominant terms gives the first-order "narrow-angle" GT-PE defined by

∂ x f = -ik c 1 c 0 f + i 2k 1 + h ′2 ∂ 2 z f -2ikh ′ + h ′′ ∂ z f .
(3.17)

The GT-NPE model derived in Section 3.1.2 is thus the nonlinear time-domain counterpart of the GT-PE.

Nonlinear Parabolic Equation model for sound propagation in rigidly-framed porous media

The problem with time-domain sound propagation models is that an absorbing ground surface cannot be taken into account by an impedance condition. Complex impedance models can be used with frequency-domain models, for example the PE or the Fast Field Program (FFP) models, but cannot directly be used with the NPE model developed. The ground layer has hence to be included into the computational system as a propagation medium.

In this section it is proposed to derive a Nonlinear Parabolic Equation (NPE) model for porous layers similar to Eq. (2.19). This model uses a minimal parametrization in which the layer is assumed to be equivalent to a continuous fluid medium: a wave causes a vibration of air particles contained in the ground pores, while the ground frame does not vibrate.

Model derivation

In this section a Nonlinear Parabolic Equation (NPE) model for high amplitude wave propagation within rigidly-framed porous media is developed. The Nonlinear Parabolic Equation (NPE) model for sound propagation in porous ground media is based on a nonlinear extension of the Zwikker-Kosten (ZK) model (see [START_REF] Zwikker | Sound absorbing materials[END_REF] and Section 1.1.3). It is characterized by a set of 4 parameters:

◮ the static air flow resistivity σ 0 ,

◮ the porosity Ω 0 , ◮ the tortuosity Φ,

◮ and the Forchheimer nonlinearity parameter ξ.

The flow resistivity σ 0 characterizes the visco-inertial effects at low frequencies and the porosity Ω 0 is defined as the ratio of the fluid volume occupied by the continuous fluid phase to the total volume of the porous material. The tortuosity Φ is defined as the ratio of a curved path length to the distance between its end points and the Forchheimer nonlinearity parameter ξ characterizes the nonlinear part of the flow resistivity (see Section 1.2.3). These quantities are assumed constant within the ground layer and with respect to time.

Considering these assumptions equations of continuity and conservation of momentum are [see for example Krylov et al. , 1996a,b;Védy, 2002b;[START_REF] Umnova | High amplitude pulse propagation and reflection from a rigid porous layer[END_REF]]

∂ t ρ T = -∂ x (ρ T u) -∂ z (ρ T w) , (3.18a) Φ∂ t (ρ T u) + ∂ x p T + Φρ T u 2 + ∂ z (Φρ T uw) = -σ 0 Ω 0 (1 + ξ |u|) u , (3.18b) Φ∂ t (ρ T w) + ∂ z p T + Φρ T w 2 + ∂ x (Φρ T uw) = -σ 0 Ω 0 (1 + ξ |w|) w . (3.18c)
As one can see in Eqs. (3.18), the tortuosity Φ reduces the pressure gradients and flow resistive terms amplitude.

In order to obtain a single variable wave equation, Eq. (3.18a) is derived with respect to time and multiplied by the tortuosity Φ. This yields

Φ∂ 2 t ρ T + ∂ x [Φ∂ t (ρ T u)] + ∂ z [Φ∂ t (ρ T w)] = 0 . (3.19)
Eqs. (3.18b, 3.18c) are then used to replace the terms Φ∂ t (ρ T u) and Φ∂ t (ρ T w), so that Eq. (3.19) becomes

Φ∂ 2 t ρ T = ∂ 2 x p T + Φρ T u 2 + ∂ 2 z p T + Φρ T w 2 + 2Φ∂ x ∂ z (ρ T uw) + σ 0 Ω 0 ∂ x [(1 + ξ |u|) u] + σ 0 Ω 0 ∂ z [(1 + ξ |w|) w] .
(3.20)

Since the propagation is mainly along the x-axis, only linear terms in z-derivatives are kept in Eq. (3.20): terms ∂ x ∂ z (Φρ T uw), ∂ 2 z Φρ T w 2 and σ 0 Ω 0 ∂ z (ξ |w| w) are neglected. Moreover, only terms of order up to two in x-derivatives are retained: the quantity ∂ 2

x Φρ ′ u 2 is discarded. This leads us to

Φ∂ 2 t ρ T = ∂ 2 x p T + Φρ 0 u 2 + ∂ 2 z p T + σ 0 Ω 0 ∂ x [(1 + ξ |u|) u] + σ 0 Ω 0 ∂ z w .
(3.21)

One must now find expressions for the flow velocities u and w. The perturbation expansion method is used, in which the same scalings and expansions as in publications by [START_REF] Mcdonald | Time domain formulation for pulse propagation including nonlinear behaviour at a caustic[END_REF] and [START_REF] Mcdonald | A tutorial on the nonlinear progressive wave equation (NPE) -Part 1[END_REF] are used (however note that the moving frame speed is set to c 0 / √ Φ). The replacements are given by

x -→ x - c 0 √ Φ t , (3.22a) z -→ ǫ 1/2 z , (3.22b) t -→ ǫt , (3.22c)
in which ǫ is a scaling factor. The scaling of z by a factor of ǫ 1/2 emphasizes the predominance of the propagation in the x-direction. The partial derivatives associated with Eqs. (3.22) are

∂ x -→ ∂ x , (3.23a) ∂ z -→ ǫ 1/2 ∂ z , (3.23b) ∂ t -→ ǫ∂ t - c 0 √ Φ ∂ x . (3.23c)
The dependent variables are expanded as follows:

ρ -→ ρ 0 + ǫρ 1 + ǫ 2 ρ 2 + • • • , (3.24a) u -→ ǫu 1 + ǫ 3/2 u 2 + • • • , (3.24b) w -→ ǫw 1 + ǫ 3/2 w 2 + • • • . ( 3
ǫ∂ t - c 0 √ Φ ∂ x ρ 0 + ǫρ ′ 1 + ǫ 2 ρ ′ 2 + • • • = -∂ x ρ 0 + ǫρ 1 + ǫ 2 ρ 2 + • • • ǫu 1 + ǫ 3/2 u 2 + • • • -ǫ 1/2 ∂ z ρ 0 + ǫρ 1 + ǫ 2 ρ 2 + • • • ǫw 1 + ǫ 3/2 w 2 + • • • . (3.25)
Elimination of particle flow velocities u and w in Eq. (3.25) is performed introducing a first order approximation of these quantities. It can be obtained by equalizing terms of order ǫ and ǫ 3/2 in Eq. (3.25). One obtains

u 1 = c 0 √ Φ ρ 1 ρ 0 , (3.26a
)

w 1 = 0 . (3.26b)
Furthermore, note that

ρ ′ = ρ 1 + O ǫ 2 , (3.27a) u = u 1 + O ǫ 3/2 , (3.27b) w = w 1 + O ǫ 3/2 . (3.27c)
Substitution of u and w by u 1 and w 1 in Eq. (3.21), respectively, hence leads us to an error consistent with the assumptions made. With these replacements Eq. (3.21) becomes

Φ∂ 2 t ρ T = ∂ 2 x p T + c 2 0 ρ 0 ρ ′2 + ∂ 2 z p T + σ 0 Ω 0 c 0 ρ 0 √ Φ ∂ x 1 + ξc 0 √ Φ ρ ′ ρ 0 ρ ′ . (3.28)
The total pressure p T is then eliminated from Eq. (3.28) by using a second-order expansion in ρ ′ from an assumed adiabatic equation of state. Variable p T is hence substituted by

p T = p 0 + c 2 0 ρ ′ + c 2 0 γ -1 2ρ 0 ρ ′2 . (3.29)
Substituting Eq. (3.29) in Eq. (3.28) yields a single-variable nonlinear wave equation for propagation in porous media defined by

Φ∂ 2 t ρ ′ = c 2 0 ∂ 2 x ρ ′ + γ + 1 2ρ 0 ρ ′2 + c 2 0 ∂ 2 z ρ ′ + σ 0 Ω 0 c 0 ρ 0 √ Φ ∂ x 1 + ξc 0 √ Φ ρ ′ ρ 0 ρ ′ . (3.30)
The one-way propagation hypothesis of NPE models is introduced with a moving-frame operator D ⋆ t , defined by

D ⋆ t = ∂ t + c 0 √ Φ ∂ x .
(3.31)

The parabolic approximation (see Eqs. (2.15, 2.16)) yields

∂ t = D ⋆ t - c 0 √ Φ ∂ x , (3.32a) ∂ 2 t = -2 c 0 √ Φ D ⋆ t ∂ x + c 2 0 Φ ∂ 2 x . (3.32b)
Replacing the second time derivative in Eq. (3.30) by the expression in Eq. (3.32b) gives a first version of the NPE model for propagation in porous media, defined by

-2c 0 √ ΦD ⋆ t ∂ x ρ ′ = c 2 0 ∂ 2 x γ + 1 2ρ 0 ρ ′2 + c 2 0 ∂ 2 z ρ ′ + c 0 σ 0 Ω 0 ρ 0 √ Φ ∂ x 1 + ξc 0 √ Φ ρ ′ ρ 0 ρ ′ . (3.33)
Integrating Eq. (3.33) with respect to x, rearranging, and introducing the dimensionless overdensity variable R gives a NPE model for propagation in porous media, such that

D ⋆ t R = - c 0 √ Φ ∂ x β 2 R 2 - c 0 2 √ Φ ∂ 2 z R dx - σ 0 Ω 0 2Φρ 0 1 + ξc 0 √ Φ |R| R . (3.34)
Eq. (3.34) can be used to simulate sound propagation in a porous ground layer. However, if one wants to use this NPE model together with the NPE model for atmospheric propagation, a last operation must be done. Indeed, both models use different frame speed: c 0 and c 0 / √ Φ for atmospheric and porous medium propagation models, respectively. A last transformation is introduced to account for the moving-frame speed difference defined by c 0 (1 -1/ √ Φ ). The operator D ⋆ t in Eq. (3.34) is replaced by

D ⋆ t = D t -c 0 1 - 1 √ Φ ∂ x , (3.35)
so that Eq. (3.34) becomes

D t R = - c 0 √ Φ ∂ x 1 - √ Φ R + β 2 R 2 - c 0 2 √ Φ ∂ 2 z R dx - σ 0 Ω 0 2Φρ 0 1 + ξc 0 √ Φ |R| R . (3.36)
The NPE model described by Eq. (3.36) is able to simulate finite amplitude sound propagation within a rigidly-framed porous material described by a set of 4 parameters. Note that if one sets Φ= 1 and neglects losses in the layer, i.e. σ 0 = 0, the model exactly reduces to the usual NPE for atmospheric propagation defined by Eq. (2.19).

Terrain-following coordinates formulation.

Transforming Eq. (3.36) in the terrainfollowing coordinates system described in Section 3.1.1 is straightforward. The same assumptions as in Section 3.1.2 are made:

◮ small meteorological perturbations,

◮ and weak nonlinearities.

This allows us to get a NPE model for porous grounds with terrain-following coordinates, defined by

D t R = - c 0 √ Φ ∂ x 1 - √ Φ R + β 2 R 2 - σ 0 Ω 0 2Φρ 0 1 + ξc 0 √ Φ |R| R - c 0 2 √ Φ 1 + h ′2 ∂ 2 z R + h ′′ ∂ z R dx + c 0 √ Φ h ′ ∂ z R . (3.37)
One more time if one sets Φ= 1 and neglects losses in the layer, i.e. σ 0 = 0, Eq. (3.37) exactly reduces to the GT-NPE model for atmospheric propagation over non-flat terrains presented in Section 3.1.2 (Eq. (3.11)).

Model properties

Eq. (3.36) allows some conclusions about finite-amplitude sound propagation in porous media to be drawn (see Fig. ◮ the sound speed in the medium is inversely proportional to the square root of the material tortuosity, i.e.

c ground = c 0 √ Φ , (3.38)
◮ the attenuation in the ground layer is composed of a linear term plus a nonlinear term,

◮ and with the hypothesis used (the plane wave approximation for the flow velocity component u, see Eq. (3.27b)), the dynamic material flow resistivity is linearly proportional to the overdensity absolute value |R|, such that

σ = σ 0 1 + ξ c 0 √ Φ |R| .
(3.39)

Nonlinear Parabolic Equation model for high-amplitude wave propagation over complex surfaces

In the previous sections two different NPE models have been derived: the first can be used for atmospheric propagation while the second is for finite-amplitude sound propagation within rigidly-framed porous media. Both can eventually be adapted to handle non-flat topographies.

D t R = - c 0 √ Φ ∂ x 1 - √ Φ R + β 2 R 2 - c 0 2 √ Φ ∂ 2 z R dx - σ 0 Ω 0 2Φρ 0 1 + ξ c 0 √ Φ |R| R ◮ Correction
As both models use the same moving-frame speed, they can be used together to simulate finiteamplitude sound propagation over a rigidly-framed porous ground layer. This section aims at establishing a first-order boundary interface condition to link these two Nonlinear Parabolic Equation (NPE) models. It will be then transformed to handle non-flat topography cases.

In the following we assume that the deformation of the interface by the wave is small (see [START_REF] Ambrosiano | Nonlinear propagation in an ocean acoustic waveguide[END_REF] and Section 2.4).

Derivation of the boundary interface condition

An air layer, whose fields are noted p ′a , u a and w a , is considered. To construct the air-ground interface condition a rigidly-framed porous ground layer is introduced; its fields are noted p ′g , u g and w g . With these notations boundary conditions across the interface are continuity of pressure and normal flow velocity, defined by

p ′a = p ′g , (3.40a) [w a ] = [w g ] , (3.40b)
where the square brackets denote the field quantity on the air-ground interface.

Expressions of w a and w g involving the pressure disturbance p ′ to the first order are sought; linearized equations are hence used. For the air layer we use the linearized Euler's equation

ρ 0 ∂ t (w a ) = -∂ z p a T .
(3.41)

The perturbation expansion method is used and the same scalings as in Section 3.2. Rewriting Eq. (3.41) with these scalings and expansions yields

ρ 0 (ǫD t -c 0 ∂ x ) ǫw a 1 + ǫ 3/2 w a 2 + . . . = -ǫ 1/2 ∂ z p 0 + ǫp ′ a 1 + . . . , (3.42)
and by equalizing terms of order 1 and 3/2 one can find

w a 1 = 0 , (3.43a) w a 2 = (ρ 0 c 0 ∂ x ) -1 ∂ z p ′ a 1 . (3.43b)
Note that

w a = w a 1 + w a 2 + O ǫ 5/2 , (3.44)
and hence, to the order of accuracy sought in this work one can write

w a = w a 1 + w a 2 , (3.45a
)

w a = (ρ 0 c 0 ∂ x ) -1 ∂ z p ′ a 1 . (3.45b)
Now that an expression for the vertical flow velocity in the air layer w a has been obtained one has to repeat the procedure for w g , the vertical flow velocity in the ground layer. To find an expression for w g we start from the following momentum equation * :

Φρ 0 ∂ t w g = -Ω 0 ∂ z p g T -σ 0 Ω 0 w g . (3.46)
The same procedure is applied; Eq. (3.46) becomes

Φρ 0 (ǫD t -c 0 ∂ x ) ǫw g 1 + ǫ 3/2 w g 2 + . . . = -ǫ 1/2 Ω 0 ∂ z p 0 + ǫp ′ g 1 + . . . -σ 0 Ω 0 ǫw g 1 + ǫ 3/2 w g 2 + . . . , (3.47)
and one hence obtains

√ Φρ 0 c 0 ∂ x -σ 0 Ω 0 w g 1 = 0 , (3.48a) √ Φρ 0 c 0 ∂ x -σ 0 Ω 0 w g 2 = Ω 0 ∂ z p ′ g 1 . (3.48b)
With the help of Eqs. (3.48) one can finally obtain an expression for the flow velocity in the * Eq. (3.46) is given as equation 10 in [START_REF] Salomons | Eulerian time-domain model for sound propagation over a finite impedance ground surface. Comparison with frequency-domain models[END_REF] ground layer w g involving the first-order pressure disturbance p ′ g 1 , such that

w g = √ Φρ 0 c 0 ∂ x -σ 0 Ω 0 -1 Ω 0 ∂ z p ′ g 1 . (3.49)
With the help of Eqs. (3.45b, 3.49) the boundary condition across the interface for the continuity of vertical velocities w a and w g can be written. It is defined by

(ρ 0 c 0 ∂ x ) -1 ∂ z p ′ a = √ Φρ 0 c 0 ∂ x -σ 0 Ω 0 -1 Ω 0 ∂ z p ′ g . (3.50)
Rearranging Eq. (3.50) and re-introducing the pressure continuity condition leads us to

p ′a = p ′g , (3.51a) √ Φ∂ z p ′ a - σ 0 Ω 0 ρ 0 c 0 ∂ z p ′ a dx = Ω 0 ∂ z p ′ g . (3.51b)
Eqs. (3.51) is the boundary interface condition that couples the atmospheric and porous ground media NPE models. As it can be seen the coupling equations only involve spatial derivatives and integrals instead of the time integral that is usually found in such time-domain impedance conditions [see for example [START_REF] Wilson | Time-domain calculations of sound interactions with outdoor ground surfaces[END_REF][START_REF] Wilson | Time-domain modeling of the acoustic impedance of porous surfaces[END_REF][START_REF] Wilson | Time-domain equations for sound propagation in rigid-frame porous media[END_REF]. This integral over time comes from the convolution that appears when the relation between flow velocity and pressure is transformed from the frequency-domain to the time-domain (see Section 1.1.3).

From a numerical point of view working with an integral over time is often synonym of fastidious implementation and long calculation times. The use of the frame-following formulation inherent to NPE models allows us to separate this time integral in two parts: a time and a spatial integral, where we finally retain only the dominant part, the spatial integral.

Discretization of the boundary interface condition

In this section the continuous form of the boundary interface condition shown in Eqs. (3.51) is discretised using the finite-difference method. The variables p ′ a i,j and p ′ g i,j are introduced to denote pressures in layer a (air layer) and layer g (porous ground layer), respectively, at range i∆x in the moving window and altitude j∆z. The air-ground interface is taken to be midway between two vertical grid points with indexes j = 0 and j = 1; auxiliary virtual points with pressures p used to discretise Eq. (3.51). For a generic layer l we use

p ′ l = p ′ l i,1 + p ′ l i,0 2 , (3.52a)
∂ z p ′ l = p ′ l i,1 -p ′ l i,0 ∆z -1 . (3.52b)
Replacing these approximations into Eqs. (3.51) gives expressions for unknown quantities p ′ a i,0

and p

′ g i,1
, such that

(A + G) p ′ a i,0 = (A -G) p ′ a i,1 + 2Gp ′ g i,0 + S i+1 m=Nx p ′ a m,1 -p ′ a m,0 , (3.53a)
(A + G) p ′ g i,1 = (G -A) p ′ g i,0 + 2Ap ′ a i,1 + S i+1 m=Nx p ′ a m,1 -p ′ a m,0 , (3.53b)
where N x is the number of points in the moving window in the x-direction and

A = √ Φ + 1 2 S , (3.54a
)

G = Ω 0 , (3.54b) S = σ 0 Ω 0 ∆x c 0 ρ 0 . (3.54c)
Eqs. (3.53) associated with Eqs. (3.54) give expressions for the unknown pressures p a i,0 and p g i,1 . Used together with the atmospheric and porous ground NPE models, these expressions allow weakly nonlinear sound propagation over a finite-impedance ground surface to be simulated.

Model properties

Limitations

First-order formulations of the constitutive equations have been used to derive the boundary interface condition. This implies that nonlinearities cannot be taken into account in the two-way coupling.

Causality

The x-integral present in NPE models (see for example Eq. (2.19)) is calculated from the right side of the calculation grid to the left side, and the same method is used for coupling (note the reversed sum indexes in Eqs. (3.53)). This ensures that no perturbation is introduced ahead of the point where the wave hits the ground, and thus implies that the interfacial condition is causal.

Consistency to simple boundary conditions

Classical boundary conditions can be obtained by setting specific values to the quantities Φ, σ 0 and Ω 0 .

◮ If one sets Φ= +∞ one obtains

R a i,0 = R a i,1 , (3.55)
which is the condition, with the discretization used, for an acoustically rigid surface (null pressure z-derivative).

◮ If one sets σ 0 = 0, Ω 0 = 1 and Φ= 1 (parameters for an air layer) then

A = 1 , (3.56a) G = 1 , (3.56b) S = 0 , (3.56c)
and hence

R a i,0 = R g i,0 , (3.57a) R g i,1 = R a i,1 , (3.57b)
which is the condition for a transparent interface (perfect transmission).

◮ If one sets σ 0 = 0 and Ω 0 = 1, Eqs. (3.53) become

R a i,0 = √ Φ -1 √ Φ + 1 R a i,1 + 2 √ Φ + 1 R g i,0 , (3.58a) R g i,1 = 1 - √ Φ √ Φ + 1 R g i,0 + 2 √ Φ √ Φ + 1 R a i,1 , (3.58b)
which is the interface condition for two fluid layers with densities ρ 0 and √ Φρ 0 [see [START_REF] Ambrosiano | Nonlinear propagation in an ocean acoustic waveguide[END_REF].

Relation to frequency-domain impedance

The frequency-domain complex characteristic impedance corresponding to the boundary conditions shown in Eqs. (3.51) is given by (see Section 1.1.3 and [START_REF] Salomons | Eulerian time-domain model for sound propagation over a finite impedance ground surface. Comparison with frequency-domain models[END_REF])

Z c = ρ 0 c 0 Φ Ω 2 0 + i σ 0 Ω 0 ρ 0 ω . (3.59)
On the causality of the Zwikker-Kosten (ZK) model. The causality of the Zwikker-Kosten (ZK) model can be studied with the causality index proposed by [START_REF] Berthelot | Surface acoustic impedance and causality[END_REF]. In its nonlocal form the causality index C nl is defined by

C nl = HT [Re (Z c )] Im (Z c ) , ( 3

Terrain-following coordinates formulation of the boundary interface condition

The boundary interface condition developed in Section 3.3.1 can be adapted to handle nonflat terrains. Expressions for the vertical flow velocity components w a and w g are modified

ρ 0 c 0 ∂ t w a -h ′ u a = -∂ z p a T , (3.61a
)

Φρ 0 ∂ t w g -h ′ u g = -Ω 0 ∂ z p g T -σ 0 Ω 0 w g -h ′ u g . (3.61b)
Using the perturbation expansion method as described in Section 3.3.1 and a plane-wave approximation for the horizontal flow velocities u a and u g allows us to obtain transformed expressions for Eqs. (3.45b, 3.49), such that

w a = (ρ 0 c 0 ∂ x ) -1 ∂ z p ′ a 1 -h ′ ∂ x p ′ a 1 , (3.62a
)

w g = √ Φρ 0 c 0 ∂ x -σ 0 Ω 0 -1 Ω 0 ∂ z p ′ g 1 -h ′ √ Φ∂ x - σ 0 Ω 0 ρ 0 c 0 p ′ g 1 . (3.62b)
The vertical flow velocity equality condition in the transformed coordinates system is hence

√ Φ ∂ z p ′ a -h ′ ∂ x p ′ a - σ 0 Ω 0 ρ 0 c 0 ∂ z p ′ a -h ′ ∂ x p ′ a dx = Ω 0 ∂ z p ′ g -h ′ √ Φ∂ x p ′ g - σ 0 Ω 0 c 0 ρ 0 p ′ g . (3.63)
The discretization with finite-difference approximations of the continuous variables gives modified expressions for Eqs. (3.53), defined by

(A 0 + G 1 ) p ′ a i,0 = (A 1 -G 1 ) p ′ a i,1 + (G 0 + G 1 ) p ′ g i,0 + S i+1 m=Nx p ′ a m,1 -p ′ a m,0 , (3.64a)
(A 0 + G 1 ) p ′ g i,1 = (G 0 -A 0 ) p ′ g i,0 + (A 0 + A 1 ) p ′ a i,1 + S i+1 m=Nx p ′ a m,1 -p ′ a m,0 , (3.64b)
in which the coupling coefficients are

A 0 = √ Φ + σ 0 Ω 0 ∆x 2c 0 ρ 0 1 - h ′′ ∆z 2 -h ′ ∆z 2 σ 0 Ω 0 ρ 0 c 0 + √ Φ ∆x , (3.65a
)

A 1 = √ Φ + σ 0 Ω 0 ∆x 2c 0 ρ 0 1 + h ′′ ∆z 2 + h ′ ∆z 2 σ 0 Ω 0 ρ 0 c 0 + √ Φ ∆x , (3.65b) G 0 = Ω 0 -h ′ ∆z 2 σ 0 Ω 0 ρ 0 c 0 + √ Φ ∆x , (3.65c
)

G 1 = Ω 0 + h ′ ∆z 2 σ 0 Ω 0 ρ 0 c 0 + √ Φ ∆x , (3.65d)
S = σ 0 Ω 0 ∆x c 0 ρ 0 .
(3.65e)

Boundary interface condition for multilayered ground surfaces

A similar boundary interface condition can be derived for multilayered ground surfaces. We consider two ground layers g and b, whose fields are noted with a g or b superscript. It is assumed than the layer b is below the layer g and that the interface may be non-flat. We get for the vertical velocity components

w g = √ Φ g ρ 0 c 0 ∂ x -σ g 0 Ω g 0 -1 Ω g 0 ∂ z p ′ g 1 -h ′ √ Φ g ∂ x - σ g 0 Ω g 0 ρ 0 c 0 p ′ g 1 , (3.66a) w b = √ Φ b ρ 0 c 0 ∂ x -σ b 0 Ω b 0 -1 Ω b 0 ∂ z p ′ b 1 -h ′ √ Φ b ∂ x - σ b 0 Ω b 0 ρ 0 c 0 p ′ b 1 . (3.66b)
Writing the equality across the interface and using discrete expressions for the pressure values

Including Forchheimer's nonlinearities in the two-way coupling

While the flow resistivity dependence on particle velocity (Forchheimer's nonlinearities) are accounted for in the NPE model for porous ground layers (last term in Eq. (3.36)), the two-way coupling between both domains does not contain high-amplitude effects on ground properties. This would yield erroneous solutions, since an additional attenuation would be introduced in the ground layer but the increased rigidity of the interface would not be accounted for.

A solution is to artificially increase the static flow resistivity σ 0 in Eq. (3.54) (and its variants Eqs. (3.65, 3.71)) according to

σ (x, t) = σ 0 1 + ξ w i , (3.67)
where w i is the vertical particle velocity at the interface. Note that the flow resistivity is now dependent on (x, t); it is thus noted σ (x, t). We then use Eq. (3.43b) to obtain an approximation of w i , such that

w i = (ρ 0 c 0 ) -1 ∂ z p ′ i 1 dx + O ǫ 5/2 , (3.68)
where p

′ i
1 is the first-order approximation of the pressure on the interface. The flow resistivity σ 0 in the coupling parameters defined in Eqs. (3.54) and its variants Eqs. (3.65, 3.71) is thus replaced by

σ (x, t) = σ 0 1 + ξ ρ 0 c 0 ∂ z p ′ i 1 dx . (3.69)
At each time iteration the flow resistivity is updated with the help of pressure values on the interface at the previous time step. This method, although approximate, allows Forchheimer's nonlinearities to be included in the two-way coupling equations.

Chapter summary and conclusions

This chapter detailed the development of a deterministic NPE model for high-amplitude wave propagation in complex media.

First, propagation over non-flat surfaces is handled through the use of terrain-following coordinates. This method yields a NPE model with additional terms to account for ground elevation and is restricted to gentle slopes (numerical and physical limitations). This NPE model will allow, once the probabilistic model of uncertainties is developed, the use of the stochastic NPE model in hilly urban areas.

R g i,0 = G 1 -B 1 G 0 + B 1 R g i,1 + B 0 + B 1 G 0 + B 1 R b i,0 + S g G 0 + B 1 i+1 m=N x S g,1 R g i,1 -S g,0 R g i,0 - S b G 0 + B 1 i+1 m=N x S b,1 R b i,1 -S b,0 R b i,0 , (3.70a) R b i,1 = B 0 -G 0 G 0 + B 1 R b i,0 + G 1 + G 0 G 0 + B 1 R G i,1 + S g G 0 + B 1 i+1 m=N x S g,1 R g i,1 -S g,0 R g i,0 - S b G 0 + B 1 i+1 m=N x S b,1 R b i,1 -S b,0 R b i,0 . (3.70b)
with: Chapter conclusions are given in Section 4.7.

G 0 = √ Φ b Ω g 0 + σ b 0 Ω b 0 ρ 0 c 0 ∆x 2 Ω g 0 - ∆z 2 h ′′ Φ g -h ′ σ g 0 Ω g 0 ρ 0 c 0 -h ′ ∆z 2 Φ b σ g 0 Ω g 0 ρ 0 c 0 + Φ g √ Φ b ∆x + Φ g σ b 0 Ω b 0 ρ 0 c 0 , (3.71a) G 1 = √ Φ b Ω g 0 + σ b 0 Ω b 0 ρ 0 c 0 ∆x 2 Ω g 0 + ∆z 2 h ′′ Φ g -h ′ σ g 0 Ω g 0 ρ 0 c 0 + h ′ ∆z 2 Φ b σ g 0 Ω g 0 ρ 0 c 0 + Φ g √ Φ b ∆x + Φ g σ b 0 Ω b 0 ρ 0 c 0 , (3.71b) B 0 = √ Φ g Ω b 0 + σ g 0 Ω g 0 ρ 0 c 0 ∆x 2 Ω b 0 - ∆z 2 h ′′ Φ b -h ′ σ b 0 Ω b 0 ρ 0 c 0 -h ′ ∆z 2 Φ g σ b 0 Ω b 0 ρ 0 c 0 + √ Φ b Φ g ∆x + Φ b σ g 0 Ω g 0 ρ 0 c 0 , (3.71c) B 1 = √ Φ g Ω b 0 + σ g 0 Ω g 0 ρ 0 c 0 ∆x 2 Ω b 0 + ∆z 2 h ′′ Φ b -h ′ σ b 0 Ω b 0 ρ 0 c 0 + h ′ ∆z 2 Φ g σ b 0 Ω b 0 ρ 0 c 0 + √ Φ b Φ g ∆x + Φ b σ g 0 Ω g 0 ρ 0 c 0 , (3.71d) S g,1 = σ b 0 Ω b 0 ∆x ρ 0 c 0   Ω g + ∆z 2   h ′′ Φ g -h ′ σ g 0 Ω g 0 ρ 0 c 0     and S g,0 = σ b 0 Ω b 0 ∆x ρ 0 c 0   Ω g - ∆z 2   h ′′ Φ g -h ′ σ g 0 Ω g 0 ρ 0 c 0     ,(3.71e) S b,1 = σ g 0 Ω g 0 ∆x ρ 0 c 0   Ω b + ∆z 2   h ′′ Φ b -h ′ σ b 0 Ω b 0 ρ 0 c 0     and S b,0 = σ g 0 Ω g 0 ∆x ρ 0 c 0   Ω b - ∆z 2   h ′′ Φ b -h ′ σ b 0 Ω b 0 ρ 0 c 0     . (3.71f)

Chapter content

The operator splitting method

An operator splitting method is used to solve the various NPEs. The operator splitting method consist of a "divide and conquer strategy", where the original Partial Differential Equation (PDE) is split into simpler sub-problems which are treated individually with specialized numerical algorithms. The NPEs, see for example Eq. (2.19) in Section 2.1.2, are of the form

D t R = LR , (4.1)
where L is a nonlinear differential operator which can be written

L = S s=1 L s . (4.2)
In Eq. (4.2) above, L s is a differential operator that represents a physical effect such as nonlinear phenomena, thermoviscous absorption or diffraction. For example solving for Eq. (2.19) (rewritten below as Eq. ( 4.3)) requires two successive steps, such that

D t R = -∂ x c 1 R + c 0 β 2 R 2 - c 0 2 ∂ 2 z R dx (4.3) = 2 s=1 L s R , in which ◮ L 1 R = -∂ x c 1 R + βc 0 2 R 2
, which accounts for nonlinear and refraction effects,

◮ and L 2 R = - c 0 2 ∂ 2 z R
dx, which accounts for propagation in the transverse direction (transverse propagation).

Note that nonlinear and refraction effects could as well be treated separately. Nevertheless, it is numerically more efficient to gather them together. The acoustic field can then be updated from time n to n + 1 by successive application of the differential operators.

The operator splitting method implicitly assumes that the different physical effects occurring in complex problems are uncoupled. Said differently, it is supposed that integration over sufficiently small propagation steps allows physical effects to be fully decoupled. This question has been investigated by Too & Lee [1995b] and [START_REF] Too | Numerical solution of the KZK equation for pulsed finite-amplitude sound beams in termoviscous fluids[END_REF] when solving a KZK equation.

The operator splitting method appeared to have no significant effect on the convergence of the combined solution. Furthermore, [START_REF] Tavakkoli | Modeling of pulsed finite-amplitude focused sound beams in time domain[END_REF] investigated the effects of inverting the order of the operators; it was shown to have no adverse effects on the numerical solution.

In this chapter the numerical solution of Eq. (4.3) is detailed. This Nonlinear Parabolic Equation (NPE) only contains nonlinear effects, refraction effects and diffraction. Section 4.2 presents the numerical method used for solving for the refraction and nonlinear effects term while Section 4.3 presents the method for solving for linear terms. Solving for more complete NPEs like Eq. (3.11) (propagation over non-flat terrains) or Eq. (3.36) (sound propagation within rigidly-framed porous media) can be done using the same numerical methods. Sections 4.3.2 and 4.3.3 briefly present modifications to perform to the original methods to solve for the complete NPE set. The treatment of initial and boundary conditions is described in Sections 4.4 and 4.5. Section 4.6 gives some aspects of software development and finally, chapter summary and conclusions are given in Section 4.7.

Calculation method for nonlinear terms with the Flux Corrected Transport algorithm

The PDE accounting for refraction and nonlinear effects is

D t R = -∂ x c 1 R + βc 0 2 R 2 , (4.4)
which can be written as a generic scalar conservation equation, such that

D t R + ∂ x (f (R)) = 0 , (4.5) in which the flux function is f (R) = c 1 R + βc 0 2 R 2 .
Solving for Eq. (4.5) with traditional finite-difference schemes would yield erroneous solutions:

discontinuities in the solution during shock formation must be handled with specialized algorithms. These so-called shock-capturing methods [see e.g. [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF] introduce some numerical viscosity in the scheme, allowing Gibb's oscillations to be reduced.

As an example, Fig.

∂ t u + 1 2 ∂ x u 2 = 0, -∞ < x < ∞, 0 < t , (4.6a) u (0, x) =      0 if x < 0 1 if x > 0 , (4.6b)
for two different numerical schemes, the MacCormack scheme and the leapfrog scheme.

Eqs. (4.6) is the inviscid Burger's equation. The similarity between Eq. (4.4) (the NPE)

and Eq. (4.6a) (the Burgers equation) should be noted: eliminating refraction effects (c 1 = 0)

and setting c 0 = β -1 in Eq. (4.4) gives Eq. (4.6a) subjected to specific initial conditions. As it can be seen in Fig. ( 4.1) both schemes do not behave correctly near the shock. It can be noted that while the leapfrog scheme has only numerical damping, the small amount of viscosity in the MacCormack scheme allows one to keep numerical oscillations near the discontinuity. It is thus possible to design numerical schemes that can handle shock formation and propagation without introducing spurious oscillations. The challenging task is to use the least possible numerical viscosity while suppressing Gibb's oscillations.

Among the shock-capturing schemes bestiary [see for example [START_REF] Harten | High Resolution Schemes for Hyperbolic Conservation Laws[END_REF][START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF][START_REF] Colella | Modern numerical methods for fluid flow[END_REF]], the Flux Corrected Transport (FCT) algorithm is second-order accurate and shows high-quality shock resolution. The FCT algorithm accomplishes this objective by combining integration schemes with low and high orders of spatial accuracy.

The low-order scheme provides a monotone solution, usually by the introduction of diffusive numerical fluxes, while the high-order scheme provides high accuracy in regions of smooth flow. The high-order solution is obtained by "anti-diffusing" the low-order, monotone solution, but only to such an extent that no new extrema are created and no existing extrema are accentuated. This is done by limiting, or correcting, the anti-diffusive fluxes of the high-order scheme, hence the name of the algorithm. Detailed mathematical description can be found in [START_REF] Boris | Flux-corrected transport III. Minimal error FCT algorithms[END_REF]; [START_REF] Harten | High Resolution Schemes for Hyperbolic Conservation Laws[END_REF]; [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF] and applications of the FCT algorithm to wave propagation in Védy [2002b,a].

Below are presented the main calculation steps that form the FCT algorithm. The following notations will be used:

∆ + R n i,j = R n i+1,j -R n i,j , (4.7a) ∆ -R n i,j = R n i,j -R n i-1,j , (4.7b)
where R n i,j denote the field value at range i∆x in the moving window, altitude j∆z and time n∆t. We define the numerical flux of a first order Total Variation Diminishing (TVD) scheme In Eq. (4.8), a i+ 1 2 ,j is the local wave speed defined by

h n i+ 1 2 ,j , such that h n i+ 1 2 ,j =      f R n i+1,j if a i+ 1 2 ,j < 0 f R n i,j if a i+ 1 2 ,j > 0 . ( 4
a i+ 1 2 ,j =          f ′ R n i,j if R n i+1,j = R n i,j f R n i+1,j -f R n i,j R n i+1,j -R n i,j otherwise , (4.9) in which f ′ R n i,j is the flux function derivative, i.e. f ′ (R) = ∂ R f (R) , (4.10a) = c 1 + c 0 βR . (4.10b)
A temporary field value R * i,j can be calculated with

R * i,j = R n i,j - ∆t ∆x ∆ -h n i+ 1 2 ,j .
(4.11)

Eq. ( 4.11) finalizes the first stage of the FCT method: it consists of a low-order transport step. Numerical errors are introduced in the temporary solution R * and have to be corrected by an anti-diffusive step.

The final acoustic field is obtained with

R n+1 i,j = R * i,j -b i+ 1 2 ,j -b i-1 2 ,j , (4.12)
where

b i+ 1 2 ,j =      0 if ∆ + R * i,j ∆ -R * i,j < 0 or ∆ + R * i+1,j ∆ -R * i+1,j < 0 sign ∆ + R * i,j min 1 2 ∆ -R * i+1,j , d i+ 1 2 ,j ∆ + R * i,j , 1 2 ∆ + R * i+1,j otherwise , (4.13) with d i+ 1 2 ,j = 1 2 Q i+ 1 2 ,j -Q LW i+ 1 2 ,j , (4.14)
in which Q i+ 1 2 ,j and Q LW i+ 1 2 ,j are numerical viscosities of upwind and Lax-Wendroff schemes, respectively. They are defined by:

Q i+ 1 2 ,j = ∆t ∆x a i+ 1 2 ,j , (4.15)
and

Q LW i+ 1 2 ,j = ∆t ∆x 2 a 2 i+ 1 2 ,j .
(4.16)

The two-stage method used by the FCT algorithm allows the numerical errors introduced in the first step (the transport stage) to be corrected in the second step (anti-diffusion step). This makes the FCT algorithm well-suited and efficient for shock wave propagation applications.

4.3 Calculation method for linear terms with a semi-implicit scheme

Solving for diffraction

The diffraction operator is

D t R = - c 0 2 x x f ∂ 2 z R dx , (4.17)
where x f is a point in the moving window ahead of the wave where the medium is quiescent, and x is the position where the solution is sought. Note the inverted integral limits in Eq. (4.17): the integration is done from the right side to the left side of the calculation grid, taking the right-most point as initial point.

The Crank-Nicolson method is used to solve Eq. (4.17): the scheme is second order in time, stable for all ∆t and is semi-implicit. It involves taking the spatial derivatives halfway between the beginning and the end of the current time step, such that

∂ k i R ≡ 1 2 ∂ k i R n+1 + R n . (4.18)
The Crank-Nicolson method is described in Appendix E.1.

Inverting integration limits in Eq. (4.17) and substituting spatial and temporal derivatives by their discrete forms gives

Q n+1 i,j -2R n i,j = c 0 ∆t 4∆z 2 x f x Q n+1 i,j+1 -2Q n+1 i,j + Q n+1 i,j-1 dx , (4.19)
in which

Q n+1 i,j = R n+1 i,j + R n i,j . (4.20)
A trapezoidal law is then used to evaluate the integral, such that

x f x Q n+1 i,j dx = ∆x 1 2 Q n+1 i,j + Q n+1 Nx,j + i+1 m=Nx-1 Q n+1 m,j , (4.21)
in which N x is the number of grid points in the window, and x f and x are defined by

x f = N x ∆x and x = i∆x. Since x f is a point ahead of the perturbation then Q n+1 Nx,j = 0 in Eq. (4.21). Substituting Eq. (4.21) in Eq. (4.19) yields

(1 + 2α) Q n+1 i,j -α Q n+1 i,j+1 + Q n+1 i,j-1 = 2R n i,j + 2α i+1 m=Nx-1 Q n+1 m,j+1 -2Q n+1 m,j + Q n+1 m,j-1 , (4.22)
where α = c 0 ∆t∆x 8∆z 2 . Eq. (4.22) can be written in a matricial form, such that

[D] q = b , (4.23)
where

[D] is a tridiagonal matrix, q = Q n+1 i,1 , . . . , Q n+1 i,j , . . . , Q n+1 i,Nz
T is the solution vector, and b is the right-hand-side term of Eq. (4.22). The matrix [D] and vector b are given by

[D] =        1 + 2α -α -α 1 + 2α -α . . . . . . . . . -α 1 + 2α -α -α 1 + 2α        , (4.24) and b =            2R n i,1 + 2α i+1 m=Nx-1 Q n+1 m,2 -2Q n+1 m,1 + Q n+1 m,0 . . . 2R n i,j + 2α i+1 m=Nx-1 Q n+1 m,j+1 -2Q n+1 m,j + Q n+1 m,j-1 . . . 2R n i,Nz + 2α i+1 m=Nx-1 Q n+1 m,Nz+1 -2Q n+1 m,Nz + Q n+1 m,Nz-1            . (4.25)
Numerical integration is made column by column, from the right side of the window to the left side, so that the terms appearing in the sums in Eq. (4.25) are known. Note that Eqs. (4.24, 4.25) do not contain any information on boundary conditions. It is for now assumed that boundary conditions on the top and on the bottom of the domain (points with indexes j = 1 and j = N z ) are somehow contained in the values R n i,Nz+1 and R n i,0 (boundary conditions are treated in Section 4.5). The tridiagonal linear equations system Eq. (4.23) is solved with a Thomas algorithm which is described in Appendix E.2. This method for solving for linear terms has been used in the works by [START_REF] Lee | Time-domain modeling of pulsed finite-amplitude sound beams[END_REF]; [START_REF] Coulouvrat | Propagation acoustique non linéaire en milieu inhomogène: de la modélisation à la simulation numérique (Nonlinear acoustic propagation in inhomogeneous media: from modelisation to numerical simulation[END_REF]; [START_REF] Marchiano | Nonlinear focusing of acoustic shock waves at a caustic cusp[END_REF] or more recently by [START_REF] Baskar | Nonlinear reflection of grazing acoustical shock waves : unsteady transition from von Neumann to Mach to Snell-References Descartes reflections[END_REF]; [START_REF] Marchiano | Experimental evidence of deviation from mirror reflection for acoustical shock waves[END_REF].

Non-flat terrains

D t R = - c 0 2 1 + h ′2 ∂ 2 z R + h ′′ ∂ z R dx + c 0 h ′ ∂ z R . (4.26)
The modified diffraction operator can be solved using the method presented in this section.

The main difference is that the left-hand-side matrix [D] (see Eq. (4.24)) becomes rangedependent, due to the presence of the ground topography derivatives h ′ (x) and h ′′ (x) under the integral term. This does not prevent previously presented methods to be used in the same way, but this has an impact on numerical efficiency.

Indeed, for flat terrains the matrix [D] could be calculated outside the time iteration loop, and more importantly the Lower/Upper (LU) decomposition performed in the Thomas Algorithm (see Appendix E.2) could be performed once and for all, the work remaining in the iteration loop being reduced to a forward and a backward substitution. When non-flat topographies are used this LU decomposition has to be repeated at each x-position and it hence increases the computational complexity (computational time is roughly multiplied by a factor 2).

Note that the new differential equation

D t R = c 0 h ′ ∂ z R , (4.27)
the last term in Eq. (4.26), should be treated together with the diffraction term and should not be split into a new differential operator. Indeed, the additional time needed to modify the left-hand-side matrix [D] is far much shorter than the time required to solve a second system of equations.

Solving for geometrical spreading and thermoviscous effects

The procedure to solve for geometrical spreading and thermoviscous effects is very similar to the one presented in Section 4.3.1: the Crank-Nicolson method and the Thomas algorithm (see Appendix E) are used. The reader may refer to [START_REF] Leissing | Nonlinear outdoor sound propagation -A numerical implementation and study using the nonlinear progressive wave equation[END_REF] to get a complete description of the implementation of these operators.

Solving the Nonlinear Parabolic Equation for sound propagation over porous ground surfaces

Solving for propagation within the porous layer

The NPE used to simulate wave propagation within porous ground surfaces is (see Eq. (3.36) in Section 3.2.1)

D t R = - c 0 √ Φ ∂ x 1 - √ Φ R + β 2 R 2 - c 0 2 √ Φ ∂ 2 z R dx - σ 0 Ω 0 2Φρ 0 1 + ξc 0 √ Φ |R| R . (4.28)
Hydrodynamic nonlinearities and moving-frame speed correction operators can be solved with the method presented in Section 4.2. The flux function and its derivative are modified according to

g (R) = c 0 √ Φ 1 - √ Φ R + c 0 β 2 √ Φ R 2 , (4.29)
and

g ′ (R) = c 0 √ Φ 1 - √ Φ + c 0 β √ Φ R . (4.30)
The diffraction term is identical to the one for the air layer; the method presented in Section 4.3.1 can directly be used by replacing c 0 with c 0 / √ Φ.

The absorption term in Eq. (4.28) can be solved in two stages: classical methods (Crank-Nicolson method and Thomas algorithm) are used to solve for the linear part. The Forchheimer nonlinearities operator can then be written

D t R = - σ 0 Ω 0 ξc 0 2Φ 3/2 ρ 0 sgn (R)R 2 . (4.31)
Eq. (4.31) is an autonomous nonlinear PDE which can be solved with any explicit finitedifferencing scheme. In the present work a first-order explicit differencing scheme is used. Step 2: backward integration over grid columns,

Step 3: calculation of left-hand-side matrix and right-hand-side vector,

Step 4: solution of the tridiagonal equation system with the Thomas algorithm.

These four steps must be repeated for each calculation grid column and at each time step.

Note that if backward integration is performed before updating the sums in Eqs. (3.53) (

Step

2 before
Step 1), the solver is highly unstable.

The same conclusions about numerical efficiency when the ground surface is non-flat can be made. Re-evaluation of coupling parameters and left-hand-side matrix at each range position is required, hence considerably increasing computational times.

For the case where Forchheimer's nonlinearities are of importance, an additional step must be performed before going through the procedure described above. Indeed, coupling parameters must be re-evaluated with updated values of flow resistivity. As explained in Section 3.

Treatment of initial conditions

Many methods can be used to initialize the NPE calculation. However, two of them covers all initialization cases and are detailed in this section.

The Kinney-Graham (KG) model has the advantage of being fast but has some important restrictions. Indeed, the analytical model gives time waveforms for blast wave propagation in free field. If the wave amplitude is very important one has to initialize the NPE calculation at a moderate distance so that nonlinear effects are weak (weak nonlinearities limitation of the NPE, see Section 2.1.3). The KG model being unable to account for ground effects or meteorological effects, another method has to be used instead if the configuration exhibit sound speed gradients or a porous ground surface. If the wave amplitude is low enough, then the KG model can be used. Several time waveforms have to be generated at different times so that a spatial waveform can be interpolated from these signals. Once the spatial waveform is obtained it can be spherically extrapolated to obtain a starting grid for the NPE simulation.

The numerical solution of Euler equations can be used for initialization in the far-field. This methods has the advantage of being accurate: two-dimensional simulations including the effects of meteorological conditions, dissipation effects and ground effects can be performed and moreover, the method does not suffer from the weak nonlinearities limitation. The main drawback is calculation times that are often on the order of hours for two-dimensional simulations. Note that the numerical solution of Euler equations directly provides two-dimensional arrays.

To put it simply, the analytical Kinney-Graham (KG) model can be used to initialize the NPE calculations in the near-field while the solution of Euler equations can be used in the far-field. The KG model is presented in Appendix B while the numerical solution of Euler equations is described in Appendix C.

Treatment of boundary conditions

In this section the boundary conditions used on the calculation grid are detailed. Three types of boundary conditions are needed:

◮ lateral boundaries: right-most and left-most points of the grid, see Section 4.5.1, ◮ the bottom row of the domain, see Section 4.5.2, ◮ and the absorbing layer at the top of the propagation domain, see Section 4.5.3.

Lateral boundary conditions

The right-most column of the computational grid is a point ahead of the perturbation, where the medium is quiescent. This property has been used in Eq. (4.21) to eliminate the field value Q n+1 Nx,j . One hence has

Q n Nx,j = 0 , j = 1, . . . , N z , ∀n . (4.32)
On the left-most column of the computational grid we assume that the perturbation is weak,

i.e. one has

Q n 1,j = 0 , j = 1, . . . , N z , ∀n . (4.33)
Note that these boundary conditions on the lateral sides of the computational grid induce some restrictions on the choice of the moving-window speed and the resulting sound speed perturbation. Indeed the right-most point was chosen as a point ahead of the perturbation.

This implies where c win is the moving-window speed.

c 1 (x, z) ≤ 0 , ( 4
Similarly, the window has to be chosen sufficiently wide with respect to the sound speed perturbations. As a rule of thumb one must have

L > T max (|c 1 (x, z)|) (4.36)
in which L is the window width and T is the time the wave propagates. This way the signal stays into the computational grid.

One could think about the opposite solution: defining c win = min (c 1 (x, y)) gives a positive sound speed perturbation and thus a forward propagation in the window. This solution may appear simpler but is not necessarily a good choice. The signal would have to be placed somehow in the middle of the window so that diffraction effects can occur in the left side of the window and that wave propagation, due to sound speed perturbations, can occur in the right side of the window. This solution would require a larger window and would thus significantly increase computational times.

Boundary condition on the bottom of the domain

The bottom row of the ground layer is assumed to be acoustically rigid. A rigid boundary condition ∂p/∂n = 0 is incorporated by assuming a solution which is symmetric with respect to the boundary and modifying the recurrence accordingly for the boundary nodes. Note that this boundary condition only appears in the diffraction operator; the other Partial Differential Equations (PDE) do not use pressure value at location (i, j -1).

Assuming that the bottom row of the computational domain has an index j = 1, the reformulation of Eq. (4.22) including the rigid boundary condition gives an expression for the field on the bottom of the ground layer, such that

(1 + 2α) Q n+1 i,1 -2αQ n+1 i,2 = 2R n i,1 + 4α i+1 m=Nx-1 Q n+1 m,2 -Q n+1 m,1 .
(4.37)

The diffraction matrix [D] and right-hand-side vector b become

[D] =        1 + 2α -2α -α 1 + 2α -α . . . . . . . . . -α 1 + 2α -α -α 1 + 2α        , (4.38a) b =            2R n i,1 + 4α i+1 m=Nx-1 Q n+1 m,2 -Q n+1 m,1 . . . 2R n i,j + 2α i+1 m=Nx-1 Q n+1 m,j+1 -2Q n+1 m,j + Q n+1 m,j-1 . . . 2R n i,Nz + 2α i+1 m=Nx-1 Q n+1 m,Nz+1 -2Q n+1 m,Nz + Q n+1 m,Nz-1            . (4.38b)

Domain truncation

This section briefly reviews available solutions to truncate the computational domain and then details the development of a Perfectly Matched Layer (PML) for NPE models.

Absorbing Boundary Conditions (ABC) and absorbing layers

Whenever one solves a PDE numerically by a finite-difference method (or any other numerical method) one must truncate the computational domain. The key question is how to perform this truncation without introducing significant numerical oscillations into the solution. The oscillating nature of wave equation solutions makes the truncation of the computational domain a difficult task. Simply truncating the grid with hard-wall conditions (see Section 4.5.2) would produce unacceptable artifacts from boundary reflections. It means that any real coordinate remapping from an infinite to a finite domain will result in solutions that oscillate infinitely fast as the boundary is approached. Obviously, the larger is the absorbing layer the better it will be, at the expense of computational time. Therefore wave equations require something different: an absorbing boundary that will somehow absorb waves that strike it, without reflecting them, and without requiring infeasible resolution.

The first attempts at such absorbing boundaries for wave equations involved Absorbing Boundary Conditions (ABC) [see for example [START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave propagation[END_REF][START_REF] Peng | An optimal absorbing boundary condition for finite difference modeling of acoustic and elastic wave propagation[END_REF]]. An ABC tries to somehow extrapolate from the interior grid points to the edge grid points, to fool the solution into "thinking" that it extends forever with no boundary. It turns out that this is possible to do perfectly in one dimension. However, the main interest for numerical simulation lies in two and three dimensions, and in these cases the infinite number of possible propagation directions makes the ABC problem much harder. Although twodimensional and three-dimensional formulations of ABC have been developed over the past decades [see for example [START_REF] Bogey | Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation cases[END_REF][START_REF] Sparrow | Absorbing boundary conditions for a spherical monopole in a set of two-dimensional acoustics equations[END_REF][START_REF] Cao | An exact absorbing boundary condition and its application to three-dimensional scattering from thin dispersive structures[END_REF],

existing ABC restrict themselves to absorb waves exactly only at a few angles, especially at normal incidence. Moreover many standard ABC are formulated only for homogeneous materials at the boundaries. For the application of this work, propagation of blast waves through large distances, waves are highly nonlinear, propagate in an inhomogeneous medium and hit the top of the computational domain with a very small grazing angle. Despite the development of an ABC for Parabolic Equation (PE) models by [START_REF] Yevick | Nonlocal boundary conditions for finite-difference parabolic equation solvers[END_REF][START_REF] Mcdonald | High-angle formulation for the nonlinear progressive wave equation model[END_REF], the development of a ABC for NPE models is a very challenging task.

In 1994 Bérenger proposed an absorbing layer instead of an ABC to truncate computational domains. An absorbing layer is composed of artificial absorbing material of a given thickness that is placed adjacent to the edges of the grid. When a wave enters the absorbing material it is attenuated by the absorption. The problem with this approach is that the transition from the physical domain to the absorbing material generally introduces spurious numerical reflection. However, Bérenger showed that an absorbing layer with special properties can be constructed, so that waves do not reflect at the interface. The Perfectly Matched Layer (PML) was originally derived for electromagnetic wave equations in the frequency domain (Maxwell's equations) but has later been adapted and used for acoustic wave equations [see [START_REF] Hu | Development of PML absorbing boundary conditions for computational aeroacoustics: A progress review[END_REF] in the time-domain [see [START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF].

Basic principle of the Perfectly Matched Layer (PML)

Consider a wave propagating in a homogeneous and time-invariant media. The solution can be decomposed into a superposition of plane waves such that

p (x, t) = k,ω P k,ω e i(k•x-ωt) , (4.39)
in which x is the position vector, k is the wave vector, ω is the angular frequency and P k,ω are some constant amplitudes. Let's suppose that one wants to truncate the domain in the x-direction; Eq. (4.39) can be written p (x, t) = P (y, z) e i(kx-ωt) . (4.40)

Bérenger showed that the coordinate mapping (4.41) in which x 0 and δ are the the position of the layer and its thickness, respectively, allows us to transforms Eq. ( 4.40) into

x -→ x + i ω x 0 +δ x 0 σ (x) dx ,
p (x, t) = P (y, z) e i(kx-ωt) e -k ω R x σ ′ (x ′ ) dx ′ . (4.42)
One can notice that provided σ (x) > 0 waves are exponentially decaying. In theory the layer thickness can be chosen infinitely small. In practice the layer thickness is often chosen to be equal to a wavelength and σ (x) quadratically increases from zero to a given value.

A Perfectly Matched Layer (PML) for NPE models Application of the PML to Parabolic Equation (PE) models can be found in [START_REF] Collino | Perfectly Matched Absorbing Layers for the Paraxial Equations[END_REF].

We use here the results of Collino to derive a Perfectly Matched Layer (PML) for (timedomain) Nonlinear Parabolic Equation (NPE) models. A PML in the z-direction is created by introducing the following change of variable:

z -→ z + i ω z 0 +δ z 0 σ (z) dz . (4.43)
Assuming a e -iωt time dependence, the partial derivative with respect to variable z is changed according to

∂ z -→ iω iω + σ (z) ∂ z . (4.44)
We start from a simplified NPE formulation: nonlinear effects are neglected and it is assumed that there is no sound speed perturbations, i.e. c 1 (x, z)= 0. The NPE writes

D t R = - c 0 2 ∂ 2 z R dx . (4.45)
Rearranging Eq. (4.45) and substituting the change of variables in Eq. (4.44) yields

iω∂ x R + c 0 2 iω iω + σ (z) ∂ z iω iω + σ (z) ∂ z R = 0 . (4.46)
A few algebraic manipulations of Eq. ( 4.46) and a transformation back to the time-domain give a NPE/PML model formulation, defined by

D t R + c 0 2 ∂ 2 z R dx = -3σ (z) R -3σ (z) 2 R dt + c 0 2 σ (z) ∂ 2 z R -∂ z (σ (z)) ∂ z R dx dt -σ (z) 3 R dt . (4.47)
Note that if σ (z) = 0 Eq. ( 4.47) reduces to the usual NPE Eq. (4.45). Similarly, NPE/PML models can be derived for propagation in inhomogeneous media or for the terrain-following coordinates formulation of the NPE model (see Eq. (3.11)).

In the derivation of Eq. ( 4.47) nonlinear effects were neglected. Although PMLs have been developed for nonlinear wave propagation [see for example [START_REF] Appelö | Application of a perfectly matched layer to the nonlinear wave equation[END_REF], it has been found that high-amplitude effects have no adverse consequences on the performances of the PML. The following NPE/PML formulation was thus used for high-amplitude wave propagation:

D t R + βc 0 2 ∂ x R 2 + c 0 2 ∂ 2 z R dx = -3σ (z) R -3σ (z) 2 R dt + c 0 2 σ (z) ∂ 2 z R -∂ z (σ (z)) ∂ z R dx dt -σ (z) 3 R dt . (4.48)
The parameter σ (z) is set to

σ (z) = c 0 2δ z -z 0 δ 2 , z ≥ z 0 , (4.49)
in which δ is chosen to be equal to the main wavelength of the signal under interest.

The designed PML allows one to truncate the computational grid without introducing spurious numerical oscillations in the domain. Near the source, for relatively high incidence angles, the layer thickness can be chosen very thin: 5 points are often sufficient to completely absorb waves while at very low incidence angles, far from the source, the thickness has to be increased to be as large as a wavelength. From a numerical point of view the NPE/PML model shown in Eq. (4.48) introduces some difficulties: the presence of time integrals implies keeping in memory previous solution fields. Fortunately the layer thickness is seldom larger than a few dozen points, keeping the addition numerical effort to integrate the terms acceptable.

An example of a simulation using the NPE/PML model can be observed in The contours represent the wave amplitude in decibels, relative to the maximum amplitude of the wave. As it can be seen, no reflections from the PML can be observed, meaning that the amplitude of the spurious oscillations is at least 70 dB lower than the maximum amplitude of the wave. 4.6 Some aspects of software development

Choice of programming languages

Two programming languages were used to develop the simulation software. Python * was used to drive computational routines written in Fortran 95. This combination of a high-level programming language such as Python with Fortran allows one to take advantage of each language strengths and weaknesses. Python and Fortran are bound together thanks to the f2py program † which automatically generates python interfaces to Fortran routines.

Program structure

The

Chapter summary and conclusions

In this chapter the numerical methods used to solve the various NPE models have been presented. This numerical solution assumes that the different physical effects are decoupled over a time step, so that the different operators are split. Two finite-difference schemes are then used to solve for these different operators. The Flux Corrected Transport (FCT) algorithm solves for nonlinear effects while the Crank-Nicolson method is used to solve for linear terms.

The FCT method uses a two-stage method that allows us to stably propagate shocks without introducing numerical oscillations. The Crank-Nicolson method is semi-implicit and yields tridiagonal systems of equations that are solved with a Thomas algorithm.

Propagation of high-amplitude quasi-plane waves in a homogeneous medium

In this section the propagation of high-amplitude waves is simulated with the NPE and numerical solutions are compared with analytical results. The propagation medium is homogeneous and the ground surface is flat and acoustically rigid. Since the nonlinear differential operator appearing in the NPE is one-dimensional we consider the propagation of quasi-plane waves, i.e. the NPE model is reduced to

D t R = -∂ x c 0 β 2 R 2 .
(5.1)

The Fubini solution

The analytical solution used for comparison to the NPE is the solution proposed by [START_REF] Fubini | Anomalie nella propagazione di ande acustiche de grande ampiezza (Anomalies in the propagation of high-amplitude acoustic waves)[END_REF] that gives pressure waveforms for propagation of high-amplitude harmonic waves. The pressure distribution is expressed by an infinite sum of weighted Bessel functions, such that *

p (σ, τ ) = p 0 ∞ n=1 2 nσ J n (nσ) sin (nωτ ) , (5.2)
where x is the shock formation distance, σ= x/x, τ = t -x/c 0 , p 0 is the initial sinusoid amplitude and J n is the Bessel function of order n. The formulation presented in Eq. ( 5.2)

gives valid solutions for one-dimensional problems only in the pre-shock region (σ< 1). Note that losses are not included in the Fubini solution described by Eq. (5.2).

Configuration

Results of nonlinear calculations given by the NPE are compared to Fubini's analytical solutions. The sound speed is constant in the propagation domain and is set to c 0 = 343 m.s -1

and the nonlinearity coefficient is β = 1.2 (air under normal atmospheric conditions). The initial condition for starting the NPE simulation is a quasi-plane sine pulse of frequency 0.1

Hz and peak amplitude 10 kPa (174 dB re 20 µPa). The waveform is recorded at each time iteration. The spatial resolution is 100 points/λ, resulting in a time step ∆t = 100 ms. The NPE simulation ran for 500 time iterations.

Results

The wave steepening process can be observed in Results given by the implemented NPE model are close to the analytical solution. The maximum relative error is 0.4 % for the fundamental frequency amplitude. For the first and second harmonics the maximum relative errors are 2 % and 3.5 %, respectively, provided their values are much greater than numerical uncertainty. Indeed, one may note that the relative difference between the models may be rather high (for example: third harmonic amplitude, at the beginning of the simulation). This is due to the very low harmonics amplitudes at the beginning of the simulation. As a result relative errors are large but as the wave steepens it quickly decreases. Nonlinear effects are accurately simulated by the implemented NPE model.

Sound propagation in an inhomogeneous medium over a flat and acoustically rigid ground surface

This section presents validation results for low-amplitude wave propagation in an inhomogeneous medium over a flat and acoustically rigid ground surface. The NPE simulations are compared to the numerical solutions of the (linear, frequency-domain) Parabolic Equation (PE).

ATMOS sound propagation code: an implementation of the linear, frequency-domain Parabolic Equation

The PE implementation used is called Advanced Theoretical Model for Outdoor Sound propagation (ATMOS). Designed and developed by Centre Scientifique et Technique du Bâtiment (CSTB) researchers [see [START_REF] Aballéa | Propagation acoustique en milieu extérieur : Application de l'équation parabolique rapide au couplage d'effets météorologiques et de topographies complexes (Outdoor Sound Propagation: Application[END_REF], ATMOS is a software which determines the acoustic impact of all types of infrastructures, both near and far away, on their environment, taking the geometry of the site into account as well as the meteorological factors inherent to each local climate. Furthermore, ATMOS is also being used as a reference numerical model within The Transmission Loss (TL) is defined as the ratio in decibels between the acoustic intensity I (x, z) at a field point and the intensity I 0 at 1 m distance from the source, such that

T L = -10 log I (x, z) I 0 (5.3a) = -20 log | p (x, z) | p 0 [dB re 1m] .
(5.3b)

Coupling the Nonlinear Parabolic Equation model with the frequencydomain Parabolic Equation

So that ATMOS can be used as a comparison tool the initial conditions used in the NPE calculations (spatial pressure distributions) have to be transformed into frequency-domain "starters" that are used as initial conditions in the ATMOS software. The pressure signal is recorded along a vertical line and transformed into the frequency-domain with the help of a Fourier transform (see Section 2.5).

This coupling method is validated with a simple configuration. The (linear) propagation of a harmonic source is studied both with the NPE and ATMOS. The sound speed is constant with altitude and equal to 340 m.s -1 and the source, a quasi-plane wave sine pulse with a frequency of 0.1 Hz, is placed on the ground. The coupling between the models is performed at two distances: 46 km and 112 km. One could think of using the NPE to handle nonlinear effects in the near field and at moderate distances from the source and then use the coupling method and finish the calculations with ATMOS (see Section 2.5 and Fig.

(2.2)). The most convenient (and most used) criteria to initiate the coupling is the maximum amplitude of the signal. The limit under which one can consider that nonlinear effects are weak is usually chosen as 100 Pa (134 dB re 20 µPa). The advantage is that ATMOS is faster for single frequency calculations. Indeed the Parabolic Equation (PE) is a very fast method but for broad-band signals the calculation has to be repeated for each frequency. This may result in extremely long calculations for high-resolution spectra. Many cases were tested and it appeared that unless a specific and reduced frequency range is of interest it is much faster to use the time-domain NPE, even when the signal exhibits low amplitudes.

Configuration

The validation procedure for propagation in an inhomogeneous medium now follows. The

Results

Transmission Loss (TL) is calculated for both the NPE and ATMOS simulations at the central frequency of the wave (10 Hz). The TL map from the NPE calculation is shown in Fig. (5.5(a)).

One can observe on this figure downward refraction effects. The difference in TL between both As it can be seen, on most of the domain the difference between ATMOS and the NPE model is acceptable. A zone near the source presents higher differences: around 3∼4 dB. This zone of lower accuracy is due to the coupling method: as shown in Section 5.2 the parabolic approximation prevents obtaining accurate results near the coupling location. Some small and localized areas also present variations higher than 2 dB but the relative difference never exceeds 4 dB.

Sound propagation over a non-flat acoustically rigid ground surface in a homogeneous medium

This section presents validation results for low-amplitude wave propagation in a homogeneous medium over an acoustically rigid non-flat ground surface. The reference model used for 1998]. The BEM and the MICADO software are briefly presented in Appendix F.

Configuration

The Generalized Terrain -Nonlinear Parabolic Equation (GT-NPE) model (see Eq. (3.11) in Section 3.1.2) is used to study the propagation of a finite-length signal over a hill. The configuration of this example is taken from the article by [START_REF] Karle | A linearized eulerian finite difference time domain sound propagation model with terrain following coordinates[END_REF].

The ground elevation is given by

h (x) = h top 1 + (x -x 0 ) 2 l 2 -1 , (5.4)
where h top is the maximum ground elevation, equal to 5 m, l is the hill width, equal to 15 m and x 0 is the maximum elevation location, equal to 70 m. The maximum slope is max h ′ (x) = 0.22. Tab. (5.1) summarizes the characteristics of the hill. The source is positioned at x = 0 m, z = 50 m and emits a 8-period 50 Hz sine wave pulse train. The sound speed c 0 is constant through the domain and set to c 0 = 343.4 m.s -1 . Spatial steps ∆x and ∆z are both equal to 30 cm, yielding a spatial resolution of approximately 22 points/λ and a time step ∆t = 0.88 ms. The ground layer is 20 point thick (60 cm) and acoustically rigid.

The propagation domain is 100 m high and 500 m wide.

Results

Snapshots of the wave propagation can be seen in close to the BEM calculations in the far-field. The parabolic approximation inherent to the GT-NPE model prevents obtaining correct SPLs near the source. Sound pressure levels (SPL) at three different altitudes (1.5 m, 15 m and 30 m) can be seen in Fig.

(5.9). In the far-field region good agreement is found between BEM and GT-NPE calculations. The difference is on the order of 1 dB, except for receiver positions where destructive interference occurs. This is mainly due to the numerical oscillations on the trailing part of the recorded signal inherent to time-domain algorithms. These low-amplitude oscillations prevent the negative interferences to be as important as they are with frequency-domain methods.

5.4 Sound propagation over a porous ground surface in a homogeneous medium

In this section numerical examples of sound propagation over porous ground layers are presented to illustrate the coupling method derived in Chapter 3 and to evaluate its performances. Linear propagation examples are presented in Section 5.4.1 while high-amplitude

wave propagation examples appear in Section 5.4.2. Linear propagation results are compared to analytical solutions and nonlinear cases are compared to numerical solutions of an implementation of the Euler equations (see Appendix C).

Linear propagation

This section presents validation results for linear propagation over a flat, finite-impedance ground surface in a homogeneous medium.

Description of the NPE model

Since there is no high-amplitude effects and no meteorological effects in this example the NPE model is defined by

D t R a = - c 0 2 ∂ 2 z R a dx in the atmospheric layer,
(5.5a) in which R a and R g denote field values in the atmospheric layer and the ground layer, respectively. Note that there is no geometrical spreading term included in the NPE model defined by Eqs. (5.5).

D t R g = - c 0 √ Φ ∂ x 1 - √ Φ R g - c 0 2 √ Φ ∂ 2 z R g dx - σ 0 Ω 0 2Φρ 0 R g in the ground layer, (5.5b) √ Φ∂ z R a - σ 0 Ω 0 ρ 0 c 0 ∂ z R a dx = Ω 0 ∂ z R g on the interface, (5.5c)

Reference solutions

The solutions of the two-dimensional Helmholtz equation are used as references. In the case of two-dimensional wave propagation over a flat and finite-impedance ground surface the pressure is given by

p r = iπH (1) 0 (kR 1 ) + QiπH (1) 0 (kR 2 ) , (5.6)
where p r is the complex pressure at the receiver, k is the wave-number, R 1 and R 2 are the source-receiver and image source-receiver distances, respectively, and H

(1) 0

is the Hankel function of the first kind and of order zero. The cylindrical reflection coefficient is noted Q;

it can be calculated with the help of Laplace transforms [see [START_REF] Di | An exact Laplace transform for a point sorce above a ground surface[END_REF][START_REF] Salomons | Eulerian time-domain model for sound propagation over a finite impedance ground surface. Comparison with frequency-domain models[END_REF]. The normalized impedance used to calculate the reflection coefficient is [see Salomons et al. , 2002, and Section 1.1.3]

Z = Φ Ω 2 0 + i σ 0 ρ 0 Ω 0 ω . (5.7)

Configuration

The sound speed is constant within the domain (c 0 = 340 m.s -1 ) and absorption from air is not included in the model. The source is positioned at (x s , z s ) = (0.0, 1.4) m. The signal used is a sine pulse with wavelength λ = 0.27 m (f = 1259.25 Hz) and its peak amplitude is low enough for the propagation to be considered linear. A receiver is placed 10 m away from the source at (x r , z r ) = (10, 1.4) m. The source and receiver positions ensure that we are within the parabolic equation angular validity domain (the angle from source to image-receiver is θ ≈ 15 • ). Spatial steps are equal to 7.5 10 -3 m in both directions, thus giving a spatial resolution of about 36 points/λ, ensuring sufficient resolution at higher frequencies and near the air/ground boundary. The time step is ∆t = ∆x/c 0 , so that at each time step the window advances one spatial step. Since the Crank-Nicolson method is used the numerical scheme is stable (see Appendix E.1). Three different ground layers of thickness 1 meter are considered.

The first ground layer is a perfectly rigid surface (Φ ≫ 1). The second and third layers have identical tortuosity (Φ = 3) and porosity (Ω 0 = 0.3), but different flow resistivity values (σ 0 = 500 kPa.s.m -2 and σ 0 = 100 kPa.s.m -2 ). The NPE window including the ground layer is 3 meter wide and 4.125 meter high (400 by 550 points).

Results

Two modifications are done on signals at the receiver. First, in order to obtain a free field reference, time histories are cropped after the direct wave. Next the trailing part of time signals is cropped after the reflected wave to suppress the low-amplitude numerical oscillations.

Let's denote the complex pressures at the receiver and in the free field by p r and p f ree , respectively. The relative Sound Pressure Level (SPL) ∆L is then calculated with Figure 5.10: Sound Pressure Level (SPL) relative to free field at the receiver for the three different ground layers, for both NPE and analytical solutions. Markers: NPE calculations; lines: analytical solutions. Blue: rigid ground layer; red: ground layer with σ 0 = 500 kPa.s.m -2 ; green: ground layer with σ 0 = 100 kPa.s.m -2 .

∆L =

Nonlinear propagation

This section presents validation results for high-amplitude wave propagation over a flat but finite-impedant ground surface in a homogeneous medium.

Description of the NPE model

The NPE model used is similar to the one in Eqs. (5.5) with nonlinearities included (hydrodynamic and Forchheimer's nonlinearities).

Reference solutions

To obtain reference results to compare to the NPE model simulations in the case of highamplitude waves, solutions of the Euler equations are used. The computational domain is composed of an air layer and a ground layer. In a two-dimensional Cartesian coordinate system the constitutive equations for the air layer are

∂ t ρ T + ∂ x (ρ T u) + ∂ z (ρ T w) = 0 , (5.9a)
∂ t (ρ T u) + ∂ x ρ T u 2 + ∂ z (ρ T uw) = -∂ x p T , (5.9b)
∂ t (ρ T w) + ∂ x (ρ T uw) + ∂ z ρ T w 2 = -∂ z p T , (5.9c)
∂ t (ρ T e 0 ) + ∂ x (ρ T ue 0 ) + ∂ z (ρ T we 0 ) = -∂ x (p T u) -∂ z (p T w) , (5.9d)
where e 0 is the energy per unit mass. Within the ground layer momentum conservation equations write

Φ∂ t (ρ T u) + ∂ x p T + Φρ T u 2 + ∂ z (Φρ T uw) + σ 0 Ω 0 (1 + ξ |u|) u = 0 , (5.10a) Φ∂ t (ρ T w) + ∂ z p T + Φρ T w 2 + ∂ x (Φρ T uw) + σ 0 Ω 0 (1 + ξ |w|) w = 0 .
(5.10b)

The energy equation Eq. (5.11) and the ideal gas law Eq. (5.12) close the equation system.

They are defined by

ρ T e 0 = ρ T c v T + ρ T |V| 2 2 , (5.11) p T = ρ T RT , (5.12)
where T is the gas temperature, c v is the specific heat capacity at constant volume and R is the gas constant. To solve this equation system a WENO algorithm [see [START_REF] Shu | Essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]] for space discretization and a third-order TVD scheme [see Press et al. , 1996a] for time marching are used. These numerical algorithms are briefly presented in Appendix C.

Configuration

In this example standard atmospheric conditions are used (T = 293 K, ρ 0 = 1.2 kg.m -3 , p 0 = 1.03 10 5 Pa). The source is positioned at (x s , z s ) = (0, 3) m and the receiver at

(x r , z r ) = (12, 3) m.
In order to start the reference calculation, the pressure, velocity, density and energy fields need to be specified. A Gaussian pulse is propagated using a one-dimensional version of the code presented in Appendix C. By adjusting the pulse amplitude and width one can obtain a one-dimensional signal at a given distance. In this example an amplitude and signal length of approximately 4 kPa and 1.5 m, respectively, were aimed for at a distance of 3 m from the source. Spatial steps ∆x and ∆z are both equal to 0.015 m, yielding a resolution of approximately 100 points per wavelength. This signal is then spherically extrapolated to obtain a two-dimensional array. To evaluate the accuracy of the NPE model, some characteristics of the reflected wave are studied: the maximum positive and negative peak pressures and their arrival times (noted respectively p + and p -, t a+ and t a-), and the positive phase duration (noted t d ). These characteristics are summarized in Tab. (5.2). Since for the softest ground layer the negative peak on the reflected wave is very weak, values of p -and t a-for this layer are irrelevant. As one can see, arrival times differ by at most 0.3 ms. The difference is larger for the softest layer; a possible reason is that the NPE model does not smear out pulses as the reference model does, yielding erroneous positive peak position. One can thus expect that as the flow resistivity decreases the error on arrival time increases. However, in outdoor sound propagation applications the flow resistivity may seldom be lower than the one used here (σ 0 = 10 kPa.s.m -2 ), so the error on arrival time will remain weak for most cases. These remarks are also applicable to the positive phase duration t d . Positive peak amplitudes differ by 6.2 % and 5.1 % for layers with σ 0 = 100 kPa.s.m -2 and σ 0 = 10 kPa.s.m -2 , respectively. This difference does not seem to be dependent on flow resistivity and as a comparison, the relative error for the perfectly rigid layer is 1 %. Relative error for negative peaks are comparable:

2.4 % and 3.4 % for the rigid layer and the layer with σ 0 = 100 kPa.s.m -2 , respectively. [START_REF] Plewa | Adaptive Mesh Refinement -Theory and Applications[END_REF], or moving window principle, see [START_REF] Sparrow | A numerical method for general finite amplitude wave propagation in two dimensions and its application to spark pulses[END_REF]), the NPE model, thanks to the use of a single-variable one-way wave equation and a fast solver (Thomas algorithm, see Appendix E), is a very efficient tool for outdoor sound propagation simulations.

Nonlinear propagation with Forchheimer's nonlinearities

To illustrate the effects of Forchheimer's nonlinearities a simulation is performed with a Forchheimer nonlinearity parameter ξ = 2.5 s.m -1 . According to the conclusions of the previous section, a low flow resistivity yields a larger error on the positive phase duration and on the time of arrival of the positive peak pressure. A low flow resistivity has been chosen

(σ 0 = 10 kPa.m.s -2 , with Φ = 3, Ω 0 = 0.3) so that the method to include Forchheimer's nonlinearities can be fully evaluated. Simulation parameters and initialization array are identical to those used in Section 5.4.2. The source is positioned at (x s , z s ) = (0, 3) m and the receiver at (x r , z r ) = (12, 3) m.

Chapter summary and conclusions

In this chapter different validation cases have been presented.

First, the NPE model is compared to analytical solutions to assess the accuracy of nonlinear effects calculations. The implemented NPE has been proved to accurately simulate wave steepening and harmonics generation with a relative error lower than 5 % until the third harmonic. Quasi-plane waves were considered in this example, but the paraxial approximation of the NPE ensures the validity of nonlinearity calculations for two-dimensional configurations.

The computational model has then been tested regarding meteorological conditions and compared to a frequency-domain implementation of the Parabolic Equation (PE). It appeared that Sound Pressure Levels (SPL) are accurate within a 2 dB relative error range except on small localized areas, where the difference with the PE method is on the order of 3 dB.

This validation case also proved the feasibility of coupling the NPE model to other propagation codes, in this particular case the frequency-domain PE. The coupling between the two different numerical methods has been shown to be successful.

The BEM is then used to perform a simulation on a hilly terrain. It has been shown that the terrain-following coordinates formulation of the NPE, the Generalized Terrain -Nonlinear Parabolic Equation (GT-NPE) model, is in good agreement with the BEM results in the far-field. This model will allow to study wave propagation in hilly urban areas, once the probabilistic approach of uncertainties is developed.

The performances of the NPE model for propagation over porous ground layers are then assessed. For linear propagation, the results obtained with this method have shown very good agreement with analytical solutions for a wide range of ground properties. For highamplitude waves, the NPE model produces time signals comparable to those obtained by the numerical solution of Euler equations. Relative error on peak pressures has been shown to be independent on material properties while differences on positive phase duration and time of arrival increases with decreasing ground flow resistivity. However, the presented model still gives good agreement even for very low flow resistivity values and provides a simple but efficient way of taking into account ground impedances. The approximate method to include the Forchheimer nonlinearities in the two-way coupling has then been proved to give satisfactory results and does not introduce any additional source of error in the two-way coupling.

To construct the NPE model, the assumption that the ground layer is equivalent to a continuous fluid has been made. This simplified approach allows one to derive a ground model that is of the same form as the NPE model for atmospheric layer. Since the two-way coupling equations involve only spatial derivatives and integrals, the complete NPE model is able to perform simulations in very short times (about 50 times faster than the Euler equations implementation used in this work). This enables the NPE model to be used as a stochastic model solved by the Monte-Carlo method. Wave field statistics in the air layer could be determined by performing a large number of simulations of sound propagation in an environment with varying parameters (e.g. propagation over a ground layer with random flow resistivity, propagation through turbulence). However, note that for realistic simulations, a spherical spreading term should be added to the NPE used in this work.

This chapter finalizes the development of a deterministic sound propagation model for high amplitude wave propagation in complex media. The model can account for most of the features of sound propagation outdoors (nonlinearities, refraction, hilly terrain and ground impedance effects) and has been shown to provide results that are in very good agreement with other numerical models. Furthermore, the fastness of the NPE model makes it a good candidate for a basis of a stochastic sound propagation model.

6

Application to stochastic high-amplitude wave propagation over urban environments using the Nonlinear Parabolic Equation modelThis chapter concerns the development of a computational model for long-range nonlinear sound propagation over urban environments. Section 6.1 briefly summarizes the previous chapters, and Section 6.2 details the methodology for constructing the computational model.

Next the probability model of the geometrical parameters of a urban environment are determined using Information Theory and the Maximum Entropy Principle (Section 6.3) and then the reference model is described (Section 6.4). The mean propagation model is then presented: it is based on the Nonlinear Parabolic Equation (NPE) model and its extension to propagation in porous media (Section 6.5). In Sections 6.6 and 6.7 the construction of the probabilistic model of the stochastic propagation model is detailed. In Section 6.8 the method used to identify its parameter is presented. Section 6.9 details two applications using the constructed stochastic propagation model and finally, chapter summary and conclusions are given in Section 6.10. the urban city is taken into account through independent random porous layer parameters, noted Γ, Λ and Θ. The probability distributions of these random variables depend on a parameter noted w. The model output is the pressure at the receiver which is noted P r (ω).

Chapter content

A second propagation model, in which the urban city is explicitly accounted for, is used to provide reference solutions. Given a parameter u which characterizes the urban city geometry and several probability models, a urban city realization is generated, and then used in the reference model to obtain the pressure P exp r (ω) at the receiver. To construct the stochastic sound propagation model, for a given parameter u, one has:

Step 1: to construct the probability models of the urban city geometrical parameters. Since from one city to another one, geometrical parameters (e.g. building density, mean elevation, etc.) can greatly vary, the construction of a probabilistic model of these parameters should be done with measured data (for example, data from a geographical information system). Moreover, town planning and buildings themselves complies with many constraints which introduce a statistical dependence between the geometrical parameters. The determination of these dependencies being out of the scope of this work, probability models are constructed assuming that no information is available concerning the relations between geometrical parameters. The probability models are determined with the help of Information Theory [see [START_REF] Shannon | A mathematical theory of communication[END_REF] and the Maximum Entropy Principle [see [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. Once the probability models of the geometrical parameters are determined, different city realizations corresponding to a given parameter u can be generated.

Step 2: to use the city realizations generated in Step 1 and to perform simulations with the reference model in order to obtain statistical information on the model output

P exp r (ω).
Step 3: to construct the probability models of the stochastic NPE model random parameters Γ, Λ and Θ. These models are determined using the same methods as in Step 1 (Information Theory and the Maximum Entropy Principle) and depend on parameter w.

Step 4: using the outputs from the reference model P exp r (ω) (Step 2), to identify parameter w corresponding to the given parameter u previously fixed. This identification is done by solving an inverse stochastic problem: the "distance" between P r (ω; w)

and P exp r (ω) is minimized so that the optimal parameter w opt is obtained. Once w opt is determined, the stochastic NPE model can be used to study nonlinear wave propagation over urban cities.

The chapter is organized as follows. Section 6.3 presents the method to construct the probability models of a urban city (Step 1). The linear reference propagation model is detailed in Section 6.4 (Step 2) while the mean NPE model for sound propagation over urban cities and its parameters are presented in Section 6.5. In Section 6.6, the probability models of the NPE model random variables are determined (Step 3). The stochastic NPE model is then presented in Section 6.7. Section 6.8 deals with a hybrid method based on the mean-square and the maximum likelihood methods to solve the inverse stochastic problem and to identify the parameter w opt (Step 4). Finally, Section 6.9 presents an application and a validation procedure for the computational model of sound propagation over urban cities. Chapter summary and conclusions are given in Section 6.10.

Prior probabilistic model of geometrical parameters of a urban city

A two-dimensional cross-section of a urban city is considered. It is composed of n buildings of rectangular shape parametrized by a set of three parameters. For a given building B i , its height and width are noted h i and w i , and the distance between two consecutive buildings Let X i be the positive-valued real random variable representing either H i , either W i , or D i .

Let p X i be the probability density function of X i and E be the mathematical expectation.

Let m X i and δ X i = σ X i /m X i be its mean value and its coefficient of variation, in which

σ 2 X i = m 2 -m 2 X i
is the variance and m 2 the second-order moment. One has

m X i =E {X i } = +∞ 0 x i p X i dx i , (6.1a) m 2 =E X 2 i = +∞ 0 x 2 i p X i dx i . (6.1b)
In order to construct the probability distributions of random variables families {H 1 , . . . , H n }, {W 1 , . . . , W n } and {D 1 , . . . , D n-1 } Information Theory [START_REF] Shannon | A mathematical theory of communication[END_REF] and the Maximum Entropy Principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] are used. One then has to define the available information for these random variables, which is the following: For all i, denoting by X i either H i , either W i or D i one has:

(1) X i is a random variable with values in ] 0, +∞ [ ,

(2) the mean value m X i = E {X i } is given and is equal to m X , independent of i, The urban city parameter u is given. The reference, linear propagation model is presented on the left, while the stochastic, simplified NPE model is shown on the right.

Simplified wave propagation model

h 1 w 1 d 1 h 2 w 2 d i h n w n d n-1 w i h i B i B 1 B 2 B i+1 B n-1 B n Figure 6
.2: Sketch of a representative urban city.

(3) the inverse X -1 i of X i is a second-order random variable that is satisfied if

E {log (X i )} = c i , with |c i | < +∞ . (6.2)
The constraint defined by Eq. ( 6.2) introduces an arbitrary constant c i which does not have any physical meaning and which is then rewritten as a function of the coefficient of variation of X i which is δ X , independent of i.

The use of the Maximum Entropy Principle yields [see for example [START_REF] Soize | Random matrix theory for modeling uncertainties in computational mechanics[END_REF] p H 1 ,...,Hn,W 1 ,...,Wn,D 1 ,...,D n-1 (h 1 , . . . , h n , w 1 , . . . , w n , d 1 , . . . , d n-1 ) =

n i=1 p H i (h i ) n i=1 p W i (W i ) n-1 i=1 p D i (d i ) , (6.3)
with

p H i (h i ) = p X (h i ; m H , δ H ) , (6.4a
)

p W i (w i ) = p X (w i ; m W , δ W ) , (6.4b
)

p D i (d i ) = p X (d i ; m D , δ D ) , (6.4c)
where

p X (x ; m X , δ X ) = ✶ ]0,+∞[ (x) 1 m X 1 δ 2 X 1 Γ 1/δ 2 X x m X 1 δ 2 X -1 exp - x δ 2 X m X . (6.5)
In Eq. (6.5), ✶ ]0,+∞[ (x) = 1 if x > 0 and 0 otherwise, and Γ (z) is the gamma function of argument z defined by

Γ (z) = +∞ 0 t z-1 e -t dt .
(6.6)

Since no available information concerning the statistical dependence between the families of random variables {H i } i , {W i } i and {D i } i is used, the Maximum Entropy Principle yields independence of all the random variables as a result, as can be seen in Eq. ( 6.3).

The following vector u of the parameters of the probabilistic model is introduced

u = (m H , m W , m D , δ H , δ W , δ D ) . ( 6

Reference model: linear propagation over urban cities

This section presents the stochastic linear propagation model used to obtain reference solutions of the problem. The problem setting (source and receiver, frequency range, etc.) is first given in Section 6.4.1 and the stochastic propagation model itself is then described in Section 6.4.2.

Setting the problem

Output P exp r (ω) from the reference model is analyzed in the frequency-domain on a frequency band defined by B = ] 0, ω max ], in which ω max is such that B is a low frequency band. Hence, the buildings surfaces are assumed to be acoustically rigid and reflections are assumed to be specular. The source and receiver are placed above the urban city layer and the source is placed sufficiently far from the first building B 1 to assume that the wave impinging on this first building is a plane wave. The receiver is placed at the same altitude than the source and at the right of the last building B n . The quantity under interest is the normalized spectrum L exp (ω) at the receiver such that L exp (ω) = 10 log 10 P exp r (ω) p free (ω) 2 , (6.8)

where p free (ω) and P exp r (ω) denote free field pressure and the pressure in the presence of the urban city at the receiver, respectively.

Stochastic Boundary Element Method (BEM) solver for constructing a reference solution

The model used to construct the reference solution is an implementation of the Boundary Element Method (BEM) [see [START_REF] Ciskowski | Boundary Element Methods in Acoustics[END_REF][START_REF] Jean | A variational approach for the study of outdoor sound propagation and application to railway noise[END_REF], in which the city geometry is explicitly entered into the computational model. For a fixed parameter u, the outputs L exp (ω) are calculated using the Monte Carlo method [see [START_REF] Hammersley | Monte Carlo Methods[END_REF] with ν exp independent realizations η 1 , . . . , η νexp of urban cities generated with the probability models presented in Section 6.3. The realization of the experimental observation for the urban city ν p is L exp ω, η νp , ω ∈ B .

Mean nonlinear parabolic propagation model for sound propagation over urban cities

This section describes the mean parabolic propagation model that will be later used, associated with a probabilistic model of uncertainties, to model sound propagation over urban environments. Section 6.5.1 formally defines the NPE propagation model and Section 6.5.2 explains the choice of the propagation model parameters and details their algebraic properties. multiple media is written as

∂ ⋆ t r a + ∂ x c 1 r a + β c 0 2 r 2 a + c 0 2 ∂ 2 z r a dx = g (x a , t)
in Ω a , (6.9a)

∂ ⋆ t r u + ∂ x c 0 (µ -1) r u + β µ c 0 2 r 2 u + µc 0 2 ∂ 2 z r u dx + αr u = 0 in Ω u , (6.9b) ∂ z r a + ǫ ∂ z r a dx = γ ∂ z r u on Γ , (6.9c)
r a = 0, r u = 0 on Γ r , (6.9d) r a , r u satisfy Sommerfeld radiation condition on Γ l , (6.9e)

∂ z r u = 0 on Γ u , (6.9f)
r a satisfies Sommerfeld radiation condition on Γ a , (6.9g)

where g (x a , t) is an external pressure field applied in the atmospheric layer. Eqs. (6.9a, 6.9b) handle propagation in the atmospheric and urban layers. Eq. (6.9c) is the interface condition between the domains and Eq. (6.9d) is a standard boundary condition applied on the boundary to the right of the moving window Γ r . Eq. (6.9e) states that the wave must not reflect at the left boundary of the domain, and hence the acoustic field must satisfy Sommerfeld radiation condition on Γ l . It is supposed that the bottom of the domain is perfectly rigid, hence ∂ z r u = 0 on Γ u (see Eq. (6.9f)). Since the propagation problem is infinite in the +z direction the acoustic field must satisfy Sommerfeld radiation condition on Γ a (Eq. (6.9g)). For the and the time derivative is changed accordingly, such that

∂ ⋆ t -→ λ ∂ ⋆ t .
(6.12)

The change of amplitude of the reflected wave is accounted for using Eq. (6.9c). Parameter γ can be used to tune the amount of reflected and transmitted waves. The time delay being already taken into account with the time stretching, there is no need to keep parameter ǫ in Eq. (6.9c). Parameter ǫ is set to ǫ = 0. One advantage of using such two parameters instead of four in the mean model, is that these two parameters γ and λ do not introduce coupling effects between wave amplitude and time delay. This means that changing parameter λ does not change the reflected wave amplitude, and changing parameter γ does not change the time delay of the reflected wave. This property will later help designing the probabilistic model of these two parameters.

Note that parameter λ cannot take the value 0, which would correspond to an infinite speed of sound. Parameter γ is used to tune the amount of reflected wave in the atmospheric layer. Setting γ = 0 yields ∂ z r a = 0 on coupling interface Γ, which is the condition for perfect reflection. Setting γ = 1 yields ∂ z r a = ∂ z r u on the coupling interface Γ, which is the condition for perfect transmission; parameter γ thus belongs to [0, 1]. Introducing underlined quantities related to the mean model, the mean NPE model for propagation in multiple media with the two above parameters is rewritten as

λ ∂ ⋆ t r a + ∂ x c 1 r a + β c 0 2 r a 2 + c 0 2 ∂ 2 z r a dx = g (x a , t)
in Ω a , (6.13a)

λ ∂ ⋆ t r u + ∂ x β c 0 2 r u 2 + c 0 2 ∂ 2 z r u dx = 0 in Ω u , (6.13b) ∂ z r a = γ ∂ z r u on ∂Ω , (6.13c)
r a = 0, r u = 0 on Γ r , (6.13d) r a and r u satisfy Sommerfeld radiation condition on Γ l , (6.13e) r u satisfies Sommerfeld radiation condition on Γ u , (6.13f) r a satisfies Sommerfeld radiation condition on Γ a . (6.13g) Furthermore, introducing a third parameter θ in order to control the output of the mean model, the relative Sound Pressure Level (SPL) at the receiver is given by L (ω) = 10 log 10 ρ 0 c 2 0 r a (x r , z r , ω + θ) p free (x r , z r , ω + θ) (6.16) in which the superscript 0 refers to the nominal value of s. The value of the initial parameter s 0 is chosen so that the nominal mean model simulates a wave reflecting on a plane which is an acoustically rigid surface.

s 0 = (γ 0 , λ 0 , θ 0 ) = (0 , 1 , 0) ,
6.6 Construction of the probabilistic model of random variables Γ, Λ and Θ Let Γ, Λ and Θ be the random variables associated with the mean model parameters γ, λ and θ, respectively. Information Theory [START_REF] Shannon | A mathematical theory of communication[END_REF] and the Maximum Entropy Principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] are used to determine the random variables probability distributions (see Section 6.3).

Construction of the probability distribution of random variable Γ

The available information for random variable Γ is the following.

◮ Γ is a random variable with values in [0, 1] (see Section 6.5.2),

◮ Its mean value m Γ = E {Γ} is given, ◮ and its coefficient of variation δ Γ is given.

It should be noted that the upper bound γ = 1 corresponds to perfect reflection at the interface and that the lower bound γ = 0 corresponds to perfect transmission. Since the neighborhoods of these two bounds can be reached with a non-zero probability, it is not necessary to introduce an available information related to the behavior of the probability distribution in the neighborhood of these two bounds. With such an available information, the Maximum Entropy Principle yields (6.17) where the constants µ 0 , µ 1 and µ 2 depend on m Γ and δ Γ and are the solutions of the equations

p Γ (γ) = ✶ [0,1] (γ) e -µ 0 -γ µ 1 -γ 2 µ 2 ,
1 0 γ e -γ µ 1 -γ 2 µ 2 dγ -m Γ 1 0 e -γ µ 1 -γ 2 µ 2 dγ = 0 , (6.18a) 1 0 γ 2 e -γ µ 1 -γ 2 µ 2 dγ -m 2 Γ + σ 2 Γ 1 0 e -γ µ 1 -γ 2 µ 2 dγ = 0 , (6.18b) 1 0 e -γ µ 1 -γ 2 µ 2 dγ -e µ 0 = 0 . (6.18c)
In the equations above, integrals are numerically evaluated with the Monte Carlo method and the equations are solved using a nonlinear multidimensional root finding algorithm.

Construction of the probability distribution of random variable Λ

Taking into account the algebraic properties given in Section 6.5.2, random variable Λ is with values in ] 0, +∞ [ . Since λ cannot take the value 0 which would correspond to an infinite speed of sound, it is necessary to write that the probability distribution goes sufficiently fast to zero when λ goes to zero with superior values. Such a property is satisfied in introducing the following condition, (6.19) which implies that the inverse Λ -1 of Λ is a second-order random variable. Therefore, the available information for random variable Λ is the following.

E{log Λ} = c , |c| < +∞ ,
◮ Λ is a random variable with values in ]0, +∞[, ◮ its mean value m Λ = E{Λ} is given, ◮ and Eq. (6.19) is satisfied.

With the constraints defined by the above available information, the use of the Maximum Entropy Principle yields (6.20) in which δ Λ is the coefficient of variation of Λ.

p Λ (λ) = ✶ ]0,+∞[ (λ) 1 m Λ 1 δ 2 Λ 1 Γ 1/δ 2 Λ λ m Λ 1 δ 2 Λ -1 exp - λ δ 2 Λ m Λ ,

Construction of the probability distribution of random variable Θ

The available information relative to random variable Θ is the same as the available information defined in Section 6.6.2 for random variable Λ. Consequently, the probability density function of random variable Θ is written as

p Θ (θ) = ✶ ]0,+∞[ (θ) 1 m Θ 1 δ 2 Θ 1 Γ 1/δ 2 Θ θ m Θ 1 δ 2 Θ -1 exp - θ δ 2 Θ m Θ , (6.21)
in which m Θ and δ Θ are the mean value and the coefficient of variation of random variable Θ.

Stochastic Nonlinear Parabolic Equation model for highamplitude wave propagation over urban cities

The following vector w of the parameters of the probabilistic models of random variables Γ, Λ and Θ is introduced

w = (m Γ , m Λ , m Θ , σ Γ , σ Λ , σ Θ ) , (6.22)
in which

σ Γ = m Γ δ Γ , (6.23a) σ Λ = m Λ δ Λ , (6.23b) σ Θ = m Θ δ Θ . (6.23c)
Parameter w belongs to the admissible set W = ( ] 0, +∞ [ ) 6 . The stochastic model for nonlinear sound propagation over urban cities is defined by

Λ ∂ ⋆ t R a + ∂ x c 1 R a + β c 0 2 R 2 a + c 0 2 ∂ 2 z R a dx = g in Ω a , (6.24a) Λ ∂ ⋆ t R u + ∂ x β c 0 2 R 2 u + c 0 2 ∂ 2 z R u dx = 0 in Ω u , (6.24b) ∂ z R a = Γ ∂ z R u on Γ , (6.24c) R a = 0, R u = 0 on Γ r , (
L (ω, w) = 10 log 10   ρ 0 c 2 0 R a (x r , z r , ω + Θ) p free (x r , z r , ω + Θ) 2   .
(6.25)

In Eqs. (6.24, 6.25), the probability distributions of random variables Γ, Λ and Θ depend on w, which in turn depends on u, the parameter that describes the urban city geometry.

Hence, to complete the construction of the computational model and to obtain observations of L (ω, w), one has to express parameter w as a function of parameter u (it should be noted that no explicit expression can be constructed but the corresponding mapping will be numerically constructed).

6.8 Identification of parameter w of the stochastic model

Identification strategy

The identification of the parameter w is performed by solving an inverse stochastic problem.

Two methods are generally used to solve such problems: the mean-square method [see [START_REF] Spall | Introduction to Stochastic Search and Optimization[END_REF][START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF]] and the maximum likelihood method [see [START_REF] Soize | Probabilistic model identification of uncertainties in computational models for dynamical systems and experimental validation[END_REF].

Both techniques were previously tested with different urban environments parameters u; the conclusions of this comparative study are the following.

(i) Since the mean-square method introduces a mean-square distance between the experimental data and the random response of the stochastic model, this type of method is equivalent to a minimization of the sum of the variance for the stochastic model response with the bias between the experimental mean value and the mean value of the random response. If these two mean values are significantly different, the bias can only be reduced in increasing the variance of the model. In this case, the distance between the experimental mean value and the mean value of the random response of the stochastic model is effectively reduced but in counterpart, the confidence region of the random response increases.

(ii) For the application analyzed it has been seen that, for certain values of parameter u, the maximum likelihood method under-estimates the width of the confidence region which means that a significant number of experimental paths cross the upper and the lower envelopes of the confidence region.

In order to provide an accurate and robust identification method for all values of parameter u, a hybrid method is used; the identification of the optimal parameter w opt is done in two successive steps.

Step 1. The mean model is first updated with experimental data. This step consists in finding the optimal mean model parameter s opt that minimizes the norm between the mean value of the reference model outputs and the mean model output. As only the mean value of the experimental observations is used, this first step only allows us to identify the parameter that optimizes the mean response of the stochastic model.

Step 2. Next, to take into account the dispersion of the experimental observations, a multiobjective optimization problem is solved: both the mean-square and the maximum likelihood methods are used concurrently. This multi-objective optimization problem aims at finding the optimal parameter w opt that maximizes the log-likelihood function between the experimental observations and the stochastic model outputs, while minimizing the areas where experimental observations do not belong to the confidence regions of the stochastic model. This problem is solved using an evolutionary algorithm [see [START_REF] Deb | A Fast Elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF][START_REF] Srinivas | Multiobjective optimization using nondominated sorting in References genetic algorithms[END_REF].

Step 1 allows us to obtain a rough approximation of the mean values m Γ , m Λ and m Θ which are the three first components of vector w and hence helps initializing the optimization problem in Step 2. Section 6.8.2 presents the procedure to update the mean model with experimental data (Step 1). The mean-square method and the maximum likelihood method are presented in Sections 6.8.3 and 6.8.4.

Updating the mean model with experimental data

The observation from the mean model depending on the choice of parameter s it is rewritten as L s (ω). The performance level of the nominal model can be measured in estimating the norm 6.8.3 Identification of the stochastic model optimal parameter with the mean-square method

||E {L exp } -L s || B = ω∈B |E {L exp (ω)} -L s (ω)| 2 dω 1/2 , ( 6
This section presents the mean-square method [see [START_REF] Spall | Introduction to Stochastic Search and Optimization[END_REF][START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF] with non-differentiable objective function [see [START_REF] Soize | Probabilistic model identification of uncertainties in computational models for dynamical systems and experimental validation[END_REF] for the identification of the parameter w opt of the stochastic model.

The objective function is defined by writing that the η νexp experimental observations {L exp (ω ; η j ) , ω ∈ B, j = 1, . . . , ν exp } (6.28) must belong to the confidence region of the stochastic model with a probability level P c fixed in ] 0, 1 [ . Before giving the formal definition of the mean-square method one must address the problem of the construction of confidence regions.

Construction of confidence regions. The confidence region is constructed by using the quantiles. Let F L(ω) be the cumulative distribution function of random variable L (ω), such that F L(ω) (ℓ) = Proba {L (ω) ≤ ℓ}. For 0 < p < 1, the p th quantile (or fractile) of F L(ω) is defined as

ζ (p ; ω) = inf ℓ F L(ω) (ℓ) ≥ p . (6.29)
Then, the upper and lower envelopes ℓ + (ω) and ℓ -(ω) of the confidence region are given by

ℓ + (ω) = ζ 1 + P c 2 ; ω , (6.30a) ℓ -(ω) = ζ 1 -P c 2 ; ω . (6.30b)
The estimation of ℓ + (ω) and ℓ -(ω) is performed by using the sample quantile [see [START_REF] Serfling | Approximation Theorems of Mathematical Statistics[END_REF]. Consider ν independent realizations of the random variable L (ω) noted ℓ 1 (ω) = L (ω ; η 1 ) , . . . , ℓ ν (ω) = L (ω ; η ν ) and let ℓ 1 (ω) < . . . < ℓ ν (ω) be the ordered statistics associated with ℓ 1 (ω) , . . . , ℓ ν (ω). One has the following estimations for the upper and lower 6.8.4 Identification of the stochastic model optimal parameter with the maximum likelihood method and statistical reduction of information

This section deals with the maximum likelihood method to identify the optimal parameter w opt . Let {ω 1 , . . . , ω m } ⊂ B be a sampling of frequency band B and let w→L be the log-likelihood function from W into |, defined by

L (w) = νexp j=1
log 10 p (L exp (ω 1 ; η j ) , . . . , L exp (ω m ; η j ) ; w) , (6.36) in which p (L exp (ω 1 ; η j ) , . . . , L exp (ω m ; η j ) ; w) is the joint probability density function of random variables L (ω 1 , w) , . . . , L (ω m , w) for the values L exp (ω 1 ; η j ), . . ., L exp (ω m ; η j ).

The maximum likelihood method [see [START_REF] Spall | Introduction to Stochastic Search and Optimization[END_REF] consists in finding w opt as the solution of the following optimization problem,

w opt = arg max w∈W L (w) . (6.37)
The standard method requires a direct evaluation of the joint probability density function appearing in Eq. ( 6.36) which involves a very important computational effort. In order to decrease this effort the method of statistical reduction of information introduced by Soize et al. [2008] is used. The idea is to proceed to a statistical reduction of information using a principal component analysis (see for instance [START_REF] Jolliffe | Principal Component Analysis[END_REF]) and then to use the maximum likelihood method in the space of the uncorrelated random variables related to the reduced statistical information.

Statistical reduction of information. For all w fixed in W, let be L (w) = (L (ω 1 , w) , . . . , L (ω m , w)). Let m (w) = E {L (w)} be its mean value and let C L(w) be its (m × m) covariance matrix defined by for which the first largest q ≤ m positive eigenvalues are e 1 (w) ≥ e 2 (w) ≥ . . . ≥ e q (w)

C L(w) = E (L (w) -m (w)) (L (w) -m (w)) T , ( 6
and the associated eigenvectors are x 1 (w) , . . . , x q (w), in | m . The approximation L q (w) of L (w) is written as 6.40) in which Y 1 (w) , . . . , Y q (w) are q real-valued random variables such that, for all α = 1, . . . , q Y α (w) = 1 e α (w)

L q (w) = m (w) + q α=1 e α (w) Y α (w) x α (w) , (
x α (w) T (L (w)m (w)) . (6.41)

It can easily be proved that Y (w) = (Y 1 (w) , . . . , Y q (w)) is a second-order random variable such that, for all α and β in {1 . . . q} E {Y α (w)} = 0 , (6.42a) (6.42b) which means that the centered random variables Y 1 (w) , . . . , Y q (w) are uncorrelated. The order q of the statistical reduction is calculated in order to get an approximation with a given accuracy ǫ which has to be chosen such that

E {Y α (w) Y β (w)} = δ αβ ,
max w∈W 1 - q α=1 e α (w) tr [C L (w)] ≤ ǫ . (6.43)
From Eq. (6.41) one can deduce that random variables Y exp α (w ; η j ) associated with the experimental realization η j are given, for all α = 1, . . . , q and j = 1, . . . , ν exp by

Y exp α (w ; η j ) = 1 e α (w) x α (w) T (L exp (η j ) -m (w)) . (6.44)
Finally, taking into account that the random variables Y 1 (w) , . . . , Y q (w) are mutually independent, one introduces the following approximation for the reduced log-likelihood function 6.45) where p Yα(w) (y ; w) is the probability density function of the real-valued random variable Y α (w). The optimal value w opt of w is then given by w opt = arg max w∈W L red (w) , (6.46) and the objective function for the maximum likelihood method is -L red .

L red , L red (w) = νexp j=1 q α=1 log 10 p Yα(w) (Y exp α (w ; η j ) ; w) , (
Eqs. (6.47, 6.48).

Step 2, non-dominated sorting and crowding distance assignment. Once the families of individuals are generated, the population is sorted against the multi-objective function F (w) defined by (see Sections 6.8.3 and 6.8.4)

F (w) = -L red (w) , J (w) . (6.50)
Each individual is assigned a rank: individuals in the first front (rank 1) dominate individuals in the second front (rank 2) and so on. In addition to the rank value a crowding distance is calculated for each individual. The crowding distance is a measure of how close is an individual from its neighbors: large crowding distances result in better diversity in the population.

Step 3, selection process. Parents are selected by using binary tournament selection based on their rank and crowding distance. Individuals with smaller rank will be selected first and between individuals of the same rank, individuals with large crowding distances are favored in order to diversify the population.

Step 4, genetic operators. Genetic algorithms use a simulated binary crossover operator [see [START_REF] Beyer | On self-adaptive features in real-parameter evolutionary algorithm[END_REF][START_REF] Deb | Simulated binary crossover for continuous search space[END_REF] for offspring generation. Two parents p 1 and p 2 (p 1 and p 2 different) are selected randomly within the population (Step 3) and children are generated according to

c 1,k = 1 2 [(1 -β k ) p 1,k + (1 + β k ) p 2,k ] , (6.51a) c 2,k = 1 2 [(1 + β k ) p 1,k + (1 -β k ) p 2,k ] , (6.51b)
in which p i,k denotes the k th component of the i th parent, c i,k denotes the k th component of the i th children and β k is a sample from a random number having the probability distribution

p (β) = 1 2 (η c + 1) β ηc , if 0 ≤ β ≤ 1 , (6.52a) p (β) = 1 2 (η c + 1) 1 β ηc+2 , if β > 1 .
(6.52b)

The parameter η c determines how well spread are the children from the parents. If the parents p 1 and p 2 are identical, i.e. p 1,k ≡ p 2,k = p k , polynomial mutation is used instead; the children component c k is then generated with (6.53) where p u k and p l k are the upper and lower bound of the parent component, respectively, and δ k is a dispersion parameter defined by (6.54b) in which r k is a sample from a uniform distribution in [0, 1] and η m is a mutation distribution index, controlling how well spread is the children from the parent.

c k = p k + p u k -p l k δ k ,
δ k = (2r k ) 1 ηm+1 -1, if r k < 0.5 , (6.54a) δ k = 1 -(2 (1 -r k )) 1 ηm+1 , if r k ≥ 0.5 ,
Step 5, recombination and selection. The offspring population is combined with the current generation and selection is performed based on non-domination to set the n i individuals of the next generation. The genetic operators, recombination and selection steps are then repeated until the number of generations wanted is reached and finally, the NSGA-II algorithm yields a final generation of individuals that best minimizes the objective function F (w).

Application and experimental validation

It is proposed in this section to use the stochastic parabolic propagation model to study sound propagation over given urban environments.

Summary of previous sections -Stochastic model validation procedure

In order to construct and validate the stochastic propagation model one has

Step 1. to choose a parameter u = (m H , m W , m D , δ H , δ W , δ D ) that describes the urban city geometry studied. With this parameter the ν exp urban environment realizations can be generated with the probability distributions given in Section 6.3 (see Section 6.9.2),

Step 2. using the ν exp realizations of the probabilistic model of the given city, and using the stochastic reference model (see Section 6.4), to produce the family {L exp (ω ; η j )} j (see Section 6.9.2),

Step 3. once the experimental realizations are obtained, to identify the optimal parameter w opt of the stochastic model with the help of the method described in Section 6.8 (see Section 6.9.3), seen that as the population evolves the values of the objective functions -L red (w) and J (w)

decrease.

For the case where the urban environment is parametrized by u 1 one can see that the first generation of individuals (after one selection/mutation/recombination process) shows very high dispersion on the objectives space and low performances (about half the individuals are out of the range of the figure). The parameters chosen to initialize the population (see Eqs. (6.47, 6.48)) are not adapted to this urban environment morphology. A ten-generation evolutionary cycle allows us to obtain a population with acceptable performances: the loglikelihood is reduced to about 23 000 and the mean-square norm is below 200. The population at the final generation is diversified: the lowest log-likelihood is 20 020 (with a mean-square norm close to 300) and the lowest mean-square norm is 84 (with a log-likelihood near 24 000).

The log-likelihood increases at the expense of the mean-square norm, and inversely. To select the final individual one has thus to chose between the two objectives which one should be favored, or accept a trade-off between the log-likelihood and the mean-square norm. Note that despite the fact that the initial parameters are not adapted to the city morphology studied the evolutionary algorithm achieves good performances.

For the case where the urban environment is parametrized by u 2 one can see that the initial population shows better results: mean-square norm is below 10 for 42 individuals. Indeed as experimental observations from parameter u 2 exhibit high dispersion and a low amplitude range (mean value between -6 dB and +6 dB), and as the initial population is uniformly distributed on W red , it is more likely that the experimental realizations belong to the confidence regions of the stochastic model, hence the low mean-square norm values.

The chosen criteria to select the final individuals is the mean-square norm: individuals that show the lowest mean-square norm are selected first and within the remaining individuals, the individual that exhibits the highest log-likelihood value is selected. This way, the least possible under-estimation of the experimental observations is achieved.

The corresponding optimal parameters w opt 1 and w opt 2 are The stochastic propagation model contains three entities: two nonlinear parabolic propagation equations for the air and the porous layer and a boundary interface condition used to couple these two domains. The model is then adapted for the application of sound propagation over urban environments. Three parameters that allow the output of the model to be controlled are selected and their algebraic properties are described. With this information, the Maximum Entropy Principle is used to determine their probability distributions. The NPE model is numerically solved using the finite-difference method. Specifically, the operator splitting method is used to solve for the different operators with specialized algorithms.

w opt 1 = m opt Γ , m opt Λ , m opt Θ , σ opt Γ , σ opt Λ , σ opt Θ = (0.
Nonlinear terms are solved using the Flux Corrected Transport (FCT) algorithm, which allows discontinuities to be stably propagated. Linear terms are handled with the Crank-Nicolson method, yielding tridiagonal systems of equations, solved with a fast algorithm (the Thomas algorithm).

The validation procedure is separated in successive tasks. First, the NPE model is compared to analytical solutions to assess the accuracy of nonlinear effects calculations. The implemented NPE has been proved to accurately simulate wave steepening and harmonics generation with a relative error lower than 5 % until the third harmonic. The computational model has then Three parameters that allow the output of the stochastic model to be controlled are then selected and their algebraic properties are described. With this information, the Maximum Entropy Principle is used to determine their probability distributions. Once the probability distributions are determined the stochastic model parameter is identified solving an inverse stochastic problem. This problem involves two different methods: the mean-square method and the maximum likelihood method. In order to reduce the computational effort associ-Table A.1: Overpressure values and corresponding damage to structures from [START_REF] Lannoy | Analyse des explosions air-hydrocarbure en milieu libre : Etudes déterministes et probabilistes du scénario d'accident[END_REF] and [START_REF] Lees | Loss Prevention in the Process Industries[END_REF].

Overpressure (kPa) Damage The waveform is given by

p (t) = p 0 1 -t/t + e αt/t + , (B.1)
where p (t) is the overpressure at time t, p 0 is the peak overpressure at time t = 0 (t = 0 at the shock arrival time) and t + is the positive phase duration. For a chemical explosion, the ratio of the overpressure to the ambient pressure p a can be calculated with

p 0 p a = 808 1 + Z 4.5 2 1 + Z 0.048 2 1 + Z 0.32 2 1 + Z 1.35 2 , (B.2)
in which Z is a reduced distance which units is m/kg 1/3 . The positive phase duration t + is calculated with

t + W = 980 1 + Z 0.54 10 1 + Z 0.02 3 1 + Z 0.74 6 1 + Z 16.9 2 . (B.
3)

The arrival time t a (t=0 in Eq. (B.1)) is defined by

t a = 1 c 0 r rc 1 1 + 6p 0 7pa dr , (B.4)
where r c is the charge radius and r is the distance from the source. The shape parameter α

C

Numerical solution to Euler Equations for propagation over a porous ground layer

This appendix briefly presents the numerical solution and the associated algorithms used to solve the Euler equations in Section 5.4.2.

The computational domain is composed of an air and a ground layer. For the air layer and in a two-dimensional Cartesian coordinate system the constitutive equations are such that

∂ t ρ T + ∂ x (ρ T u) + ∂ z (ρ T w) = 0 , (C.1a) ∂ t (ρ T u) + ∂ x ρ T u 2 + ∂ z (ρ T uw) = -∂ x p T , (C.1b) ∂ t (ρ T w) + ∂ x (ρ T uw) + ∂ z ρ T w 2 = -∂ z p T , (C.1c) ∂ t (ρ T e 0 ) + ∂ x (ρ T ue 0 ) + ∂ z (ρ T we 0 ) = -∂ x (p T u) -∂ z (p T w) , (C.1d)
where e 0 is the energy per unit mass. Within the ground layer momentum conservation equations write

Φ∂ t (ρ T u) + ∂ x p T + Φρ T u 2 + ∂ z (Φρ T uw) + σ 0 Ω 0 (1 + ξ |u|) u = 0 , (C.2a) Φ∂ t (ρ T w) + ∂ z p T + Φρ T w 2 + ∂ x (Φρ T uw) + σ 0 Ω 0 (1 + ξ |w|) w = 0 . (C.2b)
The energy equation Eq. (C.3) and the ideal gas law Eq. (C.4) close the equation system.

They are defined by

ρ T e 0 = ρ T c v T + ρ T |V| 2 2 , (C.3)
Although the combination of WENO and Runge-Kutta (RK) schemes allow discontinuities to be stably propagated, it is unable to propagate waves of infinite slope. A shock smearing will occur where the slope is too steep, resulting in small deviations from physical solutions for very high amplitude waves.

D Nonlinear Parabolic Equation model derivation without the use of the perturbation expansion method

The derivation of the constitutive equations that yields the Nonlinear Parabolic Equation (NPE) model is rewritten in a Cartesian, two-dimensional coordinate system. The procedure presented here does not use the perturbation expansion method as in Section 2.1, but rather makes use of specific assumptions about the flow movement. The aim of this section is to provide the reader a different insight to the NPE model. The derivation is taken from the article by [START_REF] Caine | A tutorial on the nonlinear progressive wave equation (NPE) -Part 2. Derivation of the three dimensional cartesian version without use of perturbation expansions[END_REF].

The Cartesian formulation of the continuity equation is The integral in Eq. (D.12) can be removed by differentiating Eq. (D.12) with respect to time, such that

∂ t ρ T + ∂ x (ρ T u) + ∂ z (
∂ 2 t ρ ′ -∂ 2 x p T + ρ T u 2 + c 0 ∂ 2 z ∂ t ρ ′ dx = 0 . (D.13)
The term ∂ t ρ ′ can be simplified by using a one-dimensional linearized momentum equation since one has only to maintain first-order accuracy in this term. The linearized momentum equation is Eq. (D.13) becomes

∂ t ρ ′ = -c 0 ∂ x ρ ′ , (D.
∂ 2 t ρ ′ -∂ 2 x p T + ρ T u 2 -c 2 0 ∂ 2 z ρ ′ = 0 (D.16)
Note that the classical (linear) wave equation can easily be derived from Eq. (D.16). If the ambient medium pressure is constant then p T = p 0 + p ′ = p 0 + c 2 0 ρ ′ in the linear case. If high order terms are neglected one gets ∂ 2 t ρ ′ -c 2 0 ∇ 2 ρ ′ = 0 , (D.17) which is the well-known wave equation.

A one-way propagation is now assumed, i.e. there is no backward propagation. Thanks to this new assumption Eq. (D.16) can be rewritten in a wave-following coordinate system. A "moving window" operator is introduced, such that

D t = ∂ t + c 0 ∂ x . (D.18)
It can be proved that D 2 t ρ ′ ≃ 0 [see Caine & West, 1995, and Section 2.1.2]; one then obtains (D.19) and Eq. (D.16) becomes One hence obtains the following expression for the total pressure p T .

∂ 2 t ρ ′ = -2c 0 ∂ x D t ρ ′ -c 2 0 ∂ 2 x ρ ′ ,
∂ x D t ρ ′ + 1 2c 0 ∂ 2 x p T + ρ T u 2 + c 2 0 ρ ′ + c 0 2 ∂ 2 z ρ ′ = 0 (D.
p T ≃ p 0 + c 2 ρ ′ + 1 2 (γ -1) c 2 ρ 0 ρ ′2 (D.24a) u ≃ c ρ ′ ρ 0 . (D.24b)
Carrying out the replacements in the equations above, integrating with respect to x and replacing the density perturbation ρ ′ by the dimensionless variable R= ρ ′ ρ 0 lets one finally obtain an expression for the NPE in a two-dimensional Cartesian coordinates system, defined by

D t R + ∂ x c 1 R + c 0 β 2 R 2 + c 0 2 ∂ 2 z R dx = 0 (D.25)
The equation above is the original formulation of the Nonlinear Parabolic Equation (NPE), as defined by Eq. (2.19).

to solve the system [B] u = d where [B] is tridiagonal and of size N × N such that

[B] =               α 1 β 1 γ 1 α 2 β 2 . . . . . . . . . γ N -2 α N -1 β N -1 γ N -1 α N               . (E.2)
The Thomas algorithm is decomposed in three steps:

◮ a Lower/Upper (LU) decomposition, ◮ a forward substitution, ◮ and a backward substitution.

Step 1, LU decomposition. The matrix [B] can be decomposed in a lower bidiagonal matrix and an upper bidiagonal matrix, such that

[L] =           1 l 1 1 . . . . . . l N -1 1           , (E.3) and [R] =           m 1 r 1 . . . . . . m N -1 r N -1 m N          
.

(E.4) One can note that r i = β i for all i; the coefficients m i and l i can be obtained as follows.

m 1 := α 1 For i = 1, 2, • • • , N do:

l i := γ i /m i m i+1 := α i+1 -l i β i (E.5)
Step 2, forward substitution. The system [L] y = d is solved. One has

y 1 := d 1
For i = 2, • • • , N do:

y i := d i -l i-1 y i-1 (E.6)
Step 3, backward substitution. The system [R] u = y is solved. One has

F

The Boundary Element Method

The Boundary Element Method (BEM) is a numerical method that arose from the theory of Boundary Integral Equations (BIE) during the sixties. The reader may refer to the book by [START_REF] Ciskowski | Boundary Element Methods in Acoustics[END_REF] for a complete description of the BEM. This method has been developed alternatively to another family of numerical methods, the Finite-Element Method (FEM), particularly for propagation problem involving infinite domains. Indeed, the BEM is well suited for infinite domains since only the boundary has to be discretized while for the FEM, the whole domain has to be meshed. The BEM allows the dimensionality of the problem to be reduced, since the acoustic field in the domain is calculated with its boundary radiation.

Moreover, the Sommerfeld radiation condition is satisfied by integral formulations.

The numerical implementation MICADO used in this work is based on the direct integral formulation [see [START_REF] Jean | A variational approach for the study of outdoor sound propagation and application to railway noise[END_REF]

Validation du modèle

La validation du modèle déterministe est réalisée en plusieurs étapes successives.

Le modèle NPE est d'abord comparé aux solutions analytiques pour la propagation d'ondes de grandes amplitudes en une dimension. Les résultats du modèle numérique sont comparés aux solutions de Fubini [1935]. L'erreur relative sur l'amplitude de la fréquence fondamentale et

Perspectives

Ce travail de recherche correspond à une première étape pour laquelle l'ensemble des modèles a été développé en 2D. La principale perspective est l'extension de ces modèles en 3D, sachant qu'il n'y a pas de difficultés particulières pour une telle extension. Concernant la validation expérimentale il y aurait alors lieu de réaliser des mesures sur maquettes pour obtenir les solutions de référence à la place de simulations numériques 3D avec la BEM.

L'outil développé peut être utilisé pour la cartographie dans le cadre de la prévention et de la maîtrise des risques. Une simulation complète en 2D, une fois que les solutions de référence sont obtenues, nécessite environ 1 heure de calcul sur un processeur. Cela permet de réaliser xxxii une cartographie complète en environ deux jours de calcul pour un processeur (50 directions angulaires). De plus, les informations contenues sur ces cartes peuvent être améliorées de deux manières.

1. Des cartes "conditionnelles" peuvent être obtenues. On peut imaginer réaliser plusieurs cartes du même site, en fonction des conditions météorologiques, permettant ainsi l'évaluation des plages horaires critiques.

2. Des cartes statistiques peuvent également être réalisées. Compte-tenu de la grande variabilité du système, on est en droit de se demander quelle est la variabilité répercutée sur les niveaux de pression. Le modèle développé permet de répondre à cette question, c'est-à-dire de donner un niveau de surpression maximal et un niveau de probabilité associé.

Sur le plan de l'utilisation pratique des modèles proposés dans le cadre ci-dessus, une amélioration consisterait à construire "point par point" (approche polynômiale, approche par réseau de neurones, etc.) la fonction qui permettrait d'obtenir la valeur du paramètre w opt en fonction des paramètres géometriques qui décrivent les environnements urbains. La longue et fastidieuse étape d'obtention des solutions de référence et d'identification du paramètre w opt serait alors évitée.

xxxiii

  un modèle NPE pour la propagation en milieu complexe -Aspects déterministes . . . . . . . . . . . . . . . . . xxviii 3 Développement d'un modèle probabiliste d'incertitudes -Application à la propagation stochastique non linéaire en milieu urbain . . . . . . . . xxxi 4 Synthèse du travail et perspectives . . . . . . . . . . . . . xxxii List of Publications as Author xxxv x List of Figures 1 Satellite view of the AZF factory site taken in 2009 . . . . . . . . 2 Photograph of damaged cars and buildings near the Buncefield depot taken approximately 400 meters from the main explosion site . . . . . . . 1.1 Atmospheric absorption in dB/100 m for an atmospheric pressure of 1 atmosphere, temperature of 20 ℃ and relative humidity of 70 % . . . . . 1.2 Sketch of the configuration for wave propagation above a flat surface . . 1.3 Measured excess attenuation for propagation from a point source over asphalt 1.4 Variation of temperature with altitude: examples of lapse and inversion conditions . . . . . . . . . . . . . . . . . . . . . . 1.5 Variation of average wind velocity in the vicinity of a flat ground surface . 1.6 Illustration of multiple reflections resulting from a downwind situation . . 1.7 Illustration of sound rays resulting from an upwind situation . . . . . 1.8 Illustration of combined effects of temperature gradient and wind gradient on sound propagation . . . . . . . . . . . . . . . . . . 1.9 Illustration of wave steepening . . . . . . . . . . . . . . 1.10 Illustration of harmonics generation during shock formation . . . . . 1.11 Illustration of shock coalescence . . . . . . . . . . . . . . 1.12 Time signals and Sound Pressure Levels (SPL) for an explosion of 30 T of TNT at 200 m and 1000 m from the source calculated with the Kinney-Graham (KG) model . . . . . . . . . . . . . . . . . . . . 1.13 Sketches and diagrams detailing the basic principles of each propagation model 2.1 Graphical representation of the NPE operators . . . . . . . . . 2.2 Schematic overview of the method for coupling Euler equations method with the NPE model and the PE . . . . . . . . . . . . . . . 3.1 Illustration of the coordinate transformation used to take into account the ground topography within the NPE model . . . . . . . . . . . 3.2 Graphical representation of additional terms accounting for ground topography in the GT-NPE model . . . . . . . . . . . . . . . . . xi 3.3 Graphical representation of the differential operators of the NPE model for propagation in porous ground layers . . . . . . . . . . . . 3.4 Illustration of the discretization of the boundary between the atmospheric layer and the porous ground layer . . . . . . . . . . . . . . . 3.5 Causality index of the ZK model . . . . . . . . . . . . . . 4.1 Analytical and numerical solution of the Burger equation with the MacCormack and leapfrog schemes . . . . . . . . . . . . . . . . . 4.2 Snapshots of the propagation of a blast wave in a PML . . . . . . . 4.3 Simplified flow chart of the NPE solver . . . . . . . . . . . . 5.1 Steepening of a sine wave calculated with the Fubini solution and the NPE model108 5.2 Fundamental and harmonics amplitudes evolution during wave steepening calculated with the Fubini solution and the NPE model . . . . . . . 5.3 Transmission loss (TL) calculated with a single NPE simulation and two coupled NPE/ATMOS simulations . . . . . . . . . . . . . . . . 5.4 Sound speed profile used for ATMOS and NPE calculations . . . . . 5.5 Transmission Loss (TL) map calculated with the NPE model and TL difference between NPE and ATMOS models . . . . . . . . . . . . . 5.6 Transmission Loss (TL) at three different altitudes for ATMOS and NPE simulations . . . . . . . . . . . . . . . . . . . . . 5.7 Contour plots of the field of the overpressure R taken at times 180 ms, 265 ms, 310 ms and 355 ms. . . . . . . . . . . . . . . . . . . 5.8 Sound Pressure Level (SPL) maps for the BEM and the GT-NPE model for sound propagation over an acoustically rigid hilly ground surface . . . 5.9 Sound Pressure Level (SPL) at three different altitudes for BEM and NPE simulations . . . . . . . . . . . . . . . . . . . . 5.10 Sound Pressure Level (SPL) relative to free field at the receiver for the three different ground layers, for both NPE and analytical solutions . . . . 5.11 Pressure waveform used to start the reference and NPE simulations . . . 5.12 Pressure field at time t = 33 ms for the solution from the NPE and Euler equations . . . . . . . . . . . . . . . . . . . . . 5.13 Time signals at the receiver for NPE and reference calculations for the three ground layers considered . . . . . . . . . . . . . . . . 5.14 Time signals at the receiver for a Forchheimer parameter ξ = 2.5 s.m -1 . 6.1 Sketches and diagrams detailing the basic principles of each propagation model 6.2 Sketch of a representative urban city. . . . . . . . . . . . . 6.3 Example of a city realization . . . . . . . . . . . . . . . 6.4 Sketch of the different domains with their respective boundaries of the NPE model for propagation in multiple media . . . . . . . . . . . xii 6.5 Sketch for the reference model (numerical experiment). . . . . . . 6.6 Comparisons of the reference model with the nominal model . . . . . 6.7 Comparisons of the reference model with the updated mean model . . . 6.8 Values of the objective functions for the 50 individuals at generations 1, 10 and 50 . . . . . . . . . . . . . . . . . . . . . . . 6.9 Comparisons of the reference model with the confidence region calculated with the stochastic propagation model . . . . . . . . . . and corresponding commonly accepted damage to structures . . . . . . . . . . . . . . . . . . . . . 4 3 Commonly accepted figures for direct harm to people from blast overpressure 5 4 Summary of the overall cost of the Buncefield incident, by main category . 7 1.1 TNT equivalence factors for some common explosives . . . . . . . 29 1.2 Overpressure values at 200 m and 1000 m from the source for an explosion of 30 T of TNT calculated with the Kinney-Graham (KG) model . . . . 35 1.3 Phases durations at 200 m and 1000 m from the source for an explosion of 30 T of TNT calculated with the Kinney-Graham (KG) model . . . . . 35 5.1 Ground topography characteristics . . . . . . . . . . . . . 115 5.2 Reflected wave characteristics for reference and NPE calculations. . . . 126 5.3 Reflected wave characteristics for reference and NPE calculations with Forchheimer's nonlinearities. . . . . . . . . . . . . . . . . . 127 5.4 Differences in reflected waves characteristics with and without Forchheimer's nonlinearities . . . . . . . . . . . . . . . . . . . . 128

  Cette thèse a été effectuée dans le cadre d'une Convention Industrielle de Formation par la Recherche (CIFRE) * , subventionnée par l'Association Nationale de la Recherche et de la Technologie (ANRT) † et le Centre Scientifique et Technique du Bâtiment (CSTB) ‡ . Ce travail de recherche a été effectué au laboratoire de Modélisation et Simulation Multi Echelle (MSME) § de l'Université Paris-Est (UPE) Marne-la-Vallée ¶ en étroite collaboration avec l'établissement de Grenoble du Centre Scientifique et Technique du Bâtiment (CSTB), au sein du Département Acoustique et Eclairage (DAE) , pôle acoustique environementale et urbaine.This thesis was partly financed by the ANRT (French agency for research and technology) through a partnership between Université Paris-Est (UPE) and CSTB.

  st 2001 a huge explosion occurred in the AZF fertilizer factory in Toulouse, belonging to the Grande Paroisse branch of the Total group. Three hundred tons of ammonium nitrates exploded in hangar 221. The whole factory was destroyed making a crater of depth 20 to 30 m, with a diameter of 200 m. Steel girders were found 3 km away from the explosion. Fig. (1) shows a satellite view of the AZF site eight years after the explosion.

Figure 1 :

 1 Figure 1: Satellite view of the AZF factory site taken in 2009. Eight years after the explosion, the crater is still clearly visible at the top of the picture.

Figure 2 :

 2 Figure 2: Photograph of damaged cars and buildings near the Buncefield depot taken approximately 400 meters from the main explosion site. This photograph is taken from Buncefield 2008.

◮

  and vibrational relaxation of molecules of oxygen and nitrogen. The atmospheric absorption in decibels per hundred meters due to the different physical effects is shown in Fig. (1.1). Note that the attenuation by absorption is constant for a given difference in propagation path lengths unlike geometrical spreading, where it is constant for a given ratio of propagation path lengths. Thus attenuation tends to be more and more important with increasing distance between source and receiver. Fig. (1.1) also shows the frequency dependency of the absorption coefficient. The sound attenuation is more and more important as the frequency increases; only low-frequency noise is able to propagate through large distances.

Figure 1 . 1 :

 11 Figure 1.1: Atmospheric absorption in dB/100 m for an atmospheric pressure of 1 atmosphere, temperature of 20 ℃ and relative humidity of 70 %. This figure is taken from Piercy & Embleton [1977].

Figure 1 . 2 :

 12 Figure 1.2: Sketch of the configuration for wave propagation above a flat surface. The source and receiver heights are h s and h r , respectively. The grazing angle is noted ψ.

Figure 1 . 3 :

 13 Figure 1.3: Measured excess attenuation for propagation from a point source over asphalt, with h s = 0.3 m, h r = 1.2 m. This figure is taken from Piercy & Embleton[1977].

  1.5) where T 0 is the temperature for z z 0 . The constant K t is determined by the roughness of the surface and the temperature above the boundary layer. During night the ground surface cools due to radiation in the air, a phenomenon known as an inversion condition. Examples of such temperature profiles are shown in Fig. (1.4).

Figure 1 . 4 :

 14 Figure 1.4: Variation of temperature with altitude: examples of lapse and inversion conditions. This figure is taken from Munn [1966].

Figure 1 . 5 :

 15 Figure 1.5: Variation of average wind velocity in the vicinity of a flat ground surface. This figure is taken from Piercy & Embleton [1977].

Figure 1 . 6 :

 16 Figure 1.6: Illustration of multiple reflections resulting from a downwind situation. This figure is taken from Embleton [1996].

Figure 1 . 7 :

 17 Figure 1.7: Illustration of sound rays resulting from an upwind situation. This figure is taken from Piercy & Embleton [1977].

Figure 1 . 8 :

 18 Figure 1.8: Illustration of combined effects of temperature gradient and wind gradient on sound propagation. This figure is taken from Ingård [1953].

  Fig. (1.10) shows the evolution of the first three Fourier components normalized amplitude with reduced distance.

Figure 1 .

 1 Figure 1.10: Illustration of harmonics generation during shock formation. The figure shows the evolution of the first three Fourier components of a sinusoidal wave with normalized distance σ.As the observation points moves away from the source, the energy in higher frequencies increases, while the energy contained in the initial frequency decreases, so that total energy is conserved. This figure is taken from[START_REF] Wochner | Numerical simulation of multi-dimensional acoustic propagation in air including the effects of molecular relaxation[END_REF].

  Fig. (1.11) illustrates

Figure 1 .

 1 Figure 1.11: Illustration of shock coalescence. A high-amplitude Gaussian modulated sine pulse burst is shown at three different distances. The initial waveform is given at distance x 1 . At distance x 2 , the larger shocks which have a higher velocity catch up the weaker shocks, and finally, at distance x 3 , zero-crossings disappear and shocks coalesce into a single one. This figure is taken from Wochner [2006].

[

  1990]. The KZK equation and NPE model are based on the same assumptions and are closely related formulations of a nonlinear paraxial equation. The similarity between Eq. (1.27) and Eq. (1.28) should be noted: the roles played by time and distance are reversed. The main difference between the KZK equation and the NPE model is the methods used for the treatment of initial and boundary conditions.

Fig. ( 1

 1 Fig. (1.1)) per ten kilometers at 100 Hz. One could argue that dissipation effects are not fundamental compared to other physical phenomena. However, a 6 dB reduction on the Sound Pressure Level (SPL) means an overpressure reduced by a factor 2 (assuming linear propagation). Considering the possible damage listed in Tab. (2), such a reduction of the overpressure can be of great importance.

Figure 1 . 12 :

 112 Figure 1.12: Time signals (top) and Sound Pressure Levels (SPL) (bottom) for an explosion of 30 T of TNT at 200 m (blue line) and 1000 m (green dashed line) from the source calculated with the Kinney-Graham (KG) model.

Task 1 .

 1 Development of a deterministic Nonlinear Parabolic Equation (NPE) model for highamplitude wave propagation over porous ground layers and in complex environments. Task 2. Development of a probabilistic model of uncertainties associated with the NPE model for propagation over urban environments. The paragraph below explains and details the methodology to construct such a deterministic sound propagation model and the probabilistic model of its parameters, yielding a stochastic NPE model. Methodology for the development of a deterministic NPE model for sound propagation over porous layers and in complex media. The problem with time-domain sound propagation models is that an absorbing ground surface cannot be taken into account by an impedance condition. Complex impedance models can be used with frequency-domain models, for example the PE or the Fast Field Program (FFP) models, but cannot directly be used with the NPE model. The ground layer has hence to be included into the computational system as a propagation medium. It is hence proposed to derive a Nonlinear Parabolic Equation (NPE) model for propagation in porous ground layers similar to Eq. (1.28).

  Fig. (1.13) shows a sketch and a diagram that detail the basic principle of each model.

  Random parameters H, W , D modeling a urban city Probabilistic model depending on a parameter u. u depends on the mean values and dispersions of random variables H,W and D.

  Figure 1.13: Sketches and diagrams detailing the basic principles of each propagation model. The reference propagation model is presented on the left, while the stochastic, simplified NPE model is shown on the right.

2

  2 describes the original Nonlinear Parabolic Equation (NPE) model derived in 1987 by McDonald et al. . Section 2.1 presents the original derivation of the model and in the second and third sections (Sections 2.2 and 2.3) the different extensions and re-formulations of the NPE model are summarized. Section 2.4 details the available literature concerning the use of the NPE model for sound propagation in multiple media. Section 2.5 deals with the procedure to couple the NPE model with other sound propagation models. Finally, significant applications using the NPE model are listed and briefly summarized in Section 2.6. Chapter summary and conclusions are given in Section 2.7. Chapter 3 details the development of a deterministic Nonlinear Parabolic Equation (NPE) model for high-amplitude wave propagation over porous ground layers and in complex media. Propagation over non-flat surfaces is handled through the use of the terrain-following coordinates method and described in Section 3.1. The NPE model for propagation within porous ground layers is detailed in Section 3.2 and propagation over porous layers in Section 3.3. Chapter summary and conclusions are given in Section 3.4. Chapter 4 concerns the numerical implementation of the NPE models developed in the previous chapters. The main principle of the numerical solution is the separation of the differential operators with the operator splitting method, described in Section 4.1. The numerical solution of nonlinear terms is detailed in Section 4.2 while the numerical solution of linear terms appears in Section 4.3. Initial conditions and boundary conditions are treated in Sections 4.4 and 4.5, respectively, and Section 4.6 gives some notes about software development. Chapter conclusions are given in Section 4.7. Chapter 5 presents the validation of the deterministic Nonlinear Parabolic Equation (NPE) model. The full validation of the computational model is separated in successive tasks and simulation results are compared to analytical and numerical solutions. In Section 5.1 nonlinear effects calculations are assessed using quasi-plane waves and analytical solutions. Simulations of propagation in an inhomogeneous medium are compared with a linear, frequency-domain implementation of the Parabolic Equation (PE) and presented in Section 5.2. In Section 5.3 propagation over a hilly ground is studied with the terrain-following coordinates version of the NPE which is compared to solutions given by a BEM implementation. Section 5.4 presents the validation of the NPE model for propagation over a finite-impedance ground surface in which both linear and nonlinear examples are given. Finally, chapter summary and conclusions are given in Section 5.5. Chapter 5 finalizes the development of the deterministic NPE model. Chapter 6 concerns the development of a computational model for long-range nonlinear sound propagation over urban environments. Section 6.1 briefly summarizes the previous chapters, Existing work on the Nonlinear Parabolic Equation model for sound propagation on plane and acoustically rigid surfaces This chapter describes the original Nonlinear Parabolic Equation (NPE) model derived in 1987 by McDonald et al. , its extensions and variants, together with main applications using this model. The original derivation is given in Section 2.1. In the second and third sections (Sections 2.2 and 2.3) the different extensions and re-formulations of the NPE model are reviewed. Section 2.4 details the available literature concerning the use of the NPE model for sound propagation in multiple media, Section 2.5 deals with the procedure to couple the NPE model with other sound propagation models, and finally, significant applications using the NPE model are listed and briefly summarized in Section 2.6. Chapter summary and conclusions are given in Section 2.7. Chapter content 2.1 Derivation and properties of the original Nonlinear Parabolic Equation model . . 44 2.2 Extensions to the original Nonlinear Parabolic Equation model . . . . . . . 50 2.3 Variants of the original Nonlinear Parabolic Equation model . . . . . . . 52 2.4 Nonlinear Parabolic Equation model for propagation in multiple media . . . . 53 2.5 Coupling the Nonlinear Parabolic Equation model with near-field and far-field sound propagation methods . . . . . . . . . . . . . . . . . . . . 55 2.6 Significant applications using the Nonlinear Parabolic Equation model . . . . 56 2.7 Chapter summary and conclusions . . . . . . . . . . . . . . . . 58 2.1 Derivation and properties of the original Nonlinear Parabolic Equation model Two different derivations that yield the original Nonlinear Parabolic Equation (NPE) model have been published. The first, published in 1987, makes use of the perturbation expansion method [see the papers by McDonald & Kuperman

  dx is the diffraction operator. It accounts for propagation in the transverse direction (the z-direction). These operators are graphically identified in Fig. (2.1). NPE models have the unique feature of segregating physical effects into separate differential operators. This makes the NPE a very modular and flexible model. One can, for example, obtain a one-dimensional nonlinear propagation equation by omitting the diffraction operatorc 0 2 ∂ 2 z R dx , or obtain a two-dimensional linear parabolic propagation equation by omitting the nonlinear term

  23) differs from the inviscid Burgers equation Eq. (2.22) only by the presence of the term c 1 , which accounts for sound speed variations, and a scaling factor in front of the nonlinear term, since the unknown variables used are different. Relation to the KZK equation Blanc-Benon et al. [2002] used a modified KZK equation which accounts for sound speed variations in order to study propagation of finite-amplitude sound waves through turbulence.

  .39) in which R w and R b are adimensional overpressures in water and bottom fluid layers, respectively, and ρ 0w and ρ 0 b are their respective ambient densities. Ambrosiano et al. estimated

2. 5

 5 Coupling the Nonlinear Parabolic Equation model with near-field and far-field sound propagation methods The NPE model has previously been coupled to various propagation methods. Attenborough et al. [2005, 2004] and Van der Eerden & Védy [2005] coupled an Euler equations implementation * to the NPE model. The NPE was then coupled to a (linear) frequency-domain Parabolic Equation (PE) implementation. Fig. (2.2) shows a schematic of the principle of coupling Euler equations, NPE and PE models.

Figure 2

 2 Figure 2.2: Schematic overview of the method for coupling Euler equations method with the NPE model and the PE.

3

  Development of a Nonlinear Parabolic Equation model for sound propagation in complex environments -Deterministic aspectsThis chapter details the development of a deterministic Nonlinear Parabolic Equation (NPE) model for high-amplitude wave propagation over porous ground layers and in complex media.Propagation over non-flat surfaces is handled through the use of the terrain-following coordinates method and is described in Section 3.1. The NPE model for propagation within porous ground layers is detailed in Section 3.2 and propagation over porous layers in Section 3.3.Chapter summary and conclusions are given in Section 3.4.Chapter content3.1 Nonlinear Parabolic Equation model for propagation over rigid non-flat ground surfaces 60 3.2 Nonlinear Parabolic Equation model for sound propagation in rigidly-framed porous media . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.3 Nonlinear Parabolic Equation model for high-amplitude wave propagation over complex surfaces . . . . . . . . . . . . . . . . . . . . . . . 70 3.4 Chapter summary and conclusions . . . . . . . . . . . . . . . . 79 3.1 Nonlinear Parabolic Equation model for propagation over rigid non-flat ground surfaces In this section a Nonlinear Parabolic Equation (NPE) model for high-amplitude wave propagation over non-flat terrains is developed.

  .1b) in which H is the height of the top of the model domain. Using the transformation in Eqs. (3.1) introduces model levels that change from coordinate-following near the ground to plane at the top of the domain [see for example Pielke, 2002]. For the simulation of acoustic wave propagation a rather shallow atmospheric layer is involved. It is therefore more convenient to use a coordinate transformation that conserves the spatial resolution over the calculation grid. The spatial coordinates are hence changed according to x -→ x , (3.2a) z -→ z + h(x) , (3.2b) in such a way that model levels are following the ground elevation from the bottom to the top of the model domain (see Fig. (3.1) for an illustration).

  Figure 3.2: Graphical representation of additional terms accounting for ground topography in the GT-NPE model.

  .24c) Substituting Eqs. (3.23) and Eqs. (3.24) in Eq. (3.18a) yields

  (3.3) for a graphical representation of the NPE operators):

  Figure 3.3: Graphical representation of the differential operators of the NPE model for propagation in porous ground layers.

  1 and in publications by McDonald & Kuperman and McDonald et al. are used.

Figure 3 . 4 :

 34 Figure 3.4: Illustration of the discretization of the boundary between the atmospheric layer and the porous ground layer. The fluid-fluid interface is taken to be midway between two vertical grid points with indexes j = 0 and j = 1. Auxiliary virtual points (red circles) p ′ a i,0 and p ′ g i,1 are created.

  Fig. (3.5) shows the causality index C nl for a porous material defined by the ZK model with toruosity Φ= 1.2, porosity Ω= 0.3 and flow resistivity σ 0 =300 kPa.s.m -2 , between 0 and 10 000 Hz. As it can be seen the causality index C nl is equal to unity on the whole frequency range (departure from causality near 0 Hz is due to the numerical evaluation of the Hilbert transform), indicating that the Zwikker-Kosten (ZK) model is causal.

Figure 3 . 5 :

 35 Figure 3.5: Causality index (nonlocal form) of the Zwikker-Kosten (ZK) model.

  their derivatives) gives expressions for unknown quantities R g i,0 and R b i,1 and the corresponding coupling coefficients (see Eqs. (3.70, 3.71) on page 81).

4

  Discretization of the NonlinearParabolic Equations with the finite-difference method This chapter concerns the numerical implementation of the NPE models derived in the previous chapters. The main principle of the numerical solution is the separation of the differential operators with the operator splitting method, described in Section 4.1. The numerical solution of nonlinear terms is detailed in Section 4.2 while the numerical solution of linear terms appears in Section 4.3. Initial conditions and boundary conditions are treated in Sections 4.4 and 4.5, respectively, and Section 4.6 gives some notes about software development.

  (4.1) shows the analytical and numerical solutions of the problem defined by

Figure 4 . 1 :

 41 Figure 4.1: Analytical (solid line) and numerical (circles) solution of the Burger equation obtained with the MacCormack (top) and leapfrog (bottom) schemes. This figure is taken from Sjögreen [1990].

  When simulating wave propagation over non-flat terrains one has to use the terrain-following coordinates version of the NPE (see Eq. (3.11) in Section 3.1.2). The equation shows a new differential operator and a modified diffraction term. These changes are shown below in Eq. (4.26), without refraction and nonlinear effects:

  Incorporating the boundary interface condition in the diffraction operator Nonlinear sound propagation over porous ground surfaces involves the boundary interface condition developed in Section 3.3. It is included by forcing values on corresponding grid points with the help of Eqs. (3.53, 3.54) (or eventually with Eqs. (3.64, 3.65) for non-flat surfaces and Eqs. (3.70, 3.71) for multi-layered ground surfaces) which from a numerical point of view can be tricky. To summarize the procedure, solving for diffraction with the boundary interface condition consists in four steps: Step 1: update of the sums appearing in Eqs. (3.53) (or Eqs. (3.64) or Eqs. (3.70)),

  3.6, the value of σ 0 appearing in Eqs. (3.54) (or, again, Eqs. (3.65) or Eqs. (3.71)) must be updated at each time iteration with the help of the field values on the interface at the previous time step.

  .34) so that the signal can only move backward in the window and does not reach the right boundary. One hence has c win = max (c (x, z)) , (4.35)

  figure shows snapshots of the propagation of a blast wave in a PML (the complete wave is shown in the top figure; zoomed version in the bottom figure). The wave is shown at times

Figure 4 . 2 :

 42 Figure 4.2: Snapshots of the propagation of a blast wave in a PML (complete waves shown in the top figure; zoomed versions in the bottom figure). The wave is shown at times 0.11, 0.17 and 0.23 s (from left to right). The PML is 5 m thick (50 points, dashed white line). The contours represent the wave amplitude in decibels, relative to the maximum amplitude of the wave. 101

  complete simulation software, called NOnLInear Time domain Acoustics (NOLITA) implements solvers of each NPE model described in this document and comes with Python scripts for pre and post-processing. The whole program flow is far too complex to be reproduced here. Indeed the number of NPE versions (non-flat terrains, multi-layered surfaces, refraction effects, thermoviscous effects, PMLs, etc.) yields multiple and complex branching in the program structure. However, one can see in Fig. (4.6.2) a simplified flow chart for solving for the NPE.

5

  Two initialization methods have been reviewed: the Kinney-Graham (KG) model and the numerical solution of Euler equations. The KG model gives time waveforms for propagation of high-amplitude waves in free field. Hence, it cannot be used to determine waveforms in the far-field, where meteorological effects and ground effects are of importance. Instead, the numerical solution of Euler equations has to be employed. It provides accurate waveforms but suffers from its numerical complexity.Boundary conditions on the computational grid are then detailed. On lateral and bottom boundaries classical conditions are used. A Perfectly Matched Layer (PML) for NPE models is developed for domain truncation. The PML allows us to truncate the domain without introducing spurious oscillations while keeping the additional numerical effort low.The validation of the implemented numerical solutions of the various NPE models is presented in the next chapter.Validation of the deterministic Nonlinear Parabolic Equation modelThis chapter presents the validation of the deterministic nonlinear parabolic equation (NPE) model. The full validation of the computational model is separated in successive tasks. In Section 5.1 nonlinear effects calculations are assessed using quasi-plane waves and analytical solutions. Simulations of propagation in an inhomogeneous medium are compared with a linear, frequency-domain implementation of the parabolic equation (PE) and presented in Section 5.2. Propagation over a hilly ground is studied with the terrain-following coordinates version of the NPE, the Generalized Terrain -Parabolic Equation (GT-NPE), and compared to solutions of a BEM implementation in Section 5.3. Section 5.4 presents the validation of the NPE model for propagation over a finite-impedance ground surface in which both linear and nonlinear examples are given. Chapter summary and conclusions are given in Section 5.5. Chapter content 5.1 Propagation of high-amplitude quasi-plane waves in a homogeneous medium . . 106 5.2 Sound propagation in an inhomogeneous medium over a flat and acoustically rigid ground surface . . . . . . . . . . . . . . . . . . . . . . 107 5.3 Sound propagation over a non-flat acoustically rigid ground surface in a homogeneous medium . . . . . . . . . . . . . . . . . . . . . . . . 112 5.4 Sound propagation over a porous ground surface in a homogeneous medium . . . 117 5.5 Chapter summary and conclusions . . . . . . . . . . . . . . . . 128

  Fig. (5.1), which presents three snapshots of the wave at time 10, 16.4 and 22.8 seconds. As it can be seen, the NPE and the Fubini solutions are very close: the wave steepens with the expected rate. The criteria to evaluate the performances of the NPE is the harmonics amplitudes evolution with time. During wave steepening, the amplitude of the fundamental frequency decreases and the amplitude of higher harmonics increases (see Section 1.2.1). Fig. (5.2) presents the evolution of the fundamental frequency together with the first and second harmonics amplitudes over time. The first 120 time steps are shown in Fig. (5.2): this corresponds to the shock formation distance, beyond which the Fubini solution is no longer valid.

Figure 5 . 1 :

 51 Figure 5.1: Steepening of a sine wave calculated with the Fubini solution (red dashed line) and the NPE model (blue line). The three plots show the wave at times 10, 16.4 and 22.8 seconds (from top to bottom).

Figure 5 . 2 :

 52 Figure 5.2: Fundamental and harmonics amplitudes evolution during wave steepening. The curves present solutions given by the Fubini formula Eq. (5.2) (solid lines) and NPE calculations (dashed lines). Fundamental amplitude: blue line; first harmonic amplitude: red line; second harmonic amplitude: green line.

  Fig. (5.3) shows the TL at an altitude of 2 m, calculated with ATMOS, the NPE coupled with ATMOS at 46 km, and the NPE coupled with ATMOS at 112 km. As it can be seen all calculations lead us to the same TL except very close to the coupling distance. The PE and NPE methods being parabolic approximations field values near the source, or in this case near the coupling location, cannot be accurate.

Figure 5 . 3 :

 53 Figure 5.3: Transmission loss (TL) calculated with a single NPE simulation and two coupled NPE/ATMOS simulations. No coupling: blue line; coupling at 46 km: red line; coupling at 112 km: green line.

  propagation of a low frequency wave in an inhomogeneous medium is considered. The source is a Gaussian spatial distribution whose central frequency is 10 Hz. It is set at altitude z = 400 m close to the right edge of the computational domain (the signal propagates backward in the window). The sound speed profile chosen is shown in Fig.(5.4) (the sound speed only varies with altitude). Although this profile may not be extremely realistic it has the advantage of presenting a lot of variations and thus allows the accuracy of refraction effects calculations to be fully estimated. The window is moving with the maximum sound speed in the domain which is 363.2 m.s -1 and the maximum sound speed perturbation is max c 1 (x, z) = 22.2 m.s -1 and occurs at an altitude of 1550 m. The medium properties at rest are those found for a temperature of 20 • C under 1 atmosphere (the nonlinearity coefficient β is here set to zero, since we only consider linear propagation). The NPE calculation window is 2001 points large in both directions and the spatial discretization steps have both been chosen equal to 1 m. Considering the wave central frequency (10 Hz), this yields a spatial resolution of 34 points/λ and a time step of approximately 2.9 ms. The simulation ran for 30 000 iterations and resulted in the calculation over a domain which is 2 km high and 25 km wide; the wave has propagated for approximately 75 seconds. In order start the ATMOS simulation which is used as a reference, the time signals are recorded along a vertical line at a distance of 1 km from the source and then Fourier transformed and used as initial condition for ATMOS.

Figure 5 . 4 :

 54 Figure 5.4: Sound speed profile used for ATMOS and NPE calculations.

  (a) Transmission Loss (TL) value calculated with the NPE model. (b) Transmission Loss (TL) difference between NPE and ATMOS calculations.

Figure 5 . 5 :

 55 Figure 5.5: Transmission Loss (TL) calculated with the NPE model and ATMOS. The top figure shows the TL calculated with the NPE model. Eight contour lines regularly distributed between -60 dB and +10 dB are shown. The bottom figure shows the difference between the values calculated with the NPE model and ATMOS. Five contours are shown; blue areas denote a difference lower than 2 dB.

Figure 5 . 6 :

 56 Figure 5.6: Transmission loss (TL) at three different altitudes (from top to bottom: 1.5 km, 1 km and 500 m) for ATMOS (blue line) and NPE (dashed red line) simulations.

  Fig. (5.7): the figure shows contour plots of the field of the overpressure R taken at times 180 ms, 265 ms, 310 ms and 355 ms. The thick black line represents the rigid ground. Sound Pressure Level (SPL) maps have been calculated for both the GT-NPE model and MICADO and are shown in Fig. (5.8). As one can see the results from the NPE model are

Figure 5 . 7 :

 57 Figure 5.7: Contour plots of the field of the overpressure R taken at times 180 ms, 265 ms, 310 ms and 355 ms (time evolution is from top to bottom). The thick black line represent the (rigid) ground.

  (a) Sound Pressure Level (SPL) map calculated with the Boundary Element Method (BEM). (b) Sound Pressure Level (SPL) map calculated with the Generalized Terrain -Nonlinear Parabolic Equation (GT-NPE) model.

Figure 5 . 8 :

 58 Figure 5.8: Sound Pressure Level (SPL) maps for the BEM (top figure) and the GT-NPE model (bottom figure) for sound propagation over an acoustically rigid hilly ground surface. Eight contours are shown regularly distributed between 40 dB and 95 dB.

Figure 5 . 9 :

 59 Figure 5.9: Sound Pressure Level (SPL) at three different altitudes (from top to bottom: 1.5 m, 15 m and 30 m) for BEM (dashed blue line) and NPE (solid red line) simulations.

  et al. [2002] found very little difference on relative SPLs when comparing sources with different geometrical decay rates. Analytical solutions for cylindrical line sources are thus valid references for comparison with the NPE model used in this work. Relative SPLs at the receiver are shown in Fig. (5.10), for both analytical and NPE calculations. Very good agreement can be observed independently of the ground properties. Even for the softest layer (σ 0 = 100 kPa.s.m -2 ) the difference between analytical and NPE calculations is at most 1 dB. The frequencies where negative interference occurs are 1325 Hz, 1273 Hz and 1246 Hz for the rigid case, the ground layer with σ 0 = 500 kPa.s.m -2 and the ground layer with σ 0 = 100 kPa.s.m -2 , respectively. As one can see in Fig. (5.10) the NPE model presented does not only accurately recreate reflected wave amplitude decrease, but does account for the change of least reflective frequencies due to the additional delay given during reflection.

  Fig. (5.11) shows the one-dimensional signal and its twodimensional extension used to start both reference and NPE calculations.

Figure 5 . 11 :

 511 Figure 5.11: Pressure waveform used to start the reference and NPE simulations. Left: onedimensional waveform; right: two-dimensional extrapolation.

Figure 5 . 12 :

 512 Figure 5.12: Pressure field at time t = 33 ms. Left: σ 0 = 100 kPa.s.m -2 ; Right: σ 0 = 10 kPa.s.m -2 ; Color map: solution from NPE model; contour lines: solution from Euler equations. Ten contour lines equally spaced from -800 to 800 Pa are shown. Contours corresponding to negative values are represented by dashed lines, positive ones by solid lines.

Figure 5 . 13 :

 513 Figure 5.13: Time signals at the receiver for NPE and reference calculations for the three ground layers considered (from top to bottom: perfectly rigid, σ 0 = 100 kPa.s.m -2 and σ 0 = 10 kPa.s.m -2 ). Solid line: Euler; Dotted line: NPE.

Fig. ( 5

 5 Fig. (5.14) shows time signals at the receiver and Tab. (5.3) summarizes their characteristics for both reference and NPE calculations. One can see that compared to the same ground layer with no Forchheimer's nonlinearities (bottom plot in Fig. (5.13)), the obtained reflected wave has a larger positive peak amplitude and a shorter time of arrival. The relative error

  Fig. (6.1) shows a sketch and a diagram that detail the basic principle of each model.

B

  i and B i+1 is noted d i (see Fig. (6.2)). The prior probability model of such a urban city is then introduced with the help of the random variables {H 1 , . . . , H n }, {W 1 , . . . , W n } and {D 1 , . . . , D n-1 }.

  Random parameters H, W , D modeling a urban city Probabilistic model depending on a parameter u. u depends on the mean values and dispersions of random variables H,W and D.

Figure 6 . 1 :

 61 Figure 6.1: Sketches and diagrams detailing the basic principles of each propagation model. The reference, linear propagation model is presented on the left, while the stochastic, simplified NPE model is shown on the right.

  .7) Parameter u belongs to an admissible set U = ( ] 0, +∞ [ ) 6 . Fig. (6.3) shows an example of a realization with n = 20 and u= (10, 20, 30, 0.2, 0.2, 0.2).

Figure 6 . 3 :

 63 Figure 6.3: Example of a city realization with n = 20 and u= (m H , m W , m D , δ H , δ W , δ D ) = (10, 20, 30, 0.2, 0.2, 0.2).

Figure 6 . 4 :

 64 Figure 6.4: Sketch of the different domains with their respective boundaries of the NPE model for propagation in multiple media.

  θ is real and positive and has no physical meaning. Its only intent is to help controlling the mean model output. The mean propagation model parameter s can now be defined as s = (γ , λ , θ) . (6.15) Vector s belongs to the admissible set S = ([0, 1] × ] 0, +∞ [ × [ 0, +∞ [ ) and the initial value of s is arbitrarily chosen as

  .26) and the nominal value of the mean model parameter s 0 can be updated in a vector s opt such thats opt = arg min s∈S ||E {L exp } -L s || B .(6.27)Hence for a fixed value of u= (m H , m W , m D , δ H , δ W , δ D ), the parameter s opt that minimizes the norm between the mathematical expectation of the experimental observations and the observation from the mean model can be determined.

  .38) in which T superscript stands for transposition. The following eigenvalue problem is introduced C L(w) x (w) = e (w) x (w) ,(6.39)

6. 9 . 4

 94 Fig. (6.9) shows the experimental observations {L exp 1 (ω ; η j )} j and {L exp 2 (ω ; η j )} j , their mean values estimates and the confidence region calculated using the stochastic model and the mean value of the stochastic model. As one can see the experimental observations belong to the constructed confidence regions (except for the very lowest frequencies outside the range of interest). For the configuration corresponding to u 2 , the confidence region calculated using the stochastic model shows a broad confidence region in the very low frequency range. This means that the underlying deterministic model used is not robust in this very low frequency range with respect to statistical fluctuations generated by the probabilistic model.

Figure 6 . 6 :

 66 Figure 6.6: Comparisons of the reference model with the nominal model (initial mean model). Reference model (made up of 500 outputs (thin lines) and the mean values estimates m L exp 1 and m L exp 2 (thick blue lines)). Nominal model output L s 0 (dashed lines). The value of parameter u is u 1 (top figure) and u 2 (bottom figure).

Figure 6 . 7 :

 67 Figure 6.7: Comparisons of the reference model with the updated mean model. Reference model (made up of 500 outputs (thin lines) and the mean values estimates m L exp 1 and m L exp 2 (thick blue lines)). Updated mean model output L s opt (dashed lines). The value of parameter u is u 1 (top figure) and u 2 (bottom figure).

Figure 6 . 8 :

 68 Figure 6.8: Values of the objective functions -L red (w) and J (w) for the 50 individuals at generations 1 (circles), 10 (squares) and 50 (final generation, stars). The value of parameter u is u 1 (left figure) and u 2 (right figure). Note that some individuals are out of the figure ranges.

Figure 6 . 9 :

 69 Figure 6.9: Comparisons of the reference model with the confidence region calculated with the stochastic propagation model. Reference model (made up of 500 outputs (thin lines) and the mean values estimates m L exp 1 and m L exp 2 (thick blue lines)). The confidence region calculated with the stochastic model is represented by the yellow area, delimited by thick black lines. The value of parameter u is u 1 (top figure) and u 2 (bottom figure). 163

  been tested regarding meteorological conditions and compared to a frequency-domain implementation of the Parabolic Equation (PE). It appeared that Sound Pressure Levels (SPL) are accurate within a 2 dB relative error range except on small localized areas. The Boundary Element Method (BEM) is then used to perform simulation on a hilly terrain. It has been shown that the terrain-following coordinates formulation of the NPE model is in good agreement with the BEM results in the far-field. The performances of the NPE model for propagation over porous ground layers are then assessed. For linear propagation, the results obtained with this method have shown very good agreement with analytical solutions for a wide range of ground properties. For high-amplitude waves, the NPE model produces time signals comparable to those obtained by the numerical solution of Euler equations. The approximate method to include the Forchheimer nonlinearities in the two-way coupling has then been proved to give satisfactory results and does not introduce any additional source of error in the two-way coupling. The model validation finalizes the development of the deterministic sound propagation model for high-amplitude wave propagation in complex media. The model can account for most of the features of sound propagation outdoors (nonlinearities, refraction, dissipation, topography and ground impedance effects) and has been shown to provide results that are in very good agreement with other numerical models. Furthermore, the speed of the NPE model makes it a good candidate for a basis of a stochastic sound propagation model. The deterministic model is then associated a probabilistic model of uncertainties to handle sound propagation over urban environments. First, the urban city geometrical parameters distributions are determined using Information Theory and the Maximum Entropy Principle.

  ρ T w) = 0 . (D.1) Notations and variables are identical to the ones used in Chapter 2. It is assumed that the flow is irrotational in the z-direction; this allows us to get an expression for the z-component of the flow velocity vector, such that∂ z u = ∂ x w , in Eq. (D.1) yields ∂ t ρ T + ∂ x (ρ T u) + ∂ z ρ T ∂ z u dx = 0 . (D.4)One has now to determine an expression for the quantity u (ρ T ). The momentum equation along the x-axis is used. This equation is defined by∂ t (ρ T u) + ∂ x p T + ρ T u 2 + ∂ z (ρ T uw) = 0 . (D.5)Integrating with respect to time leads us to an expression for ρ T u, such thatρ T u = -∂ x p T + ρ T u 2 + ∂ z (ρ T uw) dt . (D.6)Substituting Eq. (D.6) in Eq. (D.4) yields∂ t ρ T -∂ x ∂ x p T + ρ T u 2 dt -∂ x ∂ z ρ T u ∂ z u dx dt + ∂ z ρ T ∂ z u dx = 0 . (D.7)The last term in equation Eq. (D.7) is associated with diffraction in the z-direction. Since the term ∂ x ∂ z ρ T u ∂ z u dx dt represents contributions not predominantly in the propagation direction and is of second order, one chooses to discard it. One obtains∂ t ρ T -∂ x ∂ x p T + ρ T u 2 dt + ∂ z ρ T ∂ z u dx = 0 . (D.8) The pressure p T and the flow velocity component u are now developed with Taylor series, such that p T = p (ρ 0 ) + ∂p ∂ρ 0 (ρ T -ρ 0 in Eq. (D.8) yields ∂ t ρ T -∂ x ∂ x p T + ρ T u 2 dt + c 0 ∂ 2 z ρ ′ dx = 0 . (D.12)

  14)which yields∂ t ρ ′ dx = -c 0 ρ ′ .(D.15)

  20)Assuming a spatially-varying sound speed c (x, z), we define a sound speed perturbationc 1 (x, z), such that c 1 (x, z) = c (x, z) -c 0 . (D.21)So that a second order accuracy can be retained in Eq. (D.20), p T is transformed with the help of Eqs. (D.9, D.10) with the derivatives coefficients corresponding to an isentropic flow,

u

  N := y N /m N For i = N -1, N -2, • • • , 1 do: u i := (y i -β i u i+1 ) /m i (E.7)

  . The pressure p (M ) at any point M in the domain Ω must satisfy the Helmholtz equation, defined by ∆ + k 2 p (M ) = f (M ) , M ∈ Ω , (F.1) in which f (M ) is a source distribution and k is the wave number. Using the Green function G and the Sommerfeld radiation condition and after a few algebraic manipulations one obtains c (M ) p (M ) = p 0 (M ) + σ [p (M ) ∂ n S G (S, M ) -G (S, M ) ∂ ns p (M )] dS , M ∈ Ω , (F.2) where n s is the outward normal to the surface σ and c (M ) is a coefficient depending on the receiver position M such that c (M ) = 1 for M in Ω, except on its boundary , c (M ) = 1 -θ 4π if M is an angular point . (F.3c) Résumé long en français 1 Introduction générale et travail proposé Ce travail intitulé Propagation d'ondes acoustiques non linéaires en milieu complexe -Application à la propagation en environnement urbain trouve son origine dans des évènements comme l'explosion de l'usine AZF à Toulouse en septembre 2001. Cette explosion fut responsable de trente victimes, de milliers de blessés et coûta plus de 2 milliards d'euros à la communauté. Il y a encore aujourd'hui certains faits de cet évènement qui restent incompris. La présence d'un double écho entendu par certains témoins reste notamment sans explication. Ce travail de recherche tente d'apporter une réponse à cette problématique en proposant le développement d'un outil de simulation numérique pour l'étude de la propagation d'ondes de choc en milieu complexe. Cette application rentre dans le cadre général de l'acoustique et plus particulièrement de l'acoustique environnementale. Propagation acoustique en milieu complexe La propagation d'ondes acoustiques en milieu extérieur tel que sur un environnement urbain et sur de longues distances implique la prise en compte de nombreux phénomènes du milieu de propagation : -l'absorption atmosphérique, pouvant être importante compte tenu des distances considérées, -les effets de refraction dus aux conditions météorologiquies (gradients de vitesses de vent ou de température) -les effets du sol (par sa topographie et son impédances) modifiant la façon dont les ondes sont réfléchies. Compte tenu des grandes amplitudes en jeu lors de l'explosion, les effets non linéaires sont importants et provoquent la formation d'un choc, ayant pour conséquences : -la génération d'harmoniques, -la coalescence des chocs, -et une dissipation anormale au front de choc. Les effets combinés des propriétés du milieu de propagation et des effets des grandes amplitudes impliquent l'utilisation de modèles et méthodes numériques avancées. xxvii Résumé long Travail de recherche proposé Un phénomène de la propagation du son en milieu extérieur non traité jusqu'à maintenant concerne la diffusion liée à de multiples diffractions en milieu ouvert causée par l'environnement urbain sur le trajet de l'onde. Plus précisément, on s'intéresse à l'effet de la canopée urbaine sur le champ acoustique dans la couche atmosphérique planétaire. Cette problèmatique pourrait être étudiée à l'aide de modèles numériques dans lesquels la géométrie du site est rentrée explicitement dans le modèle. Par exemple, les méthodes de type lancer de rayons ou la solution des équations d'Euler sont techniquement viables pour ce genre d'applications. Ces méthodes souffrent cependant de temps de calculs prohibitifs, et compte tenu des incertitudes sur les paramètres du modèle (dimensions des bâtiments) le modèle pourrait être amélioré par l'utilisation d'un modèle probabiliste d'incertitudes.On propose donc dans ce travail une approche différente de la solution purement numérique.On choisit de modéliser la propagation dans le milieu composé de l'atmosphère et de la couche urbaine dont le modèle moyen est simplifié. La propagation dans la couche atmosphérique est prise en compte par un modèle d'Équation Parabolique Nonlinéaire (NPE). La propagation dans la couche urbaine est prise en compte par un modèle NPE de propagation dans un milieu poreux.Ce modèle permet de construire une première approximation qui sera améliorée en prenant en compte la complexité du système réel. C'est pour cela que l'on introduit un modèle probabiliste d'incertitudes au modèleNPE. Le travail de recherche présenté dans ce document peut être séparé en deux tâches distinctes :1. Développement d'un modèle déterministe NPE pour la propagation d'ondes de grandes amplitudes en milieu complexe au-dessus d'un milieu poreux.2. Développement d'un modèle probabiliste associé au modèle NPE pour la prise en compte des effets de l'environnement urbain sur la propagation de l'onde. Le développement du modèle déterministe NPE est résumé dans le paragraphe 2 et la construction du modèle probabiliste d'incertitudes est détaillée dans le paragraphe 3. Une synthèse de ce travail et quelques perspectives de recherche sont présentées en paragraphe 4. 2 Développement d'un modèle NPE pour la propagation en milieu complexe -Aspects déterministes Méthodologie L'effet de sols non plans sur la propagation est pris en compte grâce à la méthode de transformées conformes. Cette transformation des coordonnées permet de prendre en compte la xxviii topographie du site de façon analytique (et non numérique) dans le modèle NPE. Il faut cependant noter que, pour que la méthode soit exacte et à cause de l'approximation paraxiale inhérente aux modèles NPE, la pente du terrain doit rester faible. On développe ensuite un modèle NPE pour la propagation dans les matériaux poreux. Ce modèle est basé sur une extension non linéaire du modèle de Zwikker-Kosten (ZK) [voir Zwikker & Kosten, 1949], dans lequel le sol est caractérisé par quatre paramètres : -la tortuosité Φ, -la porosité Ω 0 , -la résistivité à l'air σ 0 , -et le coefficient de Forchheimer ξ. Le modèle développé est similaire au modèle NPE pour la propagation atmosphérique, mais contient des termes supplémentaires pour prendre en compte les propriétés du sol poreux. Afin de développer un modèle NPE de propagation au-dessus de couches poreuses, on développe ensuite une condition d'interface permettant de coupler les deux domaines (air et sol poreux). Ces équations de couplages impliquent uniquement des dérivés et intégrales spatiales, rendant son implémentation au sein du modèle existant directe et naturelle. La formulation de cette condition d'interface pour les sols non plans et multi-couches est également donnée. Implémentation numérique du modèle La méthode des différences finies et le principe de séparation des opérateurs différentiels sont utilisés pour la résolution numérique du modèle NPE. La séparation des opérateurs consiste à utiliser des algorithmes spécialisés pour chaque terme. Ceci implique, de façon implicite, que les effets des différents opérateurs sont découplés sur le pas d'intégration temporel. Pour la solution des termes non linéaires, on utilise l'algorithme Flux-Corrected Transport (FCT). Cette méthode permet de propager des discontinuités de façon stable et sans introduire d'oscillations numériques. Pour les termes linéaires, on utilise un schéma de différences finies semi implicite : la méthode de Crank-Nicolson. Le schéma est inconditionnellement stable et donne des systèmes d'équations tri diagonaux, résolus par l'algorithme de Thomas.

  harmoniques supérieures est toujours inférieure à 5 %. Le modèle NPE implémenté permet donc de prendre en compte les effets non linéaires de façon précise.Le modèle est ensuite testé par rapport aux conditions météorologiques. On utilise une implémentation de l'équation parabolique fréquentielle comme référence. Il apparaît que les niveaux de pression calculés sont précis à 2 dB près, à l'exception de quelques zones locales où le niveau diffère par au plus 4 dB. Ce cas de validation prouve aussi la faisabilité du couplage entre le modèle NPE et d'autres codes de propagation, ici l'équation parabolique fréquentielle.Pour évaluer la précision du modèle NPE pour la propagation sur sols non plans, une implémentation de la méthode des éléments de frontière (BEM) est utilisée comme référence. Les résultats du modèle GT-NPE sont en accord avec les prévisions de la méthode BEM.La dernière étape de la procédure de validation du modèle concerne la propagation au-dessus de sols poreux. Pour la propagation d'ondes de faibles amplitudes, on utilise la solution des équations d'Helmholtz en deux dimensions comme référence. Les niveaux de pression sont en parfait accord avec les prévisions analytiques. Pour la propagation d'ondes de grandes amplitudes, une implémentation des équations d'Euler est utilisée. Les signaux issus des deux méthodes sont comparables. L'erreur relative sur le pic de pression positive est indépendante des propriétés du sol, alors que l'erreur relative sur la durée de phase positive croît lorsque la résistivité du matériau diminue. Cependant, même pour des valeurs très faibles de résistivité, le modèle NPE produit des résultats en accord avec la solution de référence. Il est également prouvé que les effets des non linéarités de Forchheimer sur les signaux sont correctement reproduits par le modèle de propagation développé. Ce paragraphe finalise la construction du modèle de propagation déterministe d'ondes de grandes amplitudes en milieu complexe. Le modèle permet de prendre en compte la majorité des phénomènes intervenant lors de la propagation du son en milieu extérieur de façon précise, et, de part sa simplicité, permet de réaliser des simulations en des temps très courts (de l'ordre de quelques minutes sur un seul processeur). Cette propriété permet d'utiliser le modèle comme base pour la construction d'un modèle de propagation stochastique résolu par la méthode de Monte-Carlo. Dans le paragraphe suivant, on détaille la construction du modèle probabiliste d'incertitudes pour la prise en compte de l'environnement urbain sur le trajet de l'onde. xxx 3 Développement d'un modèle probabiliste d'incertitudes -Application à la propagation stochastique non linéaire en milieu urbain Le modèle moyen est donc basé sur l'Équation Parabolique Nonlinéaire (NPE) et son extension à la propagation au-dessus de sols poreux. Ce modèle moyen permet de réaliser des simulations dans des temps très courts, mais en contrepartie ne donne qu'une première approximation pour la simulation de la propagation au-dessus d'environnements urbains. On cherche donc à améliorer le modèle de propagation moyen à l'aide d'un modèle probabiliste d'incertitudes.Cette section en décrit les principales étapes de construction.On utilise en premier lieu la Théorie de l'Information et le principe du maximum d'entropie afin de déterminer les densités de probabilité des paramètres géométriques de l'environnement urbain. Ces paramètres géométriques sont caractérisés par un vecteur de paramètres que l'on note u. Le principe du maximum d'entropie, qui maximise l'incertitude du système sous les contraintes définies par l'information disponible, conduit à l'indépendance des variables aléatoires décrivant l'environnement urbain. Différentes réalisations de canopées urbaines peuvent donc être générées, et ensuite utilisées dans un modèle de propagation exact (ici une implémentation de la BEM) afin d'obtenir des solutions de référence.Le modèle de propagation NPE stochastique est paramétré par trois variables aléatoires. Leurs propriétés algébriques permettent l'utilisation du maximum d'entropie pour la détermination de leurs distributions de probabilités respectives. Ces distributions dépendent d'un vecteur de paramètres noté w.Le paramètre w est ensuite identifié en inverse par rapport au paramètre u. En d'autres termes, on cherche à connaître (de façon numérique) la relation entre les paramètres de l'environnement urbain (par exemple hauteur moyenne des bâtiments et dispersion associée) et les paramètres de la couche poreuse du modèle stochastique. On utilise deux méthodes pour résoudre ce problème inverse : -la méthode dont la fonction coût est la norme de la moyenne d'ordre deux, qui permet de minimiser les zones où les observations expérimentales n'appartiennent pas à la région de confiance du modèle stochastique, -et la méthode du maximum de log-vraisemblance, qui permet de maximiser la vraisemblance entre les observations expérimentales et les observations du modèle stochastique. Concernant la méthode du maximum de log-vraisemblance, on utilise une méthode de réduction de l'information, afin de réduire l'effort de calcul. La méthode consiste à procéder à une analyse en composantes principales des observations. Cela permet d'utiliser la méthode du maximum de log-vraisemblance dans un espace réduit de variables aléatoires décorréllées, et ainsi d'éviter l'évaluation fastidieuse de leur densité de probabilité jointe. Ces deux xxxi Résumé long méthodes sont ensuite utilisées conjointement dans un algorithme d'optimisation génétique multi-objectifs. Une fois que le paramètre optimal w opt est identifié il est utilisé pour résoudre les équations stochastiques par la méthode de Monte-Carlo. Les régions de confiance sont ensuite construites et comparées aux données de référence. La région de confiance enveloppe les observations expérimentales sur la quasi totalité du spectre, à l'exception d'une zone en très basses fréquences, et ce quel que soit la dispersion du système réel. 4 Synthèse du travail et perspectives Synthèse du travail L'objectif de ce travail de recherche était de développer un modèle de simulation numérique pour la propagation des ondes de grandes amplitudes en milieu complexe incluant un environnement urbain. La solution proposée dans ce document utilise un modèle de propagation simplifié, l'Équation Parabolique Nonlinéaire (NPE). Ce modèle est d'abord étendu pour pouvoir prendre en compte les effets de sol et d'environnements urbains dans la propagation (topographie et diffusion par le milieu urbain). Les solutions numériques de ce modèle sont obtenues grâce à la méthode des différences finies. On a montré que ces solutions sont en accord avec les résultats donnés par d'autres méthodes numériques. On améliore ensuite les performances du modèle pour la propagation au-dessus d'environnements urbains en lui associant un modèle probabiliste d'incertitudes. Ce modèle est construit en utilisant les méthodes générales à la théorie des probabilités. On montre enfin que ce modèle de propagation stochastique NPE présente de très bons résultats, quelque soit le niveau de dispersion sur la diffusion induite par les multiples diffractions dans le milieu urbain.

  Je tiens ensuite à remercier Jérôme Defrance, mon encadrant au sein du CSTB, ainsi que Philippe Jean et Christian Soize, respectivement co-directeur et directeur de cette thèse.

	H i ω		Random variable -height of building B i Pulsation
	Heav (x)	Heaviside function
	I/A		Impulse per unit area
	J (w) J n Greek uppercase letters Objective function for the mean square method Bessel function of order n
	L (ω, w) L exp (ω) Γ (z) L Γ, Λ, Θ L red ∆t -L red ∆x N x ∆z P r (ω) Φ P exp r (ω) Φ ij Q Ω 0	Relative SPL at the receiver. Output from the stochastic NPE model Gamma function of argument z Relative SPL at the receiver. Output from the reference model Random parameters of the ground layer Log-likelihood function Time step Reduced log-likelihood function Spatial step in the x-direction Objective function for the maximum likelihood method Spatial step in the z-direction Number of points in the moving window in the x-direction Material tortuosity Pressure at the pressure -stochastic model observation Rate of shear tensor Pressure at the pressure -experimental observation Cylindrical reflection coefficient Ground porosity
	R			Dimensionless density perturbation
	R p V Mathematical signs Plane wave reflection coefficient Flow velocity vector
	h ′′ k m X i p ′ Z c Z 0 W i W ∂ i Z D t Z δ	Ground elevation second derivative with respect to x Pressure perturbation Material characteristic impedance Mean value of random variable X i Characteristic impedance of air Wavenumber Random variable -width of building B i Partial derivative with respect to variable i Equivalent charge of TNT Frame-following operator Reduced distance in m/kg 1/3 Surface acoustic impedance Kronecker delta
		p free (ω)	Free field pressure at the receiver
		p 0	Peak overpressure
	p 0 p T Greek lowerwase letters Ambient pressure Total pressure
	α β γ	p X i s t +	Probability density function of random variable X i Shape parameter Mean propagation model parameter Coefficient of hydrodynamic nonlinearities Positive phase duration Ratio of specific heats
	t a u δ th u δ X i w ǫ w κ µ	. A.2 Overpressure values and corresponding damage to structures from TNO [1989] 170 170 A.3 Overpressure values and corresponding injuries to people from various sources 170 Shock arrival time Sound diffusivity Flow velocity in the x-direction Coefficient of variation of random variable X i Urban environment parameter Scaling factor Flow velocity in the z-direction Thermal conductivity Stochastic model parameter Shear viscosity
	x µ B		Shock formation distance Bulk viscosity
	ξ			Forchheimer's nonlinearity parameter
	ρ ′ ρ 0	Density perturbation Latin uppercase letters Ambient medium density
	ρ T		Total density
	σ	D i	Random variable -space between buildings B i and B i+1 Nonlinear flow resistivity
	σ	E		Mathematical expectation Shock formation distance
	F (w) σ 0	Multi-objective function Linear, static flow resistivity
	G H σ 2 X i ψ	(1) 0	Green's function Hankel function of the first kind and of order zero Variance of random variable X i Grazing angle
	xvi xx xxii		xv xix xxi

. . . . . . . . . . . . . . . .

Table 1 :

 1 List of main industrial explosions since 2001. Various web sources were used to compile this list.

	Year Place	Company	Casualties and injured	Cost
	2008 USA	Georgia sugar refinery	13 deaths, 40 injured	$15.5 million
	2008 USA	Goodyear plant	1 death, 6 injured	unknown
	2007 USA	Little General Store	4 deaths	unknown
	2006 USA	Falk Corporation	3 deaths , 47 injured	$40 million
	2005 England Buncefield oil depot	no casualties	£873 million
	2004 Algeria	LNG liquefaction plant 27 deaths	unknown
	2004 Denmark Jest fireworks	1 death	unknown
	2001 France	AZF factory	30 deaths,	e2 billion
			several thousands injured
	Three types of effects are generally considered for industrial explosions.

Table 2 :

 2 Blast overpressure values and corresponding commonly accepted damage to structures. This table is taken from[START_REF] Merrifield | Simplified calculations of blast induced injuries and damage[END_REF].

	Overpressure (kPa) Damage
	0.1	Annoying noise, if of low frequency
	0.7	Breakage of small windows under strain
	1	Typical pressure for glass failure
	2.1	Damage to some ceilings
	2.8	Minor structural damage
	3.4-7	Large and small windows shattered, damage to window frames
	5.2	Minor damage to houses, 20-50 % of tiles displaced
	6.3	Roof damage to oil storage tanks
	7	Houses made uninhabitable
	7-14	Asbestos cladding shattered, fastenings of corrugated steel panels fail,
		tiled roof lifted and displaced
	14	Partial collapse of walls and roofs of houses, 30 % of trees blown down
	21	90 % of trees blown down, steel framed buildings distorted
		and pulled away from foundations
	21-28	Rupture of oil tanks
	28-35	Severe displacement of motor vehicles
	35	Wooden utility poles snapped
	49-63	Collapse of steel girder framed buildings
	56-70	Brick walls completely demolished
	>70	Complete destruction of all un-reinforced buildings
	living environment.

Table 4 :

 4 Summary of the overall cost of the Buncefield incident, by main category. This table is taken from Buncefield 2008.

	Sector

  Blast wave propagation . . . . . . . . . . . . . . . . . . . 28 1.4 Conclusions and outline of the document . . . . . . . . . . . . . . 34

1.1 Features of sound propagation outdoors . . . . . . . . . . . . . . 12 1.2 Nonlinear effects in sound propagation . . . . . . . . . . . . . . 22 1.3

Table 1 .

 1 2: Overpressure values at 200 m and 1000 m from the source for an explosion of 30 T of TNT calculated with the Kinney-Graham (KG) model.

	Quantity	Distance [m] Value [kPa]
	Positive overpressure	200	18.9
	Negative overpressure	200	-5.7
	Positive overpressure	1000	2.6
	Negative overpressure	1000	-0.236
	Table 1.3: Phases durations at 200 m and 1000 m from the source for an explosion of 30 T
	of TNT calculated with the KG model.
	Quantity	Distance [m] Value [ms]
	Positive phase duration	200	2.9
	Negative phase duration	200	34
	Positive phase duration	1000	4.1
	Negative phase duration	1000	91

  first step to derive the NPE model is to reduce Eqs. (2.2) to a single equation. Eq. (2.2a) is differentiated with respect to time and Eqs. (2.2b, 2.2c) are substituted in this new expression, yielding a fully nonlinear wave equation, such that

  [START_REF] Ambrosiano | Nonlinear propagation in an ocean acoustic waveguide[END_REF] and[START_REF] Castor | Long-range propagation of finite-amplitude acoustic waves in an ocean acoustic waveguide[END_REF] used the NPE to study high-amplitude wave propagation in shallow-water waveguides. More specifically the effects of nonlinearities in the propagation, transmission and reflection of waves at a fluid-fluid interface were under interest. The NPE was used to compare linear and nonlinear propagation. It appeared that hydrodynamic nonlinearities have significant effects:◮ the critical grazing angle at the ocean bottom is reduced, resulting in an increased transmission of energy to the bottom, ◮ and due to shock formation, the energy loss near the source is enhanced compared to linear propagation (see Section 1.2.2).In this chapter we detailed the existing work on the Nonlinear Parabolic Equation (NPE) model. Two different derivations are found in the literature; the first uses the perturbation expansion method (see Section 2.1.2) while the second makes some assumptions about the particles motion (see Appendix D). The NPE model, closely related to similar propagation models (the frequency-domain PE, the Burgers equation and the KZK equation), makes use of three assumptions: propagation along a main direction, weak sound speed perturbations and weak nonlinearities. The extension and variants of the model were then presented in

	"No satisfactory means has yet been found to include the absorption provided by
	porous ground surfaces."
	2.7 Chapter summary and conclusions
	Sections 2.2 and 2.3.
	Piacsek [2002] used a NPE model to study the effects of atmospheric turbulence on sonic
	boom wave fronts.
	The NPE has also been used to study propagation of explosion waves [see van der Eerden &
	Védy, 2005; Attenborough et al. , 2005, 2004]. The waves propagated had an amplitude of
	approximately 5 kPa at a distance of two meters from the source. Van der Eerden & Védy
	[2005] demonstrated the possibility of coupling the NPE model to other propagation methods

The NPE was then used to identify eigenmodes of the waveguides which have been compared to normal mode calculations. For the linear regime, simulated profiles agreed with mode predictions. Finally, to quote Ambrosiano et al. 's paper: "The consistency of linear results with the existing literature, taken together with the physical plausibility of nonlinear results, lends credibility to the NPE as a suitable numerical/theoretical tool for studying nonlinear acoustic phenomena."

2.6.2 Environmental acoustics

Blanc-Benon et al.

[2002]

investigated the propagation in finite-amplitude sound through turbulence. Two different models were used: the first is based on a geometrical acoustics approach, in which a nonlinear transport equation is used, while the second is a KZK equation modified to include a random temperature field (as explained in Section 2.1.3, KZK equations and NPE models are equivalent, with role of distance and time reversed). Running several numerical simulations allowed statistics to be calculated. The KZK equation was shown to be in agreement with the geometrical acoustics method: location of caustics coincides with the maximum peak pressure locations calculated with the KZK model.

(see Section 2.5). The NPE was successfully used to propagate explosion waves up to several hundred meters. The authors noted that:

The model can easily be coupled to other near-field or far-field propagation methods. Moreover, previous studies tends to indicate that simulation results obtained thanks to the NPE model are in good agreement with more complex models or measurements.

The NPE model is promising in several ways. Its relative simplicity (single-variable, one-way wave equation) supposes low computational times while its modularity indicates straightforward implementation and manipulation of the model (numerical implementation of the NPE model is detailed in Chapter 4).

However, note that the model lacks some essential features of outdoor sound propagation applications, such as the possibility of taking into account the ground topography and the ground impedance. This point is treated in the next chapter.

  Relation to the generalized terrain parabolic equation (GT-PE): the NPE model Eq. (2.19) has been proved to be the time-domain counterpart of the frequency-domain PE (see Section 2.1.3 and McDonald & Kuperman

  Treatment of initial conditions . . . . . . . . . . . . . . . . . 93 4.5 Treatment of boundary conditions . . . . . . . . . . . . . . . . 94 4.6 Some aspects of software development . . . . . . . . . . . . . . . 102 4.7 Chapter summary and conclusions . . . . . . . . . . . . . . . . 102

			84
	4.2 Calculation method for nonlinear terms with the Flux Corrected Transport algorithm 85
	4.3 Calculation method for linear terms with a semi-implicit scheme	. . . . . .	89
	4.4

4.1 The operator splitting method . . . . . . . . . . . . . . . . .

Table 5 .

 5 1: Ground topography characteristics.

	Quantity	Maximum value Location [m]
	max |h (x)| max h ′ (x)	5 m 0.22	70 m (and 78.8 m) 61.2 m
	max h ′′ (x) 0.01 m -1	55 m (and 85 m)

Table 5 .

 5 2: Reflected wave characteristics for reference and NPE calculations.

	Model t a+	p +	t d	t a-	p -
		[ms] [Pa] [ms] [ms] [Pa]
	Rigid layer:
	Euler 37.9 756 1.4 40.2 -449
	NPE	38.0 749 1.3 40.0 -438
	σ 0 = 100 kPa.s.m -2 :
	Euler 38.1 387 1.8 40.3 -265
	NPE	38.2 411 1.6 40.2 -274
	σ 0 = 10 kPa.s.m -2 :
	Euler 39.3 202 2.4	-	-
	NPE	39	213 2.0	-	-

As a mean of comparison, calculation times for Euler and NPE models were about 3.5 hours and 4 minutes, respectively (calculations were done on a modern desktop computer). Although the Euler equations implementation could use more advanced numerical techniques (AMR methods, see

Table 5 .

 5 4: Differences in reflected waves characteristics with and without Forchheimer's nonlinearities. Results are shown for both NPE and reference calculations.

	Model t a+	p +	t d
		[ms]	[Pa]	[ms]
	Euler	-0.7 +21.78 % -0.3
	NPE	-0.6 +21.12 % -0.3

  Reference model: linear propagation over urban cities . . . . . . . . . . 137 6.5 Mean nonlinear parabolic propagation model for sound propagation over urban cities 138 6.6 Construction of the probabilistic model of random variables Γ, Λ and Θ . . . . 142 6.7 Stochastic Nonlinear Parabolic Equation model for high-amplitude wave propagation over urban cities . . . . . . . . . . . . . . . . . . . . . 144 6.8 Identification of parameter w of the stochastic model . . . . . . . . . . 145 6.9 Application and experimental validation . . . . . . . . . . . . . . 153 6.10 Chapter summary and conclusions . . . . . . . . . . . . . . . . 158

	6.1 Summary of previous chapters and introduction	. . . . . . . . . . . 132
	6.2 Construction of the computational model: principles and methodology	. . . . 132
	6.3 Prior probabilistic model of geometrical parameters of a urban city	. . . . . 134
	6.4

  213 , 0.984 , 5.136 , 0.140 , 0.023 , 0.143) ,

		(6.60)
	for u 1 and
	w opt 2 = m opt Γ , m opt Λ , m opt Θ , σ opt Γ , σ opt Λ , σ opt Θ
	= (0.262 , 1.776 , 20.575 , 0.132 , 0.091 , 5.036) ,	(6.61)
	for u 2 .

  TableA.2: Overpressure values and corresponding damage to structures from TNO[1989].[START_REF] Kinney | Explosive shocks in air[END_REF] derived analytical solutions for the propagation of blast waves for spherical charges in free field. The model gives the pressure history p (t) at a given reduced distance in an homogeneous medium, for an equivalent charge W of TNT.

	1-7	Breakage of windows	B
	3	Light damage to structures
	10	Destruction of window frames Kinney-Graham model
	14	Partial collapse of walls and tiles
	16	Lower limit for heavy damage to structures
	17	50% deterioration of brick houses
	25	Destruction of light buildings and storage tanks
	50	Total destruction of houses
	70	Total destruction of reinforced buildings
	Overpressure (kPa) Damage
	83	Total destruction of structures
	35	Heavy damage to structures
	17	Moderate damage to structures
	3.5	Light damage to structures
	Table A.3: Overpressure values and corresponding injuries to people from various sources
	Effects	Overpressure (kPa) Reference
	Lethal (lung hemorrhage, 99%) 70	TNO [1989]
	Eardrum rupture	34	TNO [1989]
	Lethal (indirect effects)	14	Lees [1996]
			Baker et al. [1983]
	Irreversible damage	5	Lees [1996]

* This charge corresponds to the median value given by experts during the inquiry concerning the AZF factory explosion.

* The expression "FCT method" is found in papers byvan der Eerden & Védy [2005];[START_REF] Attenborough | Overview of the theoretical development and experimental validation of blast sound-absorbing surfaces[END_REF][START_REF] Attenborough | Blast sound absorbing surfaces[END_REF] and is used to describe a numerical implementation of Euler equations. The FCT method not being restricted to Euler equations implementations, "FCT method" should in fact be understood as: "a numerical implementation of the Euler equations using the FCT algorithm."

* Note that the main propagation direction is now the transformed x-direction.

* Note that the same conclusion applies to thermoviscous effects, added to the GT-NPE model in Section 3.1.3

* see[START_REF] Hamilton | Nonlinear acoustics[END_REF], Chapter 4, Section 2.4.1.

* in English Integral Methods for Acoustic Calculations of Diffraction by Obstacles

Acknowledgments -Remerciements-

Appendix F: The Boundary Element Method 185

References and Section 6.2 details the methodology for constructing the computational model. Next the probability models of the geometrical parameters of a urban environment are determined using Information Theory and the Maximum Entropy Principle (Section 6.3) and then the reference model is described (Section 6.4). The mean propagation model is then presented in Section 6.5. In Sections 6.6 and 6.7 the construction of the probabilistic models of the stochastic model parameters is detailed and in Section 6.8 the method used to identify its parameters is presented. Section 6.9 details two applications using the constructed stochastic propagation model and finally, chapter summary and conclusions are given in Section 6.10.

General conclusions and perspectives of this work are finally given in page 165.

Five Appendices follow the general conclusions. Appendix A presents blast wave overpressure values and corresponding damage to structures and injuries to people from different sources.

Appendix B details the Kinney-Graham (KG) References can be found in page 187. A long summary in French and a list of personal publications appear in pages xxvii (201) and xxxv (209).

A NPE model for propagation within porous layers is then derived from a nonlinear extension of the Zwikker-Kosten (ZK) model. It is assumed that the ground layer is equivalent to a continuous fluid and with these hypothesis the model uses four parameters.

In order to develop a model for sound propagation over porous ground layers a boundary interface condition that allows one to couple the two NPE models is derived. These two-way coupling equations only involve spatial derivatives and integrals, making its implementation natural and straightforward. The terrain-following coordinate formulation of this model is then given together with expressions for multi-layered ground surfaces. that the method presented to take into account the flow resistivity dependence on particle velocity is accurate.

Summary of previous chapters and introduction

The objective of this chapter is to develop a computational model to simulate long-range nonlinear wave propagation over urban cities. We are interested in the effect of the surface irregularities (buildings) on the acoustic field above the urban layer.

This problem could be studied with deterministic numerical models where buildings geometries are explicitly accounted for. For example, ray tracing methods or models based on the Euler equations allow the environmental context (buildings, meteorological conditions, site topography, etc.) to be taken into account and would technically be suited for this application. However, for long-range propagation applications, these methods suffer from their numerical complexity and the high computational effort associated. Moreover, the suitability of these models for the application under interest can be questioned. Considering the high complexity of the real system, the model approximations (e.g. using the Uniform Theory of Diffraction (UTD) for ray-tracing methods) and the uncertainties on the model parameters (e.g. buildings geometries), the computational model could be improved introducing a probabilistic model.

In this work, a different approach is proposed. It consists in using a simplified model for the urban and atmospheric layers. Propagation is modeled with a Nonlinear Parabolic Equation (NPE) model. Since the calculation domain is limited to a small area around the signal, computational cost is generally reduced compared to Euler equations methods. The urban layer is modeled by a porous medium, in which a NPE formulation is used for analyzing shock wave propagation. Therefore such a model is composed of two domains:

◮ the atmospheric layer, where propagation is modeled by a NPE for air medium, ◮ and the urban city layer, where propagation is modeled by a NPE for porous medium.

Equations to couple these two domains are added. Using this model as a predictive model for the real system corresponds to a raw approximation of the propagation phenomenon. Hence a probabilistic approach of uncertainties is used to enhance the model capabilities. Section 6.2 explains and details the methodology to construct such a stochastic sound propagation model and the probabilistic model of its parameters, and outlines the content of the chapter.

Construction of the computational model: principles and methodology

The objective is to develop a stochastic model for sound propagation over urban cities using a NPE model originally designed for sound propagation over porous ground layers, in which sake of brevity the explicit expressions of Sommerfeld radiation conditions are not given here.

The propagation domain is truncated with the help of a Perfectly Matched Layer (PML) (see Section 4.5.3).

Wave propagation in the urban layer (Eq. (6.9b)) is characterized by parameter µ which modifies the sound speed in the urban layer so that c Ωu = µc Ωa , while α is the loss rate in the layer. The interface condition between domains Ω a and Ω u (Eq. (6.9c)) depends on parameters γ and ǫ. On the boundary Γ r a null pressure r a = 0, r u = 0 is imposed, meaning that no perturbation is introduced ahead of the wavefront.

The relative SPL at the receiver is given by

where quantities with hats denote Fourier-transformed quantities, and x r and z r are the receiver coordinates.

Construction of the mean model and description of the algebraic properties of its parameters

In this section the mean propagation model is derived from Section 6.5.1. One has to define the parameters in the propagation model (see Eqs. (6.9)) which have the capability to represent the natural variability of the real system (the urban environment). These defined parameters will be modeled by random variables as explained in Section 6.2.

When a wave is reflecting on a plane surface, its amplitude is changed and a (possibly negative) delay is given to the reflected wave. Parameters (µ, α, γ, ǫ) appearing in Eqs. (6.9b, 6.9c) could be used because the propagation model can recreate the behavior of a porous ground layer [see Leissing et al. , 2009c] and could thus be used "as is" to control the wave reflection.

However, in order to reduce the number of parameters in the mean model and consequently, to reduce the stochastic model complexity, another solution is proposed. The urban layer is now considered as semi-infinite, and one writes that no waves are transmitted back from the porous ground layer to the atmospheric layer. This implies that the urban layer must behave like an atmospheric layer, which can be simulated in writing that µ = 1 and α = 0 in Eq. (6.9b).

The time delay occurring during reflection is introduced by the use of a time-stretching. The time variable is then transformed according to (6.11) envelopes:

in which fix (x) is the integer part of the real number x.

Now introducing the dependence on parameter w, the formal definition of the mean-square method now follows. Let ℓ + (w, ω) and ℓ -(w, ω) be the upper and lower envelopes of the confidence region of the stochastic model. The functions ℓ +exp (ω), ℓ -exp (ω), z + (w , ω) and

and (6.33b) in which Heav (x) is the Heaviside function such that Heav (x) = 1 if x ≥ 0 and Heav (x) = 0 otherwise. Functions z + (w, ω) and z -(w, ω) represent selected parts of ℓ + (w, ω) -ℓ +exp (ω)

and ℓ -(w, ω)-ℓ -exp (ω) where the experimental observations do not belong to the confidence region calculated with the stochastic model.

The non-differentiable objective function J (w) is then defined by (6.34) and the optimal parameter w opt is solution of the following optimization problem:

The mean-square method hence aims at minimizing the areas where the experimental observations do not belong to the confidence region of the stochastic model: the only criteria to select the optimal parameter w opt is the amount of information not covered by the stochastic model.

Multi-objective optimization using an evolutionary algorithm

In this work the Non-dominated Sorting in Genetic Algorithms (NSGA-II) method developed by [START_REF] Srinivas | Multiobjective optimization using nondominated sorting in References genetic algorithms[END_REF]; [START_REF] Deb | A Fast Elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF] is used to solve the multi-objective optimization problem. This method is part of the family of genetic algorithms, in which an initial population evolves over several generations of individuals. For each generation a selection process selects the "best" parents from which children are generated. The NSGA-II method is composed of five major steps, summarized below.

Step 1, population initialization. The population is initialized based on the results of Section 6.8.2. The mean model updated parameters are used as initial values for the mean values m Γ , m Λ and m Θ and a 50 % variation of these quantities is allowed, i.e.

The constraints on the dispersion of the random variables Γ, Λ and Θ are arbitrarily set to σ Γ ∈ ] 0, 0.25 [ , (6.48a)

Eqs. (6.47, 6.48) define a reduced admissible set W red for the parameter w, defined by

The initial population of n i individuals is hence generated by calculating n i families of ν observations from the stochastic model, where the parameters m Γ , m Λ , m Θ , σ Γ , σ Λ and σ Θ are uniformly distributed on their respective domains defined by

Step 4. and finally, using the stochastic parabolic propagation model derived in Section 6.7, to construct the confidence region associated with the optimal parameter w opt determined in Step 3 which finally allows the stochastic propagation model to be validated (see Section 6.9.4).

6.9.2 Choice of parameter u, numerical experiment description and output from the reference model

As explained in Section 6.4, the reference model is defined as a numerical experiment consisting of two-dimensional numerical calculations using the Boundary Element Method (BEM). In this application two different city morphologies are chosen. The corresponding parameters are u 1 and u 2 such that

= (10, 20, 30, 0.2, 0.2, 0.2) , (6.55)

and

= (40, 40, 30, 0.2, 0.2, 0.2) . (6.56)

Once u 1 and u 2 are fixed the buildings dimensions are generated with the help of the probability distributions defined in Section 6.3. The number of realizations of the probabilistic model of the real urban city is ν exp = 500, inducing a convergence with a 10 -3 accuracy (at least) for the first and second-order moment of the random variable L exp . A calculation with the nominal mean model described in Section 6.5 is also performed for comparison. It is initialized with 6.9.3 Identification of the stochastic model optimal parameter w opt In this section, the procedure described in Section 6.8 is used to identify the optimal parameters w opt 1 and w opt 2 of the stochastic model corresponding to the parameters u 1 and u 2 .

Determination of the updated parameter of the mean model

In the identification of the optimal parameter w opt of the stochastic model, the first step is the determination of the updated parameter s opt of the mean model (Step 1 in Section 6.8.2), solution of the minimization problem defined in Eq. (6.27).

The effects of parameters γ , λ and θ on the reflected wave being uncoupled, the minimization problem is reduced to three one-dimensional searches. The parameters are real and we use * a parabolic interpolation and the Brent's method [see for example Press et al. , 1996b]. The optimal parameters found from this minimization problem are

Identification of the optimal parameter of the stochastic model with a genetic algorithm

The updated parameters of the mean model s opt 1 and s opt 2 are then used to define the reduced admissible set W red of w, which is used to initialize the genetic algorithm. For the evaluation of the mean-square norm the 0 Hz point is removed from the calculation and the probability level used for the construction of the confidence region is P c = 0.98. For the statistical reduction of information performed in the maximum likelihood method the order of decomposition q (see Eq. (6.39)) is q = 12. With this order of decomposition, the accuracy ǫ in Eq. (6.43) is such that ǫ < 10 -3 , meaning that at least 99.9 % of the information is contained within the random variables {Y 1 (w) , . . . , Y 12 (w)}.

Concerning the genetic algorithm, the number of individuals in the population is set to n i = 50, for each individual ν = 500 observations from the stochastic model are generated and the population evolved over 50 generations. At each generation, half the population is replaced with children. The parameters η c and η m in Eqs. (6.52, 6.54), which control how well spread are the children from the parents, are both set to 20, ensuring a diversified population. Fig. (6.8) shows the values of the objective functions -L red (w) and J (w) of the 50 individuals at generation 1, 10 and 50 (the final generation) for both parameters u 1 and u 2 . It can be Once the probability distributions are determined the stochastic model parameter is identified with an evolutionary algorithm. This inverse stochastic problem involves two different methods: the mean-square method, which aims at minimizing the areas where the experimental observations do not belong to the confidence regions of the stochastic model, and the maximum likelihood method, which aims at maximizing the likelihood between the experimental realizations and the observations from the stochastic model. In order to reduce the computational effort associated with the standard maximum likelihood method, a statistical reduction of information is performed. It consists in a principal component analysis of the observations. This allows one to use the maximum likelihood method in the space of reduced random variables, which are by construction uncorrelated, and hence do not require the fastidious evaluation of their joint probability distribution. So as to provide a robust and accurate identification method, the mean-square method and the maximum likelihood method are then used concurrently in an evolutionary algorithm. It consist of initializing a population of individuals, some of which being selected based on non-domination to be the parents of the next generation. The initial population is generated with the help of the mean model. It was shown that the evolutionary algorithm provides robustness: even if the initial population exhibits poor performances, after only 10 generations one can obtain almost optimal individuals.

Within the final generation one can find well diversified individuals: the evolutionary algorithm returns individuals that optimize one objective at the expense of the second one. The decisive criteria to select the final individual was chosen to be the mean-square norm. This way the least possible underestimation of experimental observations is achieved. Once the optimal parameter w opt of the stochastic model is determined it is used to solve the stochastic equations with the Monte Carlo method. Confidence regions are then constructed and compared to experimental observations. It was shown that the observations fall within the confidence regions, except at the very lowest frequencies, independently of the amount of dispersion in the real system. To put it simply, the highly complex surface that represent the urban environment was successfully replaced by a stochastic but flat porous surface, which allows statistics on the acoustic fields above this layer to be obtained.

The constructed model could thus be used to study nonlinear wave propagation in complex environments. Features of outdoor sound propagation, such as dissipation effects or refraction effects, are naturally present in the NPE model and can be incorporated in the stochastic model. Thanks to the low numerical effort associated with this model, large parametric studies could be performed, including most of the features of sound propagation outdoors.

Synthesis of the research and perspectives

Synthesis of the research

The objective of this research was to develop a numerical tool for the study of blast wave propagation in complex environments and over urban environments with application to the risk assessment of industrial installations. The solution proposed in this document was the following. Instead of using time consuming numerical methods it was chosen to develop a simplified model. This mean propagation model allows one to get a raw approximation of the propagation phenomenon. To enhance the model capabilities, it is associated a probabilistic model of uncertainties to take into account the real system complexity, yielding a stochastic propagation model.

The deterministic simplified model is based on the Nonlinear Parabolic Equation (NPE) model. The assumptions used for the derivation of this model are: weak nonlinearities, weak meteorological perturbations and a paraxial approximation. The original model is unable to take into account ground topography or ground impedance.

The deterministic model is then extended to handle various features of sound propagation outdoors. The terrain-following coordinates method is used to take into account ground topographies. This model allows to simulate sound propagation over hilly urban areas. A second NPE model is then developed for propagation within porous ground layers. Associated with coupling equations, the complete NPE model allows wave propagation within a moving atmosphere, over a non-flat, non-rigid ground surface to be simulated.

Perspectives ated with the standard maximum likelihood method, a statistical reduction of information is performed. So as to provide a robust and accurate identification method, the mean-square method and the maximum likelihood method are then used concurrently in an evolutionary algorithm. Once the optimal parameter of the stochastic model is determined it is used to solve the stochastic equations with the Monte-Carlo method. It was shown that the observations fall within the confidence regions, except at the very lowest frequencies, independently of the amount of dispersion in the real system. The complex surface that represents the urban environment was successfully replaced by a stochastic but flat porous surface.

Perspectives

On the deterministic NPE model. The deterministic NPE model used is based on existing work presented in Chapter 2, and is further extended in Chapter 3. This model can handle most of the features of sound propagation outdoors and can thus be used to study shock wave propagation in complex media, including the effects of meteorological conditions, topography and ground impedance, or thermoviscous effects. It hence provides a complete, accurate and fast simulation tool, perfectly suited for the study of, for example, sonic boom propagation, infra sound propagation, evaluation of noise from military tests, etc.

On the probabilistic approach of uncertainties. When comes the choice of a numerical model for the study of a particular problem, one ideally chooses the most accurate model with respect to the physical phenomenon involved. In practice, one must usually accept a trade-off between the model accuracy and its complexity (related to its numerical implementation, computational times, cost, etc.), so that the choice of a numerical model is often reduced to the question "What trade-off is acceptable? ". Without references (measured data or solutions from a reference numerical model) or a solid experience of a specific numerical model, assessing the model accuracy is a difficult task [Wilson & Petit, 2009].

We used in this work a quite different approach from classical sound propagation models.

We voluntarily choose a (too!) simple model for the physical phenomenon involved, so that the trade-off would not be acceptable. For example, back-scattering between buildings, which plays an important role in the propagation process, is totally absent from the NPE model. This physical phenomena and others (and their coupling) are all accounted for in the probabilistic model of uncertainties. This method for designing numerical models can be applied to a wide varieties of physical problems, and most of the time provides very rich and interesting answers. The key question is now "How well can the probabilistic model account for the real system complexity? ". Considering that the probabilistic model of uncertainties is not restricted by physical parameters, it is often possible to derive sufficiently robust models.

On the stochastic NPE model for the study of blast wave propagation. This research corresponds to a first step in which the propagation models are developed in 2D.

The main perspective is the extension of the different models in 3D. The method presented in this document is very general and nothing prevents this extension to be performed. However, computations for three dimensional cases seem a quite computationally intensive task. In this case, in order to obtain reference solutions, measured data or measurements on scale models should be preferred instead of numerical calculations with the BEM.

The simulation tool developed could be used to perform risk assessment of industrial installations. A complete two-dimensional stochastic simulation, once the reference solutions are obtained and the optimal parameter is identified, takes about an hour to complete on one CPU. This would enable the use of the model for overpressure mapping. A complete map could be obtained in about two days (about 50 compass directions). Moreover, the information on the maps could be improved in two ways.

1. First, "conditional maps" could be obtained. Indeed, as we saw in Chapter 1, meteorological effects can yield strong shadow zones, or can produce a double echo. Taking into account this feature would allow one to produce maps related to the local (in space and time) meteorological conditions.

2. The model could be used to produce "statistical maps". Considering the high variability of the system considered one may wonder what the induced variability on the overpressure levels is. The stochastic model allows one to provide a maximum overpressure value with an associated probability level, which is, in the frame of risk assessment, and in the author's opinion, much more pertinent than a single overpressure value.

Concerning the practical use of the models proposed in this work, an improvement would be the construction of the function (polynomial expression, neural network approach, etc.) that relates the optimal parameter w opt to the urban city geometrical parameters u.

Future research with the stochastic NPE model. The stochastic NPE model can be used to lead some research on wave propagation in random media. Specifically, the acoustical consequences of propagation through regularly or randomly spaced obstacles can be observed.

The model can also be used to study wave scattering with respect to the scatterers properties (density, characteristic size, shape, etc.). Among propagation in urban environments, room acoustics, and specifically the study of acoustic diffusers, is a field where the stochastic sound propagation model developed could provide interesting solutions.

A

Blast wave overpressure values and corresponding damage to structures and injuries to people, from different sources

This appendix presents figures found in the literature for some overpressure values and the corresponding damage to constructions and injuries to people from different sources.

A.1 Overpressure values and corresponding damage to structures Tabs. (A.1, A.2) presents overpressure values and corresponding damage to structures found in [START_REF] Lannoy | Analyse des explosions air-hydrocarbure en milieu libre : Etudes déterministes et probabilistes du scénario d'accident[END_REF]; [START_REF] Lees | Loss Prevention in the Process Industries[END_REF] and TNO [1989].

A.2 Overpressure values and corresponding injuries

Tab. (A.3) presents overpressure values and corresponding injuries found in [START_REF] Lees | Loss Prevention in the Process Industries[END_REF]; TNO [1989]; [START_REF] Baker | Explosion Hazards and Evaluation[END_REF].

is determined with the help of the impulse per unit surface I/A, defined by

With this set of equations, time waveforms (and hence spatial waveforms) can be calculated.

Examples of blast waveforms calculations using the Kinney-Graham (KG) model can be found in Section 1.3.6. and p T = ρ T RT , (C.4) in which T is the gas temperature, c v is the specific heat capacity at constant volume and R is the gas constant. To solve this equation system a Weighted Essentially Non-Oscillatory (WENO) algorithm [see [START_REF] Shu | Essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF] for space discretization and a third-order Total Variation Diminishing (TVD) scheme [see Press et al. , 1996a] for time marching are used.

The WENO scheme has fifth-order accuracy. The main principle of the WENO scheme is the use of multiple stencils to evaluate the derivative at a given point. The algorithm first determines where there is a discontinuity and then weights stencils accordingly to avoid spurious numerical oscillations. This features make the WENO scheme accurate for propagating shock waves.

The time discretization scheme is of the form It is hence an average between fully explicit and fully implicit models of PDE's. This is where the second order convergence comes from. This yields a tridiagonal linear equation system that can quickly be solved with a Thomas algorithm (see below).

E.2 Thomas algorithm

Semi-implicit schemes such as the Crank-Nicolson scheme require a linear system of equations to be solved. If centered first-order finite-difference approximations of the derivatives are used then the left-hand-side matrix is tridiagonal. The Thomas algorithm is an efficient method

The implementation of the BEM used in this work (MICADO) is based on the variational approach. The geometry of the problem is two-dimensional : the source is an infinite coherent line of sources and the configuration remains unchanged in the direction perpendicular to the vertical plane.