David Nouvel

Sébastien Travadel

Erik Hollnagel

published or not. The documents may come   L'archive ouverte pluridisciplinaire

The analysis of complex events shows that the performance of crew members and/or air traffic controllers, as well as the aircraft and its systems, may vary within different time frames. Thus, the understanding of the current situation can differ significantly within the crew and between the cockpit and the ground, as well as be at odds with the aircraft's real position and attitude over time. Further, these components interact throughout the flight.

This paper presents a proposal for how better to model these complex interdependencies, using the principles of the Functional Resonance Accident Model (FRAM). The outcomes of this modelling can conveniently be shown by means of open source software, the FRAM Visualizer.

Using an incident initiated where the Captain selected the wrong HSI track, we explore the many factors that influenced the situation when the crew put the aircraft into descent but failed to intercept the localizer. This led to a crisis that they had to manage, departing significantly from Standard Procedures, and without being fully aware of their proximity to the ground when they initiated a go-around.

Our aim is to develop new analytical tools in order to better characterize and understand how the flight evolved beyond acceptable safety margins, and how the crew recovered from a critical situation. Some of the factors and their comparative interactions identified in this manner may reveal inherent risks that are still present but also preventable.

Keywords: Functional resonance analysis, aircraft accident investigation,

The BEA investigated an incident that occurred on 23 November 1997 on approach to Paris-Orly to the MD83 registered F-GRMC. A final report 1 was issued at the time. It was determined that the incident resulted from the decision to put the aircraft into descent when, as a result of a display error, it was neither on the localizer track nor on the glide path, and with no context defined for this improvised manoeuvre. The operator's company culture directly contributed to the incident through the importance it attached to accelerated training given to new First Officer and to undertaking commercial flights. Ten years later, a re-

The Functional Resonance Accident Model

The functional resonance accident model (FRAM; [START_REF] Hollnagel | Barriers and accident prevention[END_REF] describes system failure in terms of the resonance of normal performance variability. This provides a convenient way of representing the non-linear propagation of events and also makes it possible to account for adverse outcomes in cases where there were no manifest malfunctions or failures. The principle of FRAM is to characterise individual system functions independently of how they may be connected in a specific situation. The characterisation of each function -or node -is done in terms of six aspects and the values of these aspects determine how nodes may be coupled under given conditions. To produce a description of functional variability and potential resonance, and to determine recommendations for damping unwanted variability, a FRAM analysis consists of four steps:

Step 1: Identify essential system functions, and characterise each function by six basic aspects or parameters. The six aspects are input (I, that which the function uses or transforms), output (O, that which the function produces), preconditions (P, conditions that must be fulfilled to perform a function), resources (R, that which the function needs or consumes), time (T, that which affects time availability), and control (C, that which supervises or adjusts the function).

Nodes and their aspects may be described in a table and can subsequently visualized in a hexagonal representation (cf. Figure 2 and 4 below).

Step 2: Characterize the context dependent variability of each node. For an accident analysis, the variability is known from the investigation data. In this case the analysis focuses on comparing the observed and the normal performance. For risk assessment, the variability may be derived from a characterisation of the common performance conditions (CPCs), of which there currently are eleven. These CPCs address the combined human, technological, and organizational aspects of each function. After identifying the CPCs, the variability must be determined in a qualitative way in terms of stability, predictability, sufficiency, and boundaries of performance.

Step 3: Defining the functional resonance based on possible dependencies/couplings among functions and the potential for functional variability. The output of the functional description of step 1 is a characterisation of functions and their aspects. The aspects provides the basis for identifying how functions may be coupled. For example, the output of one function may be an input to another function, or produce a resource, fulfil a pre-condition, or enforce a control or time constraint. When the couplings between functions are found, this is combined with the characterization of performance variability from Step 2. In this way the analysis will show how the variability of one function may have an impact on others. This analysis thus determines how resonance, and ultimately adverse outcomes, can result from variability across functions in the system. For example, if the output of a function is unpredictably variable, another function that depends on this output as a resource may be performed unpredictably. Many such occurrences and propagations of variability may have an effect like resonance.

Step 4: Identify barriers for variability (damping factors) and specify required performance monitoring. Barriers are hindrances that may either prevent an unwanted event from taking place, or protect against the consequences [START_REF] Hollnagel | Barriers and accident prevention[END_REF]. Barriers can be described in terms of barrier systems (the organizational and/or physical structure of the barrier) and barrier functions (the manner by which the barrier achieves its purpose). In FRAM, four categories of barrier systems are identified (each with their potential barrier functions, see [START_REF] Hollnagel | Barriers and accident prevention[END_REF]. In addition to recommendations for barriers, FRAM can also be used to specify recommendations for the monitoring of performance and variability, as a way to detect and manage undesired variability. Performance indicators may thus be developed for every function and every link between functions.

Example of application to a complex incident

Summary of the incident

During a commercial flight to Paris-Orly, the Captain decided to perform an automatic landing due to the weather conditions at destination. Since the First Officer was not qualified for this procedure, the Captain took over most of the crew actions. On approach, he selected a wrong approach track on the HSI. Due to this error, the auto-pilot did not intercept the landing track (localizer), and the aircraft departed from the published trajectory. The Captain tried to recover from this situation by selecting, unsuccessfully, various modes on the auto-pilot. The aircraft crossed the glide path and twenty seconds later, the Captain corrected the selection of the track on the HSI. At the same time, he put the aircraft into descent, heading simultaneously towards the approach track. The auto-pilot was not able to capture the glide path, so that the aircraft continued its descent without external references. The Captain did not manage to stabilize his approach, and eventually initiated a go-around. The lowest height reached by the aircraft was 67 feet from the ground.

Description of the event using the FRAM Visualizer

The FRAM Visualizer is a tool which is used to describe the event as a set of functions and their associated control parameters. A node defines the basic information unit. There is no scale restriction to describe a function. In this example, it was decided to take the elementary action, such as the action of positioning the index on the HSI, to define the resolution of the scheme. This is in line with the pilot's sequencing of actions: action -check -correction.

However, the availability of the data impeded the appropriate breaking down of the sequence and led to a more macroscopic approach for some functions.

We considered mainly the incident flight from the start of radar vectoring by Orly Approach to the go-around performed by the Captain. Nevertheless, some preliminary information was relevant to set the initial circumstances that developed throughout this phase of flight. Indeed, first node corresponds to the "task-sharing decision".

fact that the Captain departed significantly from the standard procedure, and improvised a sequence of actions, challenged the usual examination of a flight through the identification of deficiencies as compared to a standard flight. The standard flight could not be the sole reference to break down the actions performed into a comprehensive set of nodes. This justified the choice of taking into account a mix of:

• standard procedures that were followed (high level of performance),

• or standard procedures that were not or poorly followed (low level of performance),

• or improvised actions that can be seen as an attempt to adapt to an unusual situation with more or less success.

Below is an example of a node corresponding to the action by the Captain of rotating the HSI index from the interception heading to the landing track. According to the operator's operating manual, this action is undertaken once the LOC CAP mode is displayed, and triggers a check by the First Officer. Two inputs were identified: one triggered the action (LOC CAP display) and the other was a parameter that would evolve throughout the action. This choice will be discussed in Section 3.1.

The outputs include a callout that should have triggered a check by the first officer as well as the incorrect positioning of the index. In order to refer to the resonance approach, it was decided to take into account the variation of the output, between a "nominal" state and a "downgraded" state. In this latter case, the corresponding output was qualified as: "wrong", "failed", etc. The selection of track 258 provided erroneous information to position the index. This information persisted during the action and was therefore considered as a precondition.

The control of the function can be seen as a of the consistency of the respective positions of the index and the needle in accordance with the heading. The proximity of the tail of the needle with the intended landing track was misleading when the Captain positioned, as a reflex, the index on the tail of the needle. A quick did not allow him to detect the discrepancy. However the reading of the heading led the Captain to feel that something was wrong. The figure below illustrates the display seen by the Captain, as compared to the one with the correct track selection. At this stage, it should be noted that the functional approach helps to question the role and the description of every action, such as the one specified in the operating manual. For instance, while elaborating this node, in connection with the others, it was noticed that the operating manual did not specify what this check by the First Officer should refer to and more generally what was the intention of this action.

Connections of the nodes

The FRAM Visualizer automatically connects the control parameters of the nodes. This function provides graphic information regarding the influence that one parameter can have on several nodes. The automatically generated visualization produced a rather intricate scheme. It was then decided to rearrange the nodes in a time-sequenced scheme as shown in the next figure. For this representation, we decided to portray the Captain's actions on the main track, and position the other nodes along side. The nodes are according to whether they refer to:

• crew actions,

• system functioning,

• aircraft motion,

• system display.

This representation intrinsically defines some "components", by analogy with dynamic systems.

Node "CPT index positioning"

Focus on three sequences of the incident highlighted by FRAM

It should be noted that it was not possible to collect any additional data or to undertake any simulation processes and, therefore, the statements that follow are only valuable for purposes of illustration but do not constitute, stricto sensu, new investigation findings. However, using the functional approach, the data was put together with a view to placing the actors back in their interactive context, in order to bring to light developed information. We provide hereafter three examples of sequences that refer to the Captain's perception and adaptation.

Captain's sequence awareness

The functional approach brought to light the fact that the Captain had a kind of sequence awareness, in the sense that his experienced partially counter-balanced his poor situational awareness. This is probably one of the factors that led him to alter the heading by forty degrees to capture the localizer track, after he crossed the approach track, whereas the situation displayed on his HSI would rather suggest that he keep the previous interception heading. Later on, the glide deviation triggered some reactions by the Captain. In particular, he called for gear extension as the aircraft was approaching the glide path. Even though this call came late in comparison the standard procedure, it shows that the Captain acted reactively. From this perspective, the fact that the aircraft crossed the glide path the first time probably increased the time pressure on the Captain. As he did not capture the localizer, he did not put the aircraft into descent at that time, nor did he call for the flaps. Instead, he focused his attention on his horizontal situation till he understood his mistake. When he noticed his wrong track selection, he tried to catch up with his vertical position. It would be interesting to deeper analyse the parameters that determine the Captain's sequence awareness. They certainly vary in importance according to the phase of the approach. In this respect, it's worth comparing the above-mentioned circumstances to the circumstances in which the Captain initiated the go-around.

Captain's understanding of the erroneous track selection

During the fifteen second sequence preceding the Captain's correction of the track selection, the aircraft was heading straight with wings level; there was no radio communication. The situation displayed on the HSI was "consistent" with an interception phase. The Captain armed the ILS mode and then expected the autopilot to capture the localizer. He was then, for the first time, in appropriate conditions to scrutinize his HSI and discover the discrepancy.

Captain's decision to go around

As the aircraft crossed the glide path for the second time, the "glide slope" warning sounded, followed by an announcement by the First Officer. The Captain disconnected the auto-pilot, turned left to intercept the landing track as the needle became alive on his HSI. He then reconnected the auto-pilot and, within a few seconds, disengaged it and initiated the go-around. According to his testimony, the Captain was not aware of his vertical position and decided to go-around because he did not feel stabilized. The use of FRAM suggests that this feeling may be influenced by the gap between the modes displayed on the FMA and the ones expected at this stage according to the operating manual. Indeed, none of the three modes expected were displayed on the FMA.

A contribution to further developments to FRAM

States, functions and definitions

The FRAM Visualizer was initially designed to describe functions, such as actions by the crew. However, it can be seen that, for instance, the node shown as an example in section 2 is a mix of an action, the Captain positioning the index, and a change of state, the rotation of the index. In this case, the understanding of the consequences of this function was not hampered. The situation becomes more challenging regarding the evolution of the HSI display in accordance to the aircraft position. Therefore, it appeared necessary to distinguish the functions that refer to a crew action from the functions referring to change of of systems such as the HSI or the autopilot mode announcement (FMA). From this perspective, it became necessary to define the control parameters of a node in accordance with its nature.

From an "action" perspective, an input is something which triggers an action that produces an output. The control parameter is understood as an immediate check, associated with an action. A resource-consuming control action, such as a check prescribed in a standard operating procedure, is described via a dedicated node. The precondition is understood as a contextual element that influenced the result of the function, that pre-existed to the function and that varied slowly during the function's processing. The fact that the crew experienced poor CRM was a precondition of most of the functions during the approach.

From the "state" perspective, an input is something processed by the system to produce the output. The precondition is then a necessary condition for the state to evolve. For instance, the fact that the localizer track be captured was a precondition for the capture of the glide by the autopilot. The periodic check by the automatic system of the situation -the non interception of the localizer track -did not allow for the successful change of state of the autopilot, from an armed mode to a track mode. From this perspective, the resource refers to the kind of energy required. Time is a parameter which requires further work to be taken into account. During this study, the associated tag was simply used to position a node in the sequence. Time should intervene as a control parameter of the variation of the function.

Multilayer analysis

The need to take into account the interactions between different components such as systems or crew members calls for a rearrangement of these components in various layers. From this perspective, each layer would describe the evolution in time (chain of functions or changes of states) of homogeneous components. It would be then necessary to represent the interactions between the layers. This kind of distinction would facilitate comprehension of the interactions, as well as the identification of the resonance phenomenon. Indeed, components of the same nature tend to evolve on a similar time scale. For instance, the aircraft followed a quite simple track during the approach, marked by three changes of heading; at the same time, the FMA displayed a multitude of different modes.

Conclusion

The functional approach should lead us to question every action performed or missed. In this sense, the development of this method should help with the collection of data in the course of an investigation. Moreover, it highlights the deficiencies of a function itself, as well as those

Figure 1 :

 1 Figure 1: Horizontal (upper graph) and vertical (lower graph) tracks.

Figure 2 :

 2 Figure 2: Description of a node by FRAM Visualizer.

Figure 3 :

 3 Figure 3: Left side: incorrect positioning of the index on the tail of the needle, taking into account the erroneous track selection (258°). Right side: intended position of the index on the head of the needle, with the correct track selection (065°)

  Figure 4: Automatically generated visualization by FRAM Visualizer

Figure 5 :

 5 Figure 5: Horizontal track (blue dots) and FRAM nodes

Introduction of the Concept of Functional Resonance in the Analysis of a Near-Accident in Aviation.

	David NOUVEL, Sébastien TRAVADEL
	Bureau d'Enquêtes et d'Analyses pour la sécurité de l'aviation civile
	Statistics and Safety Analysis Division
	Bâtiment 153, Aéroport du Bourget
	93352, LE BOURGET Cedex, France
	Erik HOLLNAGEL
	École nationale supérieure des Mines de Paris
	Pôle Cindyniques
	Rue Claude Daunesse
	F-06904, Sophia Antipolis, France

33rd ESReDA Seminar: Future challenges of accident investigation

related to the immediate context of an action. The connection of the nodes through their control parameters shows the evolution of this context.

However, the analysis must further consider systemic deficiencies. How did an adverse outcome propagate? How were its effects dampened or heightened? To answer these questions in relation to a complex environment, it appears necessary to pursue the development of the analytical tools, and to better describe the evolution over time of the interactions. The challenge is to point out some necessary factors that condition the stability of the standard "flight structure".


J F Semblat

Modeling seismic wave propagation and amplification in 1D/2D/3D linear and nonlinear unbounded media

Keywords: earthquake engineering, engineering seismology, wave propagation, numerical modeling

To analyze seismic wave propagation in geological structures, it is possible to consider various numerical approaches: the finite difference method, the spectral element method, the boundary element method, the finite element method, the finite volume method, etc. All these methods have various advantages and drawbacks. The amplification of seismic waves in surface soil layers is mainly due to the velocity contrast between these layers and, possibly, to topographic effects around crests and hills. The influence of the geometry of alluvial basins on the amplification process is also know to be large. Nevertheless, strong heterogeneities and complex geometries are not easy to take into account with all numerical methods. 2D/3D models are needed in many situations and the efficiency/accuracy of the numerical methods in such cases is in question. Furthermore, the radiation conditions at infinity are not easy to handle with finite differences or finite/spectral elements whereas it is explicitely accounted in the Boundary Element Method. Various absorbing layer methods (e.g. F-PML, M-PML) were recently proposed to attenuate the spurious wave reflections especially in some difficult cases such as shallow numerical models or grazing incidences. Finally, strong earthquakes involve nonlinear effects in surficial soil layers. To model strong ground motion, it is thus necessary to consider the nonlinear dynamic behaviour of soils and simultaneously investigate seismic wave propagation in complex 2D/3D geological structures! Recent advances in numerical formulations and constitutive models in such complex situations are presented and discussed in this paper. A crucial issue is the availability of the field/laboratory data to feed and validate such models.

Modeling seismic wave propagation

Many various numerical methods are available to model seismic wave propagation and we will discuss them first. Afterwards, the issue of seismic wave amplification (site effects) in both linear (weak motion) and nonlinear (strong motion) ranges will then be examined. To analyze seismic wave propagation in 2D or 3D geological structures, various numerical methods are available (Fig. 1):

 the finite difference method is accurate in elastodynamics but is mainly adapted to simple geometries [START_REF] Bohlen | Accuracy of heterogeneous staggered-grid finite-difference modeling of rayleigh waves[END_REF][START_REF] Frankel | A Three-Dimensional Simulation of Seismic Waves in the Santa Clara Valley, California, from a Loma Prieta Aftershock[END_REF][START_REF] Moczo | 3D Heterogeneous Staggered-Grid Finite-Difference Modeling of Seismic Motion with Volume Harmonic and Arithmetic Averaging of Elastic Moduli and Densities[END_REF][START_REF] Virieux | P-SV Wave Propagation in Heterogeneous Media: Velocity-Stress Finite-Difference Method[END_REF]),  the finite element method is efficient to deal with complex geometries and numerous heterogeneities (even for inelastic constitutive models [START_REF] Bonilla | Computation of Linear and Nonlinear Site Response for Near Field Ground Motion[END_REF]) but has several drawbacks such as numerical dispersion and numerical damping [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF][START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF][START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF], Semblat, 2000a[START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF] and (consequently) numerical cost in 3D elastodynamics,  the spectral element method has been increasingly considered to analyse 2D/3D wave propagation in linear media with a good accuracy due to its spectral convergence properties [START_REF] Chaljub | Spectral Element Analysis in Seismology[END_REF][START_REF] Faccioli | Spectral Domain Decomposition Methods for the Solution of Acoustic and Elastic Wave Equations[END_REF][START_REF] Komatitsch | The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures[END_REF],  the boundary element method allows a very good description of the radiation conditions but is preferably dedicated to weak heterogeneities and linear constitutive models [START_REF] Beskos | Boundary elements methods in dynamic analysis: Part II (1986-1996[END_REF][START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Dangla | A plane strain soil-structure interaction model[END_REF], 2005[START_REF] Sanchez-Sesma | Seismic Response of Three-Dimensional Alluvial Valleys for Incident P, S and Rayleigh Waves[END_REF], Semblat 2008, 2000b). Recent developments have been proposed to reduce the computational cost of the method especially in the high frequency range [START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF], 2009[START_REF] Fujiwara | The fast multipole method for solving integral equations of three-dimensional topography and basin problems[END_REF],  the finite volume method was recently developed in the field of elastodynamics [START_REF] Glinsky-Olivier | A finite-volume method for the 2D seismic wave propagation[END_REF],  the Aki-Larner method which takes advantage of the frequency-wavenumber decomposition but is limited to simple geometries [START_REF] Aki | Surface Motion of a Layered Medium Having an Irregular Interface due to Incident Plane SH Waves[END_REF][START_REF] Bouchon | A Boundary Integral Equation-Discrete Wavenumber Representation Method to Study Wave Propagation in Multilayered Media Having Irregular Interfaces[END_REF]),  the scaled boundary finite element method which is a kind of solution-less boundary element method [START_REF] Wolf | The scaled boundary finite element method[END_REF],  other methods such as series expansions of wave functions [START_REF] Liao | A series solution and numerical technique for wave diffraction by a three-dimensional canyon[END_REF][START_REF] Sanchez-Sesma | Diffraction of elastic waves by three-dimensional surface irregularities[END_REF]. Furthermore, when dealing with wave propagation in unbounded domains, many of these numerical methods raise the need for absorbing boundary conditions to avoid spurious reflections. Since each method has specific advantages and shortcomings (Table I), it is consequently often more interesting to combine two methods to take advantage of their peculiarities. It is for instance possible to couple FEM and BEM [START_REF] Aochi | A complete BIEM-FDM-FEM simulation of an earthquake scenario[END_REF][START_REF] Dangla | A plane strain soil-structure interaction model[END_REF][START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF]) allowing an accurate description of the near field (FEM model including complex geometries, numerous heterogeneities and nonlinear constitutive laws) and a reliable estimation of the far-field (BEM involving accurate radiation conditions). The main advantage of the method in elastodynamics is the accurate description of the infinite extension of the medium through exact radiation conditions [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF], Semblat 2008). In contrast to other discretization methods, there is no need to consider absorbing boundary conditions. The Boundary Element Method involves singular integrals [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF]. There are three main kinds of singularities: (i) weak singularity (i.e. in the ordinary Riemann sense), (ii) strong singularity (i.e. in the Cauchy principal value sense) or (iii) hyper-singularity (i.e. in the Hadamard finite part sense). Strong and hyper-singular integrals have to be converted to regular ones in the regularization of the BEM formulations [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Sladek | Singular integrals and boundary elements[END_REF]). Indeed, weak singularity is not treated by regularization. However, from the point of view of numerical integrations, one should devote a great attention to the evaluation of these integrals because standard integration quadratures fail in accuracy [START_REF] Beskos | Boundary elements methods in dynamic analysis: Part II (1986-1996[END_REF][START_REF] Dangla | A plane strain soil-structure interaction model[END_REF]. Therefore each type of singularity has to be treated by appropriate techniques [START_REF] Bonnet | Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods[END_REF][START_REF] Dangla | A simple and efficient regularization method for 3D BEM: application to frequency-domain elastodynamics[END_REF][START_REF] Sladek | Singular integrals and boundary elements[END_REF]). It is noticed that the regularization can be performed either before or after the discretization, i.e. in the global or local (intrinsic) coordinate space, as observed in some papers mentioned above. A comprehensive review of BEM in dynamic analysis has been proposed by [START_REF] Beskos | Boundary elements methods in dynamic analysis: Part II (1986-1996[END_REF]. An example (pressurized cavity in a full space) showing the influence of the regularization process on the numerical solution is proposed in Fig. 2: the uncorrected solution is very far from the regularized one which matches well the analytical solution [START_REF] Dangla | A simple and efficient regularization method for 3D BEM: application to frequency-domain elastodynamics[END_REF].

The main drawback of the BEM is that it leads to full non symmetric matrix systems. For large models, the numerical cost and memory storage requirements may then by huge. Recent researches considered symmetric Galerkin boundary element methods more efficient for large numbers of unknowns [START_REF] Bonnet | Symmetric Galerkin boundary element methods[END_REF]. Current researches try to reduced the number of computations through Fast Multipole formulations initially developed in the field of physics [START_REF] Fujiwara | The fast multipole method for solving integral equations of three-dimensional topography and basin problems[END_REF][START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF], 2009). Such alternative approaches are discussed in the following. [START_REF] Dangla | A simple and efficient regularization method for 3D BEM: application to frequency-domain elastodynamics[END_REF].

Fast Multipole formulation for the Boundary Element Method

In the field of physics (Maxwell or Laplace equations), recent advances in boundary element methods lead to a very important decrease of the computational cost. Instead of point to point interaction as in classical boundary element methods, the fast multipole method [START_REF] Greengard | Accelerating fast multipole methods for the Helmholtz equation at low frequencies[END_REF][START_REF] Fujiwara | The fast multipole method for solving integral equations of three-dimensional topography and basin problems[END_REF] considers interactions between groups of points (cells centered on a multipole, Fig. 3) and hence avoids multiple computations of nearly identical terms corresponding to very close points. It is possible to apply this method at a single scale or even at various scales through a multilevel approach. The size of each cell around its multipole depends on the distance to the other cells: the larger the distance between cells, the larger the cells (Fig. 3). The computations of the singular integrals are then performed through this approach. They are split in an integral on the surface around the singularity and another integral on the complementary distant surface. The first one is estimated using classical regularization techniques, whereas the latter is computed with a fast multipole algorithm [START_REF] Greengard | Accelerating fast multipole methods for the Helmholtz equation at low frequencies[END_REF]. The advantage for the computational cost is very significant since it depends on N 2 for classical methods and on N or NlogN for the fast multipole method [START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF]. Furthermore, the computational cost is reduced for both the memory storage and the calculation time. This method allows the analysis of very large problems involving millions of unknowns on a singleprocessor PC (to compare to several tens of thousands previously). Current researches also investigate a fast multipole method well-adapted to the Helmholtz equation and even to elastodynamics [START_REF] Fujiwara | The fast multipole method for solving integral equations of three-dimensional topography and basin problems[END_REF][START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF], 2009). The fast multipole method then allows the computation of very large models considering a larger number of heterogeneities, a more realistic geometrical representation of geological structures (especially in 3D) as well as higher frequency values (detailed modeling of short wavelengthes amplification).

Modeling seismic wave amplification by the Boundary Element Method

Amplification of seismic waves in geological structures

The amplification of the seismic motion mainly occurs in alluvial deposits the characteristics of which (geometry, wave velocities) control the amplification process [START_REF] Bard | The two dimensional resonance of sediment filled valleys[END_REF], Dobry 1976, Duval 1998[START_REF] Sanchez-Sesma | Diffraction of elastic waves by three-dimensional surface irregularities[END_REF], 2000, Semblat 2000, 2003). Considering the Boundary Element Method, seismic wave amplification (i.e. site effects) can be analyzed numerically taking into account the infinite extension of the medium as well as complex geometrical features. The main phenomena that can be recovered by the BEM models are the following:

 amplification due to the (1D) velocity contrast between soil layers,  focusing effects due to complex layers or free-surface geometries,  basins effects (2D/3D) due the trapped surfaces waves in surficial layers. Significant wave scattering may also be found around strong topographic irregularities such as crests or hills slopes [START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF][START_REF] Paolucci | Amplification of earthquake ground motion by steep topographic irregularities[END_REF][START_REF] Reinoso | Three-dimensional scattering of seismic waves from topographical structures[END_REF][START_REF] Sanchez-Sesma | Diffraction of elastic waves by three-dimensional surface irregularities[END_REF], Semblat 2009). Such phenomena are called topographic site effects and may aso be studied by the Boundary Element Method. In both case (alluvial basins and topographic irregularities), BEM models in 2D or 3D allow the estimation of the amplification factor of the seismic motion. The BEM formulation can be considered in frequency-domain for long durations and narrow bandwidths [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Dangla | A simple and efficient regularization method for 3D BEM: application to frequency-domain elastodynamics[END_REF] or time-domain for shorter durations and larger bandwidths [START_REF] Gaul | A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains[END_REF][START_REF] Jin | A 2D time-domain boundary element method with damping[END_REF][START_REF] Manolis | Boundary Element Methods in Elastodynamics[END_REF]). Other numerical methods may also be considered as discussed in previous sections [START_REF] Faccioli | Spectral Domain Decomposition Methods for the Solution of Acoustic and Elastic Wave Equations[END_REF][START_REF] Glinsky-Olivier | A finite-volume method for the 2D seismic wave propagation[END_REF][START_REF] Komatitsch | The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures[END_REF][START_REF] Moczo | 3D Heterogeneous Staggered-Grid Finite-Difference Modeling of Seismic Motion with Volume Harmonic and Arithmetic Averaging of Elastic Moduli and Densities[END_REF][START_REF] Virieux | P-SV Wave Propagation in Heterogeneous Media: Velocity-Stress Finite-Difference Method[END_REF]).

Topographic amplification of seismic waves

The analysis of wave scattering around topographic irregularities has been performed by various authors for canonical configurations [START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF][START_REF] Reinoso | Three-dimensional scattering of seismic waves from topographical structures[END_REF][START_REF] Sanchez-Sesma | Diffraction of elastic waves by three-dimensional surface irregularities[END_REF] or for actual 2D/3D topographies [START_REF] Paolucci | Amplification of earthquake ground motion by steep topographic irregularities[END_REF][START_REF] Semblat | Seismic Site Effects in a Deep Alluvial Basin: Numerical Analysis by the boundary element method[END_REF]. As shown in Fig. 4 in the case of Caracas [START_REF] Semblat | Seismic Site Effects in a Deep Alluvial Basin: Numerical Analysis by the boundary element method[END_REF], the maximum amplification Amax along a complex topography may lead to significant changes in the seismic wavefield. When comparing the top plot at frequency 0.25 Hz to the bottom one (0.35 Hz), these topographic site effects are nevertheless found to be less than those found in alluvial basins. Detailed analyses of such effects are proposed in [START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF][START_REF] Paolucci | Amplification of earthquake ground motion by steep topographic irregularities[END_REF], Semblat 2009). [START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF] proposed empirical laws to estimate the topographic effects for 2D steep crests. For 3D models, efficient numerical methods are needed. [START_REF] Paolucci | Amplification of earthquake ground motion by steep topographic irregularities[END_REF] studied an actual steep italian site using the Spectral Element Method. Since the development of Fast Multipole BEMs in elastodynamics allows to handle large scale 3D models, several large canonical 3D topographies have been recently studied [START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF]. As shown in Figure 5 for an ellipsoidal canyon and an oblique incident plane P-wave [START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF], the FM-BEM solution is in good agreement with previous solutions obtained by [START_REF] Reinoso | Three-dimensional scattering of seismic waves from topographical structures[END_REF] for both the vertical and horizontal motion components. Several other canonical configurations were defined in the framework of the QSHA research project (Quantitative Seismic Hazard Assessment) and are used to compare the efficiency and accuracy of various numerical methods (data available at http://qsha.unice.fr/). Current researches aim at investigating actual 3D topographies by the Fast Multipole Method. example in Caracas, Venezuela [START_REF] Semblat | Seismic Site Effects in a Deep Alluvial Basin: Numerical Analysis by the boundary element method[END_REF]. (data for various canonical configurations available at http://qsha.unice.fr/)

Amplification of seismic waves in alluvial basins

The amplification of seismic waves in alluvial basins is generally larger than that observed/computed around topographic irregularities [START_REF] Baise | The effect of shallow San Francisco bay sediments on waveforms recorded during the Mw 4.6 Bolinas California Earthquake[END_REF][START_REF] Beauval | Quantification of frequency-dependent lengthening of seismic groundmotion duration due to local geology: applications to the Volvi area (Greece)[END_REF][START_REF] Frankel | A Three-Dimensional Simulation of Seismic Waves in the Santa Clara Valley, California, from a Loma Prieta Aftershock[END_REF][START_REF] Sanchez-Sesma | Seismic Response of Three-Dimensional Alluvial Valleys for Incident P, S and Rayleigh Waves[END_REF], Semblat 2009[START_REF] Sommerville | Emerging art: earthquake ground motion[END_REF]). Since extensive field measurements were performed on Volvi EuroSeisTest, Greece [START_REF] Bard | Wave propagation in complex geological structures and their effects on strong ground motion[END_REF][START_REF] Beauval | Quantification of frequency-dependent lengthening of seismic groundmotion duration due to local geology: applications to the Volvi area (Greece)[END_REF], Chavez-Garcia 2000[START_REF] Pitilakis | Site effects : recent considerations and design provisions[END_REF]), BEM models were performed for this site (Semblat 2005). The numerical results were compared to numerous seismological field data from various earthquake recordings.

The Volvi basin is 6km wide and 250m deep. As shown in Figure 6, the amplification level changes from one frequency to the other. The influence on the amplification process of both the basin geometry and the soil layering were studied in (Semblat 2005). Since the geometry of the basin is complex, strong focusing effects occur leading to large amplification levels (around 10) at some peculiar frequencies (Figure 6). However, depending on the accuracy of the geological model, the description of the soil layering has a strong influence at higher frequencies (Figure 6, bottom right). It has been quantified for a simplified (2 layers) as well as detailed (6 layers) model of the Volvi profile in (Semblat 2005). The detailed analysis leads to larger amplitude trapped surface waves at higher frequencies.

F1=0.6Hz, A1=2.9 F2=0.8Hz, A2=9.5 F4=1.2Hz, A4=7.5 F6=2.4Hz, A6=7.3 Figure 6. Seismic wave amplification in the Volvi-EuroSeisTest basin, Greece (Semblat 2005): BEM simulations at different frequencies Fi and maximum related amplifications.

3 Modeling seismic wave propagation by the Finite/Spectral Element Method

Numerical wave dispersion

Considering the finite/spectral element method to model seismic wave propagation, as for other numerical methods, two types of numerical errors can be studied [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF][START_REF] Chaljub | Spectral Element Analysis in Seismology[END_REF][START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF] : the relative period error and the algorithmic damping. The relative period error is different from one integration scheme to another [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF]. For wave propagation problems, the relative period error appears in the estimation of wave-velocities and is called numerical dispersion. Seismic wave propagation in a numerical scheme then depends on the element size, the integration scheme, the element type, etc. This phenomenon is called numerical dispersion as a reference to physical dispersion leading to a dependence of the wave-velocity on frequency [START_REF] Deraemaeker | Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[END_REF][START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF][START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF], Semblat 2000b). To model wave propagation phenomena by the finite element method or the finite difference method, one has to control the accuracy of the numerical scheme since the numerical error tends to increase during the propagation process.

Various theoretical works concern the analysis of the numerical error made in the estimation of the approximated wave number [START_REF] Deraemaeker | Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[END_REF][START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF], Semblat 2009). For instance, it was shown that there is a cutoff frequency above which there is no propagation at all [START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF]. Depending on the frequency of the excitation, the numerical wave propagates slower or faster than the theoretical solution. It is then necessary to analyze numerical dispersion of seismic waves and precisely quantify the numerical error.

Higher order finite elements vs spectral elements

Numerical wave dispersion is also influenced by the order of the shape functions of the finite elements. Higher order finite elements are for instance known to have a very good precision for elastic-plastic computations. In acoustics and elastodynamics, several theoretical works propose analytical expressions to estimate the numerical dispersion [START_REF] Deraemaeker | Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[END_REF][START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF]. Some comparisons between lower order and higher order finite elements to model wave propagation were recently proposed by [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF]. In Figure 7, the results from 1D FEM simulations (Semblat 2000b) are plotted for various shape functions orders.

The same number of degrees of freedom is considered in each case. The linear finite elements (Figure 7, top) are shown to have strong numerical dispersion (numerical wave velocity artificially increased). For quadratic finite elements (Figure 7, middle), the accuracy is satisfactory when compared to the theoretical delays (vertical dashed lines). Finally, for the higher order finite elements (Figure 7, bottom), the numerical wave dispersion is very low. These results show that, for the same number of DOFs, the accuracy of higher order finite element for wave propagation simulations is much better. In addition to "classical" higher order finite elements, spectral finite elements have increasingly been studied since their spectral convergence is of great interest for wave propagation simulations [START_REF] Chaljub | Spectral Element Analysis in Seismology[END_REF][START_REF] Faccioli | Spectral Domain Decomposition Methods for the Solution of Acoustic and Elastic Wave Equations[END_REF][START_REF] Komatitsch | The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures[END_REF]. The spectral elements are generally considered for high order only (4 to 8) since the difference with classical finite elements is not significant below. An evaluation of the cost effectiveness (and accuracy) of both Chebyshev spectral and p-finite elements has been proposed by [START_REF] Dauksher | An evaluation of the cost effectiveness of Chebyshev spectral and p-finite element solutions to the scalar wave equation[END_REF].

Modeling linear attenuation

The analysis of seismic wave propagation in attenuating media, such as soils and rocks, raises the need for suitable attenuation models. One of the most common damping model in structural dynamics is the Rayleigh formulation giving a damping matrix as a linear combination of stiffness and mass matrices [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF]. The damping ratio is then minimum for a given frequency and infinite for zero and infinite frequencies (Semblat 2009). In the field of structural dynamics, the main advantage of the Rayleigh formulation is that it leads to a diagonal damping marix in the real eigenmodes base [START_REF] Chopra | Dynamics of Structures, Theory and Applications to Earthquake Eng[END_REF].

For wave propagation phenomena, the Rayleigh damping formulation is equivalent to a particular rheological model: generalized Maxwell model [START_REF] Semblat | Rheological interpretation of Rayleigh damping[END_REF]. For weak to moderate damping values, this equivalence allows an easier and explicit estimation of the two coefficients of the Rayleigh formulation from experimental behaviour parameters. Other types of "matrix based" damping formulations are available such as the Caughey formulation [START_REF] Chopra | Dynamics of Structures, Theory and Applications to Earthquake Eng[END_REF], Semblat 2009) corresponding to a generalization of the Rayleigh formulation and allowing complex damping-frequency dependences.

Various types of mechanical or rheological viscoelastic formulations were also proposed. [START_REF] Kjartansson | Constant Q wave propagation and attenuation[END_REF] investigated frequency constant attenuation ("CQ") models for wave propagation problems. [START_REF] Day | Numerical simulation of wavefields using a Padé approximant method[END_REF] considered a Padé approximant to simulate wave propagation in attenuatingmedia. [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] proposed a rheological model leading to a nearly constant attenuation ("NCQ") in a given frequency range. [START_REF] Carcione | Time-domain Modeling of Constant-Q Seismic Waves Using Fractional Derivatives[END_REF] studied a damping formulation based on fractional derivatives and [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF] discussed time domain simulations using rheological models. All these approaches are formulated in the linear range and need, for time domain computations, to deal with memory variables which may be costly. As discussed in the following, nonlinear constitutive laws are often necessary in the case of strong earthquakes.

Radiation conditions at infinity

For finite differences, finite or spectral element methods, it is necessary to take into account the radiation conditions through efficient techniques to avoid spurious wave reflections at the mesh boundaries. It can be done considering absorbing boundary conditions or infinite elements for purely solid [START_REF] Chadwick | Diffraction of short waves modelled using new mapped wave envelope finite and infinite elements[END_REF] as well as multiphase media [START_REF] Modaressi | Paraxial Approximation for Poroelastic Media[END_REF]. It is nevertheless difficult for very heterogeneous media [START_REF] Chammas | Characterization of heterogeneous soils using surface waves : homogenization and numerical modeling[END_REF].

Absorbing layer methods rather try to model the wave attenuation in a thick layer at the medium boundaries (Semblat 2009). They have been increasingly developed in recent years since their efficiency appears good in various configurations. They are generally known as "Perfectly Matched Layer" methods. Various types of PML formulations have recently been proposed:  Classical PML [START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF]: the wavefield propagating in the layer is damped according to a onedimensional amplitude decrease law similar to that of an attenuating medium:

  i x x x ) ( ~  (1)
This complex valued law leads to an amplitude decrease of the numerical wave along the x-axis. As illustrated in Figure 8, grazing incidences are related to small horizontal wavenumbers and a limited amplitude decrease will be obtained. This uniaxial law is thus not sufficient to have an optimal absorbing effect in all configurations.

 Filtering PML [START_REF] Festa | PML absorbing boundaries[END_REF][START_REF] Festa | The Newmark scheme as a Velocity-Stress Time staggering: An efficient PML for Spectral Element simulations of elastodynamics[END_REF]: for shallow numerical models, the classical PML formulation may also amplify surface waves (Figure 8). A filtering PML formulation (F-PML) has been proposed by [START_REF] Festa | The Newmark scheme as a Velocity-Stress Time staggering: An efficient PML for Spectral Element simulations of elastodynamics[END_REF] to avoid this problem: [START_REF] Meza-Fajardo | A Nonconvolutional, Split-Field, Perfectly Matched Layer for Wave Propagation in Isotropic and Anisotropic Elastic Media: Stability Analysis[END_REF]: another alternative to ensure the numerical stability of PMLs (e.g. grazing incidences) is the multidirectional PML formulation (M-PML) recently proposed by Meza Fajardo and Papageorgiou [START_REF] Meza-Fajardo | A Nonconvolutional, Split-Field, Perfectly Matched Layer for Wave Propagation in Isotropic and Anisotropic Elastic Media: Stability Analysis[END_REF]. The original uniaxial law is generalized in the multi-axial case through the following system:

c i x x x       ) ( ~ (2)  Multidirectional PML
               z i z z y i y y x i x x x z x y x x ) ( ) ( ) ( 1 ~1 ~1 ~      (3)
It is thus possible to choose the attenuation vector  in order to optimize the absorbing effect depending on the wave type and incidence. Another interest of this formulation is its numerical stability for strongly anisotropic media (Figure 8).

Modeling strong ground motion and non linear effects

To model strong ground motion, it is necessary to combine seismological tools to geotechnical methods involving nonlinear constitutive laws [START_REF] Heuze | Estimating site-specific strong earthquake motion[END_REF]. The nonlinear soil models should reproduce one of the main features of soil dynamic/cyclic behaviour: the shear modulus decreases and the hysteretic damping increases with strain amplitude [START_REF] Seed | Moduli and damping factors for dynamic analyses of cohesionless soils[END_REF]). Such models are generally considered in the framework of the finite element method and often involve Masing type cyclic behaviour [START_REF] Bonilla | Computation of Linear and Nonlinear Site Response for Near Field Ground Motion[END_REF]. The influence of plasticity, confining pressure and pore-pressure build up can also be analyzed with more complex models [START_REF] Hujeux | A double memory model with multiple mechanisms for cyclic soil behaviour[END_REF][START_REF] Gyebi | Finite element analysis of viscoplastic soils with Q-factor[END_REF][START_REF] Lade | Elasto-plastic stress-strain theory for cohesionless soil with curved yield surfaces[END_REF][START_REF] Loret | Growth and decay of acceleration waves in non-associative elastic-plastic fluidstaturated porous media[END_REF][START_REF] Mellal | A simplified numerical approach for nonlinear dynamic analysis of multilayered media[END_REF][START_REF] Park | Soil damping formulation in nonlinear time domain site response analysis[END_REF][START_REF] Prevost | A simple plasticity theory for frictional cohesionless soils[END_REF]. It is then possible to have a deeper insight into the liquefaction processes. More and more strong motion data from well-controlled bore-holes are available and allow comparisons with nonlinear computations.

In the last decades, linear (viscoelastic) equivalent models were also extensively used to have a simplified description (i.e. few parameters) of the shear modulus decrease and the hysteretic damping increase [START_REF] Schnabel | SHAKE-A computer program for equation response analysis of horizontally layered sites[END_REF]. Simplified models are interesting since they may allow combined seismological/geotechnical computations to simultaneously take into account basin effects at large scales (2D/3D) and nonlinear local effects. Recent researches proposed new simplified models avoiding some drawbacks of equivalent models such as frequency independence [START_REF] Assimaki | Model for dynamic shear modulus and damping in granular soils[END_REF][START_REF] Kausel | Seismic simulation of inelastic soils via frequency-dependent moduli and damping[END_REF]). As shown in Fig. 9 (left), this nonlinear frequency dependent model combines a strain spectrum with the classical G() and () curves. [START_REF] Delépine | A simplified non linear model to analyze site effects in alluvial deposits[END_REF][START_REF] Delépine | Nonlinear viscoelastic wave propagation: an extension of Nearly Constant Attenuation (NCQ) models[END_REF] proposed another alternative: a simple constitutive model generalizing the classical NCQ (Nearly Constant) attenuation model [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] in the nonlinear range. This extended NCQ model leads to decreasing shear modulus and increasing damping/attenuation. This model allows the analysis of seismic waves propagation from strong earthquakes and may be used at large scales (e.g. alluvial basins). Finally, a crucial issue is the availability of the field/laboratory data to feed and validate such models.

Conclusion

This paper presents recent results in various research fields related to seismic wave propagation in geological structures. New methods or major improvements of existing methods (staggered grid FD, Finite Volumes, Fast Multipole-BEM) are promising ways to better model seismic wave propagation, amplification and interactions in 3D surficial geological structures. The major advances of current works also concern improved treatments of the radiation conditions for spectral elements/finite elements (absorbing layer methods such as F-PML or M-PML).

Another key point is strong seismic motion: recent simplified nonlinear constitutive models allow the analysis of seismic wave propagation from strong earthquakes at large scales. Nevertheless, several topics have not been discussed in this paper: fault mechanics and the influence of fluids [START_REF] Aochi | A complete BIEM-FDM-FEM simulation of an earthquake scenario[END_REF], soil-structure interaction [START_REF] Kham | Seismic Site-City Interaction: Main Governing Phenomena through Simplified Numerical Models[END_REF][START_REF] Sextos | Inelastic dynamic analysis of RC bridges accounting for spatial variability of ground motion, site effects and soil-structure interaction phenomena. Part 1: Methodology and analytical tools[END_REF]) and site-city interaction [START_REF] Bielak | Earthquake ground motion and structural response in alluvial valleys[END_REF][START_REF] Clouteau | Modifications of the ground motion in dense urban areas[END_REF][START_REF] Groby | Simulation of seismic response in a city-like environment[END_REF], Semblat 2008[START_REF] Tsogka | Seismic response of a set of blocks partially embedded in soft soil[END_REF][START_REF] Wirgin | Effects of buildings on the duration and amplitude of ground motion in Mexico city[END_REF], soil-pile interaction and liquefaction [START_REF] Brennan | Effectiveness of vertical drains in mitigation of liquefaction[END_REF][START_REF] Takahashi | Liquefaction-induced large displacement of pile-supported wharf[END_REF], etc.

Acknowledgements

The author would like to thank all his colleagues and students who contributed to some of the results presented herein and the related papers listed in the References.

Figure 1 .

 1 Figure 1. Various numerical methods to model seismic wave propagation: (a) the Finite Difference Method, (b) the Finite Element Method, (c), the Spectral Element Method, (d) the Boundary Element Method, (e) the Discrete Element Method, (f) the Finite Volume Method.

Figure 3 .

 3 Figure 3. Comparison of the principles of the standard BEM and the Fast Multipole BEM.

Figure 4 .

 4 Figure 4. Topographic (top) vs statigraphic (bottom) seismic wave amplification:example in Caracas, Venezuela[START_REF] Semblat | Seismic Site Effects in a Deep Alluvial Basin: Numerical Analysis by the boundary element method[END_REF].

Figure 5 .

 5 Figure5. Validation of the Fast Multipole BEM from a 3D configuration studied by[START_REF] Reinoso | Three-dimensional scattering of seismic waves from topographical structures[END_REF]:diffraction of an oblique plane P-wave by a semi-ellipsoidal canyon[START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF].(data for various canonical configurations available at http://qsha.unice.fr/)

Figure 7 .

 7 Figure 7. Numerical wave dispersion for various finite element orders: linear (top), quadratic (middle) and fourth-order (bottom) (Semblat 2000b).

Figure 8 .

 8 Figure 8. Various critical configurations for the efficiency/stability of PML absorbing conditions (from Semblat, 2009).

aFigure 9 .

 9 Figure 9. Recent improvements of the equivalent linear model for the analysis of strong seismic motion in soils: (a) Frequency dependent nonlinear model (Kausel 2002); (b) Nonlinear "extended NCQ" model (Delépine 2009).

Table I .

 I Features of various numerical methods for wave propagation modeling.

	Method	FDM	FEM	SEM	BEM	DEM	FVM
	Continuous	yes	yes	yes	yes	no(/yes)	yes
	Constitutive	lin./NL	lin./NL	lin./NL	linear	lin./NL	lin.
	Radiation	abs. bound.	abs. bound.	abs. bound.	intrinsic	abs. bound.	abs. bound.

2 Modeling seismic wave propagation by the Boundary Element Method

2.1 Interest of the Boundary Element Method ("BEM")


Shiyan Hu

Pierrick Coupé

Jens C Pruessner

D Louis Collins

Validation of appearance-model based segmentation with patch-based refinement on medial temporal lobe structures

This paper presents a new automatic segmentation scheme that combines active appearance (AAM) modeling and patch-based label fusion into a segmentation framework. AAM, which uses eigen-decomposition to analyze the statistical variation of image intensity and shape information over the population, is used to capture the global shape characteristics of the structure of interest with a generative model, while patch-based label fusion, which uses a non-local means method to compare the image local intensity properties, is applied to locally refine the segmentation results along the structure boundary area to improve the segmentation accuracy. The proposed segmentation scheme is used to segment human medial temporal lobe structures, which have low intensity contrast in MRI and complexity in shape. The experiments demonstrate that this new segmentation scheme is computationally efficient and robust. In a leave-one out validation with fifty-four normal young subjects, the method yields a mean Dice κ of 0.87 for hippocampus, 0.81 for amygdala, 0.73 for para-hippocampal complex cortex, and 0.73 for perirhinal cortex between manual and automatic labels.

Introduction

The medial temporal lobe, located in the inner side of each temporal lobe, is the most important part of the limbic system, and includes the hippocampus (HC), the amygdala (AG), and surrounding cortical areas (entorhinal (EC), perirhinal (PRC), and parahippocampal (PHC) cortices). These structures play important roles in learning, memory, and emotion. HC is the the most critical component of the medial temporal lobe memory system, and plays an essential role in memory and spatial navigation. AG is strongly involved in emotional and social processing, particularly fear and worry. EC is the main interface between HC and neocortex, and plays an important role in the formation and optimization of spatial memories. PRC is involved in both visual perception and memory, and PHC is involved in scene recognition and social context [START_REF] Barense | Functional specialization in the human medial temporal lobe[END_REF] [START_REF] Baxter | Involvement of medial temporal lobe structures in memory and perception[END_REF]. Recently, medial temporal lobe structures have received considerable attention due to their importance in neurological diseases and disorders [START_REF] Mori | Medial temporal structures relate to memory impairment in Alzheimers disease: an MRI volumetric study[END_REF]. For example, many researchers [START_REF] Wang | Hippocampal shape abnormalities in early AD: A replication study[END_REF] [START_REF] Duzel | Hippocampal atrophy in temporal lobe epilepsy is correlated with limbic systems atrophy[END_REF] showed that HC volume changes were an important marker of the early stage of Alzheimer's disease and temporal lobe epilepsy. Because of the importance of these structures in the neurodegeneration, there is significant interest in developing accurate, robust, and reliable segmentation techniques to automatically extract those structures from magnetic resonance (MR) imaging for volume and shape analyses.

Manual segmentation is considered highly accurate, and treated as the gold standard so far. However, manual segmentation is time consuming, requires anatomical expertise, and may have important intra and inter-rater variability. To overcome the disadvantages of manual segmentation, many automatic segmentation techniques have been proposed. Most model-based segmentation techniques can be grouped into the following three categories: deformable models [START_REF] Chupin | Anatomically constrained region deformation for the automated segmentation of the hippocampus and the amygdala: method and validation on controls and patients with Alzheimers disease[END_REF], appearance-based models [START_REF] Klemencic | Non-rigid registration based active appearance models for 3D medical image segmentation[END_REF] [8], and atlas-based techniques [START_REF] Fischl | Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain[END_REF] [START_REF] Collins | Automatic 3D modelbased neuro-anatomical segmentation[END_REF].

To avoid the bias in the segmentation by using a single template, Heckemann et al. [START_REF] Heckemann | Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[END_REF], Shattuck et al. [START_REF] Shattuck | Construction of a 3D probabilistic atlas of human cortical structures[END_REF] and Aljabar et al. [START_REF] Aljabar | Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[END_REF] proposed multi-atlas based methods to improve the segmentation efficiency. Collins and Pruessner [START_REF] Collins | Towards accurate, automatic segmentation of the hippocampus and amygdala from MRI by augmenting animal with a template library and label[END_REF] incorporated the label fusion into the multi-atlas warping, and Wang et al. in [START_REF] Wang | The Alzheimer's Disease Neuroimaging Initiative: A learningbased wrapper method to correct systematic errors in automatic image segmentation: Consistently improved performance in hippocampus, cortex and brain segmentation[END_REF] used the multi-atlas technique with error correction to yield the best published results for HC segmentation. Coupé et al. [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF] used non-local means patch-based approach to weight the expert manual segmentation in a library of templates (see Table 1 for the segmentation method review). However, neither atlas-based or patch-based methods explicitly incorporate the structure's shape information into the segmentation. [START_REF] Powell | Registration and machine learning based automated segmentation of subcortical and cerebellar brain structures[END_REF] Machine learning 0.85 Lijn 2008 [START_REF] Lijn | Hippocampus segmentation in MR images using atlas registration, voxel classification, and graph cuts[END_REF] Multi-atlas+graph cuts 0.86 Morey 2009 [START_REF] Morey | A comparison of automated segmentation and manual tracing for quantifying hippocampal and amygdala volumes[END_REF] FST/FIRST 0.79 0.73 Aljabar 2009 [START_REF] Aljabar | Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy[END_REF] Multi-atlas+template selection 0.84 0.78 Lotjonen 2010 [START_REF] Lotjonen | The Alzheimer's Disease Neuroimaging Initiative: Fast and robust multi-atlas segmentation of brain magnetic resonance images[END_REF] Multi-atlas+intensity modeling 0.82∼0.88 0.77 Collins 2010 [START_REF] Collins | Towards accurate, automatic segmentation of the hippocampus and amygdala from MRI by augmenting animal with a template library and label[END_REF] Multi-atlas+label fusion 0.89 0.83 Benavides 2010 [START_REF] Benavides | Manual validation of freesurfers automated hippocampal segmentation in normal aging, mild cognitive impairment, and alzheimer disease subjects[END_REF] FreeSurfer 0.78 Sabuncu 2010 [START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF] Label fusion 0.82∼0.87 0.8∼0.82 Coupé 2011 [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF] Non-local means 0.88 Patenaude 2011 [START_REF] Patenaude | A bayesian model of shape and appearance for subcortical brain segmentation[END_REF] Bayesian appearance 0.81 0.74 Wang 2011 [START_REF] Wang | The Alzheimer's Disease Neuroimaging Initiative: A learningbased wrapper method to correct systematic errors in automatic image segmentation: Consistently improved performance in hippocampus, cortex and brain segmentation[END_REF] Multi-atlas+error correction 0.89∼0.91 Bishop 2011 [START_REF] Bishop | Novel Fast Marching for Automated Segmentation of the Hippocampus (FMASH): Method and validation on clinical data[END_REF] FMASH 0.8∼0.82 Khan 2011 [START_REF] Khan | Optimal Weights for Local Multi-atlas Fusion using Supervised Learning and Dynamic Information (SuperDyn): Validation on Hippocampus Segmentation[END_REF] Multi-atlas+spatially-local selection 0.83∼0.85 *parahippocampal + ambient gyri.

To integrate the shape constraint into the segmentation, we combine the appearance model and patch-based technique into a general segmentation frame-work. The main contribution of this paper includes: 1) the proposed method combines the appearance model and non-local means patch-based methods where appearance models are used to capture the global shape variation, and the nonlocal means method is used to locally refine the segmentation, 2) the proposed method is applied to segment all medial temporal lobe structures. Compared to HC and AG, the other medial temporal lobe structures, like PRC and PHC, have much greater anatomical variability, and 3) nonlinear registration instead of linear registration is applied to both training and testing data to improve the alignment between subjects.

Method

Appearance-model based Global Segmentation

Appearance-model based global segmentation applies the eigen-decomposition technique on gray MR images and shape data to capture the statistical variations of the gray intensity and shape information of the training data. Based on the eigenvectors derived from the training data, using the notation from [START_REF] Hu | Joint level set shape modeling and appearance modeling for brain structure segmentation[END_REF] the final shape and gray image can be given by,

φ = φ + P φ Q s w-1 s c g t1 = ḡt1 + P g,t1 Q g,t1 c (1)
The segmentation is achieved by minimizing the difference between the test image and the one synthesized from Eq. 1. The cost function in the least square measure can be written as:

E = Np j=1 (I t1,j -g t1,j ) 2 (2)
where I t1,j and g t1,j are the intensity of the j-th voxel of the T1 test MR image and the synthesized T1 image respectively. N p is the total number of voxels in each image. The detailed appearance-based segmentation method was described in [START_REF] Hu | Joint level set shape modeling and appearance modeling for brain structure segmentation[END_REF]. In this study, all data including the test and training volumes, are nonlinearly registered into an unbiased non-linear average template [START_REF] Fonov | BDCG: Unbiased average age-appropriate atlases for pediatric studies[END_REF] to minimize the size, orientation and position differences between subjects.

Non-local Means Patch-based Segmentation

The non-local means patch-based segmentation uses the local intensity distance between patches to estimated the weighting of manual segmentations in the training data, and the local label fusion is based on the weighted average of manual segmentations. The method described here is the same as in [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF]. In particularly, for each voxel x i in the test image, the weighting of each patch w(x i , x s,j ) is calculated by non-local means filter as:

w(x i , x s,j ) = e -p(x i )-p(x s,j ) 2 2 h 2 (3)
where p(x i ) is the cubic patch centered at x i , p(x s,j ) is the cubic patch centered at x j in the sth training subject and ||.|| 2 is the normalized intensity distance between two patches. The final segmentation is the weighted labels of all labeled samples inside the search volume v i of N training subjects.

v(x i ) = N s=1 j∈vi w(x i , x s,j )y s,j N s=1 j∈vi w(x i , x s,j ) (4)
where y s,j is the manual segmentation at Voxel x j of the sth training subject, and w(x i , x s,j ) is the weighting of the patch.

Combine Appearance Modeling and Non-local Means Patch in Segmentation

Appearance-model based global segmentation is good at capturing the global shape variation but not sensitive to the local shape change; while non-local means patch-based segmentation is sensitive to the local geometry, but doesn't have global shape constraints since it only uses the local intensity property in the segmentation. To take advantage of the global shape constraints and local intensity information, we combine these two methods together to improve the segmentation accuracy. For this purpose, we first apply appearance-model based segmentation to the test subject and then use the result from the global segmentation as input of the patch-based segmentation to locally refine the segmentation along the structure boundary area. A summary of the proposed algorithm is as follows.

- Instead of using the manual labels in Eq.4, the signed Euclidean distance functions of the manual labellings are integrated into the equation, i.e.

φ(x

i ) = N s=1 j∈vi w(x i , x s,j )φ s,j N s=1 j∈vi w(x i , x s,j ) (5)
-The final segmentation is achieved by threshing φ.

Illustrative Experiments and Results

The proposed segmentation algorithm is applied to segment human medial temporal lobe structures from real MR images. The validation uses the subset of the International Consortium for Brain Mapping (ICBM) database. The medial temporal lobe structures (HC, AG, EC, PRC, PHC) on T1w MRI data for 54 subjects from the ICBM database (Philips 1.5T gyroscan, TE=10ms, TR=30ms, ang=30 o , 1mm 3 isotropic voxels) were manually segmented using the protocol defined by Pruessner [27] [28]. The automatic segmentation results are compared with the manual labels. The similarity of these two labellings are measured by Dice kappa (κ) [START_REF] Dice | Measures of the amount of ecological association between species[END_REF]

[κ = 2 * (V (M ∩ A))/(V (M ) + V (A))].

Effects of Different Patch-size on Segmentation Performance

To study the impact of different patch size on segmentation accuracy, we segmented HC, para-hippocampal complex cortex (PCC=EC+PHC), and PRC using different patch sizes. The kappa results were presented in Fig. 1 which demonstrates that the best median kappa values were obtained with a patch size of 7 × 7 × 7 for all structures. These results indicated that too small patch size might not be able to capture the local geometry, and too big patch size might fail to find the best matched patches in the training data.

Effects of Different Search-size on Segmentation Performance

The impact of different search window size on segmentation accuracy was also analyzed by HC, PCC, and PRC. The κ values of 14 test subjects were shown in Fig. 2. From Fig. 2, we could see that the best median kappa values were obtained with a search window size of 5 × 5 × 5 or 7 × 7 × 7. Compared to the results in [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF], the optimal search window size was a litter bit smaller. However, the data in [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF] were aligned by linear registration and the data in this experiment were aligned by non-linear registration. Thus, this result suggested that the better alignment among subjects could reduce the search size.

Validation of Segmentation Accuracy on MTL Structures

The proposed method was applied to segment both left and right of HC, AG, PCC, PRC, and temporopolar cortex (TPC) from 54 MRI volumes using leaveone out method. The segmentation performance in terms of κ values based on the appearance-model based method with and without local correction was shown in Table 2. From these experiments, we can see: 1) for all structures, the individual κ value from the Appearance-based segmentation with local correction was always higher than that without local correction, and 2) the non-local means patch-based local correction increased the mean κ by 1.5% ∼ 3.5%, and those improvements were statistically significant (p ≤ 0.001). The improvement of the patch-based refinement on segmentation accuracy could also be visually observed in Fig. 3 which showed segmentation examples of all structures for three test subjects. The experiments also showed that the proposed method was able to quickly segment a new subject (less than 30 seconds for appearance-based global segmentation, and 1.5 minutes for patch-based refinement) on a 1.5GHz Linux PC once data was aligned (6 minutes per subject for nonlinear registration).

Discussion and Conclusions

In this paper, we present a novel segmentation algorithm which combines the appearance modeling and non-local means patch-based segmentation into a general segmentation framework, and applied it to the cortical structures of the medial temporal lobe. During segmentation, appearance modeling is used to capture the global shape variation, and non-local means patch-based method is used to improve the local fitting of the segmentation result. The experimental results have demonstrated feasibility, good performance, and robustness of this algorithm in 3D image segmentation. In the experiments, we tested the proposed method using the ICBM data set with 54 young healthy subjects. The leave-one out validation experiments of 54 subjects in 3D volumes demonstrated the segmentation accuracy (the mean κ of 0.87 for HC, 0.81 for AG, 0.73 for EPC, 0.73 for PRC, and 0.71 for TPC). The direct comparison between our technique and others in the literature is difficult because of the different anatomic definitions of the structures, different types of input data, and different qualities of manual segmentations. However, taking these caveats into consideration, our results are among the best results of previous publications in HC and AG (Detail see Table . 1). For other medial temporal lobe structures, there is no published result available except for PHC [START_REF] Heckemann | Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[END_REF]. Those promising initial results provide the impetus for future studies where the procedure will be applied to MR data from patients with temporal lobe epilepsy or Alzheimer's disease.

  Do appearance-model based global segmentation and obtain the segmented distance function φ. -Define local refinement area Ref ring which are the voxels inside the distance range [d 1 , d 2 ] of φ. -For each voxel x i inside Ref ring , re-calculate the patch-similarity function of φ(x i ) using non-local means patch-based method described in Sec 2.2.

Fig. 1 .

 1 Fig. 1. Impact of patch size on segmentation performance. κ values of 14 test subjects under different patch sizes. (a) Segmentation for HC (b) Segmentation for PCC, and (c) Segmentation for PRC.

Fig. 2 .

 2 Fig. 2. Impact of search window size on segmentation performance. κ values of 14 test subjects under different search window sizes. (a) Segmentation for HC (b) Segmentation for PCC, and (c) Segmentation for PRC.

Fig. 3 .

 3 Fig. 3. Visualization of segmentation results in 2D of three median test subjects: each row represents one test subject, and columns shows the segmentation results from (a) Manual labeling, (b) Appearance-based segmentation without local refinement, and (c) Appearance-based segmentation with local refinement. Note that the segmented labels of different structures are at the top of corresponding T1 test MR image, and presented by different colors: purple for HC, blue for AG, sky blue for PCC, white for PRC, and green for TPC.

Table 1 .

 1 Segmentation method review

	Author	Method summary	Result in term of κ
			HC	AG	PHC
	Fischl 2002 [9]	FreeSurfer	0.8	0.75∼0.78
	Heckemann 2006 [11]	Multi-atlas based	0.82	0.8	0.81*
	Chupin 2007 [6]	Seeding+region growing	0.84	0.76∼0.8
	Powell 2008

Table 2 .

 2 κ statistical result of segmentation Appearance-based segmentation Appearance-based segmentation Matched-pairs t-test < 0.001 < 0.001 Values in columns 2-5 are mean κ and the standard deviations are in parentheses.

	without local correction	with local correction	(p-value)
	Left	Right	Left	Right	Left	Right
	HC 0.851 (0.028)	0.862 (0.02)	0.867 (0.025)	0.873 (0.019)	< 0.001 < 0.001
	AG 0.80 (0.048)	0.792 (0.055)	0.812 (0.043)	0.803 (0.053)	< 0.001 < 0.001
	PCC 0.711 (0.068)	0.697 (0.083)	0.735 (0.066)	0.714 (0.082)	< 0.001 < 0.001
	PRC 0.696 (0.067)	0.707 (0.060)	0.730 (0.048)	0.739 (0.047)	< 0.001 < 0.001
	TPC 0.684 (0.065)	0.675 (0.048)	0.712 (0.057)	0.703 (0.047)