Acoustic Scattering Models from Rough Surfaces: A Brief Review and Recent Advances

This paper proposes a brief review of acoustic wave scattering models from rough surfaces. This review is intended to provide an up-to-date survey of the analytical approximate or semi-analytical methods that are encountered in acoustic scattering from random rough surfaces. Thus, this review focuses only on the scattering of acoustic waves and does not deal with the transmission through a rough interface of waves within a solid material. The main used approximations are classified here into two types: the two historical approximations (Kirchhoff approximation and the perturbation theory) and some sound propagation models more suitable for grazing observation angles on rough surfaces, such as the small slope approximation, the integral equation method and the parabolic equation. The use of the existing approximations in the scientific literature and their validity are highlighted. Rough surfaces with Gaussian height distribution are usually considered in the models hypotheses. Rather few comparisons between models and measurements have been found in the literature. Some new criteria have been recently determined for the validity of the Kirchhoff approximation, which is one of the most used models, owing to its implementation simplicity.


Introduction
The scattering of waves from rough surfaces is an important phenomenon in diverse areas of science, such as optics, geophysics, communications, acoustics, hydraulics [1], elastodynamics for non-destructive testing, etc. For instance, in the particular case of ultrasonic waves, scattering from rough surfaces leads to many applications in non-destructive testing [2], in underwater acoustics for high resolution ultrasonic sonars [3] and for medical sciences [4,5]. The modeling of such phenomena is a century-old and still active research topic.
Scattering models of acoustic waves from rough surfaces could be classified into analytical models, numerical ones and combinations of both analytical and numerical methods. However, such a classification is rather ambiguous, since some analytical models (usually called semi-analytical) employ numerical calculations and since some numerical methods start from analytical approximations. The current review focuses on analytical or semi-analytical methods and the validity of each method will be emphasized as well as some applications examples. Some analytical methods are presented or used in the following in more detail: the Kirchhoff approximation (KA), the perturbation theory (PT), the integral equation (IE) method, the parabolic equation (PE) method and the small slope approximation (SSA). The Kirchhoff and PT methods represent early classical approaches but are still very often employed in the literature, notably the Kirchhoff one. The other methods represent more modern approaches which have larger domains of validity but are more complex and therefore more difficult to implement or can lead to a higher calculation time. All the previously cited methods have

Statistical Properties
Surface roughness is described here as the variation of height z = h(x, y) above or below a mean surface. <h>, the average of the height variation over the surface, is by definition 0. h is assumed to be random and described statistically using parameters such as the RMS height and spatial correlation lengths. In the theory presented hereafter and developed by Ogilvy [6,27], it is assumed that the height distribution is well represented by a Gaussian.
The probability of any measured height being between h and h + dh above or below the mean plane is then given by p(h) dh, where: where σ is the RMS height: Random rough surfaces show some spatial correlation, as the height at a location is not totally independent of surrounding heights. This may be described by the correlation function of a surface, C(R) where: which represents the extent to which the surface height at any one point r is, on average, related to the surface height at point r + R away. Gaussian functions are often used as an approximation of the correlation of real surfaces. Assuming the same Gaussian correlation in every direction, C can be written: where L is called the correlation length of the surface. Anisotropic correlation functions can also be used to describe surface with a direction-dependent correlation. Such surfaces often occur in practice in industrial structures, as roughness can be generated in anisotropic manner during the manufacturing process.

Generation Process of an Individual Gaussian Rough Surface
The method of generating these surfaces is described in detail in [6,27]. It first generates uncorrelated random numbers with a Gaussian distribution of zero mean and specified RMS value. Then, it performs a moving average with weights that ensure that the correlation function is Gaussian with the desired correlation length. The following weights yield a surface with two independent lengths L x and L y in the two perpendicular directions: ∆x is the interval between adjacent surface points in the x direction and ∆y is the corresponding interval in the y direction. The surface correlation function is defined for a reference point (x 0 , y 0 ) by: The weights given in Equation (5) yield the following surface correlation function: C(x, y) = e −[x 2 /L 2 x +y 2 /L 2 y ] This method allows generating different realizations of random surfaces sharing the same statistical parameters (see Figure 7d for an application example).

Implementation Example for Generation of an Individual Gaussian Rough Surface
The parameters to be chosen as inputs for the parametrization of a Gaussian rough surface are: - The RMS surface height: σ = h 2 1 2 . As the height is defined with zero mean ( h = 0), the standard deviation σ = h − h 2 2 1 2 of the height is equal to the RMS height rms(h) = h 2 1 2 : -The correlation lengths, L x and L y , along the local axes x and y.
The number of mesh points of the rough surface along the local axes x and y has then to be automatically fixed by the chosen mesh algorithm. The previously described Gaussian rough surface model is currently implemented for planar and cylindrical surfaces in 3D computations, and for their linear cross-sections in 2D, in the CIVA software [28] (platform for non-destructive evaluation simulation [29][30][31][32][33]). Examples of generated surfaces for small specimens with relatively high roughness are shown in Figure 1.
This method allows generating different realizations of random surfaces sharing the same statistical parameters (see Figure d for an application example).

Implementation Example for Generation of an Individual Gaussian Rough Surface
The parameters to be chosen as inputs for the parametrization of a Gaussian rough surface are: - The RMS surface height: The correlation lengths, x L and y L , along the local axes x and y.
The number of mesh points of the rough surface along the local axes x and y has then to be automatically fixed by the chosen mesh algorithm. The previously described Gaussian rough surface model is currently implemented for planar and cylindrical surfaces in 3D computations, and for their linear cross-sections in 2D, in the CIVA software [28] (platform for non-destructive evaluation simulation [29][30][31][32][33]). Examples of generated surfaces for small specimens with relatively high roughness are shown in Figure 1.  Figure 6) and   Figure 6) and L/σ = 4.3 (case of Figure 11); (d) planar for L/σ = 1; (e) cylindrical for L/σ = 1. Scales in mm.

Historical Approximations
This section deals with the perturbation and Kirchhoff approximations. Both have been extensively applied to rough surfaces and have their own strengths and weaknesses. Other approaches will be dealt with in later sections. Figure 2 shows the geometry used in this section for rough surface scattering. We consider a rough surface whose roughness is defined by z = h(x,y) where h is the height variation above or below the mean plane, z = 0. The xz plane contains the incident wavevector.

The Perturbation Theory
The perturbation theory (PT) assumes some restrictive hypotheses on the local height h and gradient ∇h at points of a rough surface: k|h| << 1 (8) |∇h| << 1 (9) k being the incident wavevector modulus. It is assumed that if previous conditions are satisfied, the roughness is weak compared to acoustic wavelength and slowly varies. It expresses the scattered field at an observation point r as an expansion series: where ψ inc is the incident wave field and ψ sc n is the nth-order approximation to the scattered field. The n = 0 term is the scattered field from a smooth surface. Then the main principle of the perturbation Theory consists in assuming that the boundary conditions may be expanded in a Taylor series about the mean plane z = 0. Mathematically, the boundary conditions on the rough surface expressed as: where f can be, for instance, the acoustical pressure for Dirichlet conditions-which can be expanded in a Taylor series about the mean plane z = 0 as follows: Some applications of the perturbation theory consider first scattering from canonical rough surfaces [18]. For instance for the Lamb's problem (scattering from a surface excited by a load force applied at one particular surface point) [34], the solution is expressed as convolution integrals evaluated only for a triangular irregularity of small height compared to its width and to the wavelength. The application of PT to random rough surfaces leads to consider separately a mean coherent contribution and an incoherent one. For the coherent part, reflection and transmission coefficients are determined. For the incoherent part, the mean intensity of the scattered field is calculated [35]. The determined first-order PT term can be substituted back into the Helmholtz scattering formula to iteratively deduce higher order terms.
Blakemore [36] has then presented a Perturbation Theory at the first order for a fluid/solid interface. PT simulations in backscattering for rigid and elastic solids are very close except at the vicinity of the Rayleigh angle as shown in Figure 3. coefficients are determined. For the incoherent part, the mean intensity of the scattered field is calculated [35]. The determined first-order PT term can be substituted back into the Helmholtz scattering formula to iteratively deduce higher order terms.
Blakemore [36] has then presented a Perturbation Theory at the first order for a fluid/solid interface. PT simulations in backscattering for rigid and elastic solids are very close except at the vicinity of the Rayleigh angle as shown in Figure 3.  [36]. Backscattered intensity with the perturbation theory on a steel plate immersed in water. The used parameters are respectively frequency = 2.25 MHz, angular beam width = 3°, observation distance = 300 mm, RMS surface height  = 10 µm and the two correlation lengths of the rough surface both equal to 100 µm. Reproduced with permission from [36], Elsevier, 1993. In this region a small peak in the intensity of the incoherently scattered field is predicted, consistent with experiment. Nevertheless, in the case of the proposed experimental validation in pulse echo configuration [37] (Figure 4), Kirchhoff prediction (not shown here) are about 3 dB lower than PT one at backscattering angles theta = 40°, and 5.5 dB lower at theta = 70°.  [36]. Backscattered intensity with the perturbation theory on a steel plate immersed in water. The used parameters are respectively frequency = 2.25 MHz, angular beam width = 3 • , observation distance = 300 mm, RMS surface height σ = 10 µm and the two correlation lengths of the rough surface both equal to 100 µm. Reproduced with permission from [36], Elsevier, 1993. In this region a small peak in the intensity of the incoherently scattered field is predicted, consistent with experiment. Nevertheless, in the case of the proposed experimental validation in pulse echo configuration [37] (Figure 4   Some applications and references of PT for scattering from fluid/fluid, fluid/solid interfaces and fluid/stratified layers [38][39][40] for inhomogeneous media [41,42] can be found in the literature, generally for ocean applications. An acoustic model in the time domain based on first-order perturbation theory for scattering from layered sediments has recently been proposed [43].
Only a few authors have been interested in the scattering from canonical non-planar rough targets. Notably Bjorno [44,45] recently investigated the scattering problem for an elastic rough sphere. He developed in that case a second-order perturbation model. He compared the scattering  [36]. Experimental validation of Blakemore's perturbation theory at first order: at left measurements from de Billy et al. [37] on a brass plate immersed in water; at right, Blakemore's results with frequency = 2.2 MHz, angular beam width = 3 • , observation distance = 300 mm, RMS surface height σ = 14.9 µm and the correlation lengths L x = 90 µm and L y = 105 µm. Reproduced with permission from [36], Elsevier, 1993.
Some applications and references of PT for scattering from fluid/fluid, fluid/solid interfaces and fluid/stratified layers [38][39][40] for inhomogeneous media [41,42] can be found in the literature, generally for ocean applications. An acoustic model in the time domain based on first-order perturbation theory for scattering from layered sediments has recently been proposed [43].
Only a few authors have been interested in the scattering from canonical non-planar rough targets. Notably Bjorno [44,45] recently investigated the scattering problem for an elastic rough sphere. He developed in that case a second-order perturbation model. He compared the scattering patterns obtained versus the adimensional wave-number ka and for different values of the ratio of the rms height to the sphere radius with its perturbation solution to the classical acoustic scattering solution for a smooth sphere [46,47]. Measurements of the backscattering response of a cast iron sphere versus ka led to a good experimental validation of the proposed model: the roughness yields a decrease in the backscattered amplitude notably for large ka values. The accuracy is limited to small values of the ratio rms height to wavelength.
Note also, for instance, a recent development in electromagnetism of the perturbation theory for a slightly deformed cylinder [48].
We also have to cite the phase perturbation method [49,50] in which the phase (rather than the field) is expanded in a series of different orders in kh. This theory exhibits links with both PT and Kirchhoff [51]. A validation [52] of the phase perturbation method has been proposed and has consisted in comparing its results to numerical ones in a region where PT and Kirchhoff are found ineffective. The phase perturbation method is found to work well in such a region except for low grazing observation angles. However, applications of this new theory seem not numerous and are more focused on radar cross sections calculations [53].
Although the perturbation approach can take into account of the scattering and re-radiation of Rayleigh waves, it appears possibly more restrictive than the Kirchhoff method in its range of applicability. The perturbation method is notably an approach valid for low roughness (when compared to acoustic wavelength). Several validations of the PT are provided in the next sections; they also show that PT is rather accurate far from the specular direction.

Theory and History
The main hypothesis of the Kirchhoff approximation (also referred to as the physical optics approximation) is to assume that each point of the surface behaves as if it belonged to the infinite tangent plane to this surface. Then it is practical to employ geometrical acoustics to approximate the exact fields on the local surface: the field on the illuminated part of surface is the sum of the incident field and the field reflected by the local tangential plane. This supposed surface field is then propagated far from the surface using a propagator (called Green function). Indeed the scattered far field is expressed by the way of an integral over the flaw surface (the Kirchhoff-Rayleigh integral) whose integrand is-for the example of Dirichlet conditions-the product of the Green function by the derivative of the geometrical field on the flaw surface. For complex incident fields or scatterer shapes, this integral is evaluated numerically.
The Kirchhoff approximation has been extensively used to model scattering of waves in optics (early introduced by Meecham [54]), electromagnetism [55,56], acoustics, etc. In both acoustics [57] and elastodynamics [58,59], the Kirchhoff approximation also has the advantage in dealing with anisotropic materials and impedant surfaces [57].
Note that when doing an asymptotic evaluation of the KA integral on a finite surface, the KA resulting leading terms describe the reflected waves, and the other terms describe diffracted waves by the surface edges which have the same nature (cylindrical or conical wave fronts) but not the same amplitude as those predicted by the geometrical theory of diffraction (GTD) [55,60]. The main advantage of the Kirchhoff approximation is to provide finite results for the scattered field in the whole space, contrary to GTD [61][62][63][64][65]. One of the main inconveniences of the Kirchhoff approximation (KA) is its incorrect quantitative prediction of edges diffraction, since the approximation of the surface field by the geometrical field is invalid near the edges. To handle this drawback, the physical theory of diffraction (PTD) has been proposed in electromagnetism [66], acoustics [59,67] and elastodynamics [68][69][70] and compared to other uniform scattering models [67,70,71].
In the literature, the Kirchhoff approximation has been studied by different authors. A review has been carried out by Ogilvy [18]. The first one to introduce the Kirchhoff approximation is Eckart [72].
Ogilvy [27] has well described the theory for 1D rough profile ( Figure 2) and Dirichlet conditions (soft interface) extending Eckart's works. De Billy et al. proposed a Kirchhoff-based theory extended to shadowing effects for the backscattering of acoustic waves by randomly rough surfaces [37]. De Billy et al. employed in fact the potential-method formulation developed by Welton which showed by iterating the scattering integral equation that the scattered pressure can be expressed as a series: the first term corresponds to the singly and the others to the multiply scattered field components. De Billy et al. considered only the first term. Their experimental validations of their model using Gaussian rough (liquid-solid interfaces) surfaces have shown that the experimental measurements of the intensity backscattered by randomly rough interfaces (liquid-solid) are in good agreement with the theory until 40 • backscattering angles (with respect to the normal to the mean plane). Nevertheless, an anomaly has been observed in Kirchhoff's predictions for the Rayleigh's angle of incidence. In fact, to model the impedance surfaces used in their experimental validations, they simply add a reflection coefficient as a factor inside the Kirchhoff integral, and for backscattering the local reflection coefficient is the one at normal incidence for a plane surface: this formalism does not take into account any mode conversion in the elastic medium and is unable to predict any perturbation at the different critical angles. Following de Billy's works which studied the Kirchhoff approximation using monochromatic plane waves, Imbert [73] has proposed a Kirchhoff model in time domain for 2D rough profiles and Neumann conditions (rigid interface).
We can sum up the main steps of Kirchhoff approximation application as follows: 1.
Use of the Helmholtz equation: the scattered field is: The equation is also called Green theorem and G(r, r 0 ) = e ik|r−r 0 | 4π|r−r 0 | is called the Green function, where r 0 is a current point on the scattering rough surface S 0 and r the observation point. The field in the fluid medium can be obtained thanks to the knowledge of its value on the surface S 0 .

2.
Boundary Conditions on the rough surface: Ogilvy [27] considered for instance Dirichlet boundary conditions (soft case) so that ψ(r 0 ) = 0. So one term in the previous Green integral disappears: 3.
Kirchhoff approximation: the tangential plane approximation enables to approximate the value of ∂ψ(r 0 ) ∂n 0 inside the previous integral.

4.
Far field approximation: the integrand value can be simplified by assuming kr >> 1, k being the wave number.

5.
Projection of integration on the rough surface on the smooth surface: this projection involves the gradients of the roughness. 6.
Incorporation of a shadow function to approach auto-shadowing: Wagner [74] has proposed to incorporate a shadow function inside the Kirchhoff integral to take into account shadowing effects.
Note that the Kirchhoff approximation is mainly used to predict far-field scattering; a model has been proposed to evaluate the scattering pattern in near-field [75] after the interaction of an acoustic wave with a rough interface.

Validity
Historically KA has been considered-in the point of view of its local tangential plane assumption-as accurate when the rough surface can be seen planar by the incident wave. In Ogilvy's review [7], KA is thus said valid for kr c cos 3 θi >> 1, with r c the surface local radius of curvature and θi the incident angle. This historical criterion forgets that multiple scattering and shadowing effects are not taken into account in KA leading to large errors for grazing incidences or higher roughness.
KA prediction has been compared to Integral Equation (IE) results by Thorsos [21] for 1D rough profiles. Thorsos proposed to use the surface correlation length L as main parameter for defining the Kirchhoff validity away from the low grazing angle region: KA is said to be accurate for a correlation length greater than the wavelength or even for lower correlation lengths if shadowing effects are not important. We notice this first criterion proposed by Thorsos at this step is too simple and not involve the rms height. In the low grazing angle region some examples of validation with comparison to another methods are given in next Section 4; at low grazing angles, additional validity dependence enters through the RMS slope which is for Gaussian rough surfaces: s = tan γ = √ 2 σ L i.e., the rms height is now involved. However, at shorter correlation lengths, the needed shadowing correction proposed by Thorsos becomes greater to improve Kirchhoff prediction at low grazing angles in the case of backscattering.
Note that in order to improve its prediction for low grazing angles, Chapman [76] provided an improved Kirchhoff formula (using a smoothing special function) notably by including the correlation length of the scattering surface into the reflection coefficient. The modification impact is not noticeable since shadowing and multiple scattering will dominate for low grazing angles. Chapman showed that the PT and KA results differ except in the limit of infinite correlation length (smooth surface).
Shi et al. [77] recently proposed for elastic materials (also for shear waves [78,79]) and for compressional waves a different surface parameter than the surface correlation length L used by Thorsos, a function σ 2 /L of both L and σ, to define a criterion when the KA use is valid. They found that, with a normal incidence angle θi = 0 • , KA is valid with a scattering angle −70 • < θs < 70 • and when σ 2 /L ≤ c, the constant upper bound c being in two dimensions 0.20 λp and in three dimensions 0.14 λp (λp being the compressional wavelength). Figure 5 shows a comparison between KA and a finite element (FE) method for a 2D case when σ 2 /L~c. A small incidence angle of 30 • can improve the KA accuracy when σ 2 /L exceeds c as in the case shown in Figure 6. We can however regret that the numerical evaluation of this proposed criterion was limited to small incidences angles: however, this criterion is likely to be suitable for many NDE configurations for which incidence is close to normal. A study of this criterion near grazing incidence should also be useful.
As said in the introduction, another different KA validity criterion has been proposed for radar applications [22] and is also based on numerical calculations: it is given by Equation (1) in [22] (see also their Figure 3, Kirchhoff being called the scalar approximation) and ensures a KA validity for L 2 /σ more than a constant times the wavelength and a condition kL > 6 is also associated. We find the condition on L 2 /σ very disbelieving and in contradiction with the previously cited studies and also that given by Voronovich (σ/L < C and large kL; see Figure 1.1 in [20]). The condition kL >6 is also found in other references as [52]. Notably, Figure 2 in [52] represents the stated validity of first-order perturbation theory and KA in kσ-kL space obtained for bistatic scattering cross sections simulation: PT is said valid for small kσ < 2 to 4 depending on kL value, KA valid for kL > 6 and both KA and first-order PT needs σ/L < C (C being a constant value).
In our opinion, this KA condition kL > 6 is too restrictive for many applications, such as non-destructive evaluation, which often employs near specular observation inspection. First, we previously study or carry out numerical comparisons of KA simulation for cracks of finite size L inspected with longitudinal waves in solids with numerical methods: the method of optimal truncation (MOOT) [80] and a finite element method (FEM) [81,82]. These comparisons with MOOT [83] and FEM [69,84] both show that KA is valid for near specular observation inspection for kL > 2. Shi et al. [77] did not propose a validity criterion on kL (but only the condition on σ 2 /L); they nevertheless obtained an acceptable disagreement of 2dB (see their Figure 10 [77]) for kL = π/2 (corresponding to L = λp/4) and for an important roughness σ = λp/4 = L equal to the correlation length. It is a pity that that they did not provide numerical validations below kL = π/2. 4 λp (λp being the compressional wavelength). Figure 5 shows a comparison between KA and ite element (FE) method for a 2D case when σ 2 /L~c. A small incidence angle of 30° can improve th accuracy when σ 2 /L exceeds c as in the case shown in Figure 6. We can however regret that th merical evaluation of this proposed criterion was limited to small incidences angles: however, th iterion is likely to be suitable for many NDE configurations for which incidence is close to norma study of this criterion near grazing incidence should also be useful.   Chapman showed that the PT and KA results differ except in the limit of infinite correlation length (smooth surface). Shi et al. [77] recently proposed for elastic materials (also for shear waves [78,79]) and for compressional waves a different surface parameter than the surface correlation length used by Thorsos, a function σ 2 /L of both L and σ, to define a criterion when the KA use is valid. They found that, with a normal incidence angle θi = 0°, KA is valid with a scattering angle −70° < θs < 70°and when σ 2 /L ≤ c, the constant upper bound c being in two dimensions 0.20 λp and in three dimensions 0.14 λp (λp being the compressional wavelength). Figure 5 shows a comparison between KA and a finite element (FE) method for a 2D case when σ 2 /L~c. A small incidence angle of 30° can improve the KA accuracy when σ 2 /L exceeds c as in the case shown in Figure 6. We can however regret that the numerical evaluation of this proposed criterion was limited to small incidences angles: however, this criterion is likely to be suitable for many NDE configurations for which incidence is close to normal. A study of this criterion near grazing incidence should also be useful.   Our conclusive opinion is that for many applications using rather near specular observation and non-grazing incident angles, two criteria can be used for KA validity: σ 2 /L ≤ Cλp and kL > 1 or 2. For more important observation angles to the normal the criterion on kL and the constant C value are likely to be a little more restrictive: Zhang [85] using 2D FEM calculations recently showed that KA can be considered accurate for kL > π and for incidence and scattering angles over the range from −80 • to 80 • .

Elastic Targets and Applications
The Kirchhoff approximation models bulk waves, and does not take into account surface acoustics waves (SAW) as Rayleigh waves neither head waves. Head waves led to recent researches [70,[86][87][88] and a modification [89] of the Kirchhoff approximation is henceforth proposed for head waves account. Head waves on irregular surfaces can also be modelled [90].
The Kirchhoff approximation for elastic rough targets is proposed by Ogilvy [13] and then revisited by Bjorno [91] and Dacol [92]. The Kirchhoff approximation validity have also been specifically considered in the particular cases of rigid and soft targets in a modeling study [93] of the reflected field including multiple scattering. One example of KA application to impedance surfaces [57] for non-destructive evaluation (NDE) is the modeling of ultrasonic telemetry of immersed targets [59,94].
We report here an application of KA simulation to impedance rough surfaces for NDE. Ultrasonic inspection of a planar steel component immersed in water (water path 25 mm) at frequency = 2 MHz (wavelength λ = 0.75 mm) is simulated with the Kirchhoff approximation for elastic rough entry surfaces of varying roughness using CIVA software inspection simulations tools [57,95,96]. A phased-array (represented in yellow in Figure 7a,b) is used to emit an incident beam directed at 117 • with respect to the surface thanks an emission delay law (applied in red in Figure 7a to only some left elements) and all elements receives the scattered waves without delay law in the angular domain (118 • , 58 • ).
In Figure 7a,b, the true B-scans (true-to-geometry imaging) are surperimposed to the specimen in color codes; they represent the echographic signal versus time for all elements in reception (corresponding to scattering angles in the angular domain (118 • , 58 • )). Figure 7c shows that increasing the roughness deteriorates the specular reflection echo and causes appearance of echoes at non-specular angles. Figure 7d shows that the scattering directivity patterns obtained for three different generated rough surfaces of the same statistical properties are slightly different.
A recent study [97] has compared the accuracy of KA and PT using a finite element (FE) model for fluid and elastic ocean bottoms. For a fluid bottom, KA is close to FE with maximal gap of 2 dB except for small grazing incident angles for which PT is more effective. PT behaves globally well except near specular observation direction notably for small grazing incident angles. For the elastic case, the disagreement between KA and FE increases to 5 dB and KA is inaccurate for incident angles less than the longitudinal critical angle: the interferences between transversal, longitudinal and surface waves are not taken into account in KA. The PT prediction is generally good, but another limitation appears in the elastic case for incident angles near the longitudinal critical one. However, we regret the provided validations are established only for a small roughness k σ= 0.6 favorable to PT.

Periodic Surfaces
As it will be discussed in the Section 5, Perspectives, it can also be useful to model scattering from periodic surfaces-such surfaces can notably appear in additive manufacturing parts or in oceans. For the case of periodic surfaces, one exact solution is the classical Rayleigh-Fourier method (RFM) [98]. Impinged by a plane wave, a periodic surface leads to Bragg scattering with a finite number of plane wave propagation angles in far field. A very recent study evaluates the validity of the Kirchhoff approximation (KA) for the scattering problem from a sinusoidal surface-mimicking the ocean one-with comparison to the RFM solution [99]. An asymptotic evaluation of the Kirchhoff integral which is consequently equivalent to geometrical acoustics [100] provides a ray theory which is also evaluated to better understand the Kirchhoff prediction. This notably shows that Kirchhoff has the advantage to exhibit a smooth behaviour near caustics generated by the sinusoidal surface contrary to geometrical acoustics. The KA validity is also studied versus the observation range (projected distance between the source and receiver points): for a short range, KA is very close to RFM. At moderate range, KA sometimes over-estimated the peak pressure due to interferences between several rays. For longer ranges, the KA disagreement increases strongly. This study [99] is very useful for better understanding the KA limitations due to multiple paths arrivals and grazing incidence/observation. Nevertheless, KA remains a very useful tool even for periodic surfaces when source and observations are not too far. In Figure 7a,b, the true B-scans (true-to-geometry imaging) are surperimposed to the specimen in color codes; they represent the echographic signal versus time for all elements in reception (corresponding to scattering angles in the angular domain (118°, 58°)). Figure 7c shows that increasing the roughness deteriorates the specular reflection echo and causes appearance of echoes at non-specular angles. Figure 7d shows that the scattering directivity patterns obtained for three different generated rough surfaces of the same statistical properties are slightly different.
A recent study [97] has compared the accuracy of KA and PT using a finite element (FE) model for fluid and elastic ocean bottoms. For a fluid bottom, KA is close to FE with maximal gap of 2 dB except for small grazing incident angles for which PT is more effective. PT behaves globally well

More Modern Approaches with Larger Domains of Validity
We have presented the two classical approximations: the perturbation theory (PT) and the Kirchhoff approximation (KA). In the view of main application to scattering from ocean, it is important to cite the two-scale roughness theories [101][102][103][104]. They consist in assuming that scattering can be considered as the combination of two phenomena: that produced by large facets and well approximated by KA and that due to small roughness and approximated by PT. One major inconvenient of this technique is that it requires to fix a wavenumber parameter [105] to separate the two scattering domains.
The small slope approximation, described in the next section, is a much more convenient technique to link the classical approximations and it reduces to their results in some limiting cases.
In this section, some sound propagation models which have been notably devised for underwater acoustic applications and providing a better validity notably for grazing observation angles are studied. Underwater acoustic propagation depends on many factors. When sound propagates underwater the sound intensity decreases when increasing the propagation distance. Propagation loss is a quantitative measure of the reduction in sound intensity between two points, normally the sound source and a distant receiver. Modeling the acoustic reflection loss [106] at the rough ocean surface needs simulating forward scattering for incident grazing angles and is of interest for sonar especially for shallow oceans or oceans with a surface duct (sound waves trapped in a duct layer below the sea surface) [106]. Several models described hereafter (small slope approximation-SSA, and parabolic equation-PE) were notably designed to handle grazing angles. The integral equation (IE) method is also used as reference solution for validating SSA and comparing it to other approximations as KA and geometrical optics (GO).

Small Slope Approximation (SSA)
The SSA technique was described first by Voronovich [107][108][109]. An alternative small slope approximation has originally been developed by Urusovskii and then studied by McDaniel (1994) [110]. All these methods attempt to exploit the assumed smallness of the surface slope without restricting the RMS surface height to wavelength ratio.
The SSA model from Voronovich has been investigated by Thorsos and Broschat [111,112]. SSA consists in introducing a roughness correction to the Helmholtz integral. This correction is expanded into a functional Taylor series in slope whose terms depend on the integral of the wave number spectrum of the surface height function. The approximation of each order of SSA is based on the approximations of previous orders. Owing to its principles SSA takes into account shadowing and diffraction.
Thorsos and Broschat showed that, at lowest order the small slope approximation of Urusovskii/McDaniel is limited to forward scattering [110] but that the lowest order term in Voronovich's SSA was incorrect in McDaniel's SSA.
SSA tends to the perturbation approximation for small k σ (adimensional wave numbers in terms of RMS surface height σ). The second order SSA (SSA (2)) reduces to the Kirchhoff approximation (KA) at high frequencies, hence does not take into consideration the shadowing effects. The predictions of SSA(2), SSA(3), and SSA(4) are compared [112] with Monte Carlo integral equation (IE) calculations devised earlier by Thorsos [21], and with other approximations: the perturbation approximation, the Kirchhoff approximation (KA), and geometrical optics (GO). The scattering geometry is described in Figure 8.

Small Slope Approximation (SSA)
The SSA technique was described first by Voronovich [107][108][109]. An alternative small slope approximation has originally been developed by Urusovskii and then studied by McDaniel (1994) [110]. All these methods attempt to exploit the assumed smallness of the surface slope without restricting the RMS surface height to wavelength ratio.
The SSA model from Voronovich has been investigated by Thorsos and Broschat [111,112]. SSA consists in introducing a roughness correction to the Helmholtz integral. This correction is expanded into a functional Taylor series in slope whose terms depend on the integral of the wave number spectrum of the surface height function. The approximation of each order of SSA is based on the approximations of previous orders. Owing to its principles SSA takes into account shadowing and diffraction.
Thorsos and Broschat showed that, at lowest order the small slope approximation of Urusovskii/McDaniel is limited to forward scattering [110] but that the lowest order term in Voronovich's SSA was incorrect in McDaniel's SSA.
SSA tends to the perturbation approximation for small k  (adimensional wave numbers in terms of RMS surface height  ). The second order SSA (SSA (2)) reduces to the Kirchhoff approximation (KA) at high frequencies, hence does not take into consideration the shadowing effects. The predictions of SSA(2), SSA(3), and SSA(4) are compared [112] with Monte Carlo integral equation (IE) calculations devised earlier by Thorsos [21], and with other approximations: the perturbation approximation, the Kirchhoff approximation (KA), and geometrical optics (GO). The scattering geometry is described in Figure 1.  Figure 1 [112]. Scattering geometry. Reproduced with permission from [112], AIP Publishing, 1997.
SSA and IE scattering predictions are presented in Figure 9 for respective k σ and kL 0.5 and 4.0, for 45 • incident angle and RMS slope s = 0.177 (s = √ 2 σ L for Gaussian surface), and the RMS slope angle γ = tan −1 s is 10 • . The contribution of specular reflection-leading to a peak in the IE prediction-has not been integrated in the SSA models. In Figure 10, k  and kL have both been increased by a factor of three to 1.5 and 12.0, respectively. The increase in the surface roughness (resp. correlation length) causes disappearance of the specular peak (resp. decrease of the backscattered strength). A higher SSA order improves the prediction in backscattering. According to Thorsos [21], the Kirchhoff approximation KA is accurate away from low grazing angles for rough surfaces smooth on the scale of a wavelength, assuming small surface slopes. For a  [112]. Comparison of the second-order SSA (2), third-order SSA(3), and fourth-order SSA(4) small slope approximation and Monte Carlo integral (IE) scattering strengths (=10 log of the scattering cross section) for a modest rms surface slope angle. Here, k σ = 0.5, kL = 4.0, the RMS slope angle γ is 10 • , and the incident angle θ i is 45 • . Reproduced with permission from [112], AIP Publishing, 1997.
In Figure 10, k σ and kL have both been increased by a factor of three to 1.5 and 12.0, respectively. The increase in the surface roughness (resp. correlation length) causes disappearance of the specular peak (resp. decrease of the backscattered strength). A higher SSA order improves the prediction in backscattering. In Figure 10, k  and kL have both been increased by a factor of three to 1.5 and 12.0, respectively. The increase in the surface roughness (resp. correlation length) causes disappearance of the specular peak (resp. decrease of the backscattered strength). A higher SSA order improves the prediction in backscattering. According to Thorsos [21], the Kirchhoff approximation KA is accurate away from low grazing angles for rough surfaces smooth on the scale of a wavelength, assuming small surface slopes. For a Gaussian roughness spectrum, the KA is said to be generally accurate for kL >6, that is when 1 L   According to Thorsos [21], the Kirchhoff approximation KA is accurate away from low grazing angles for rough surfaces smooth on the scale of a wavelength, assuming small surface slopes. For a Gaussian roughness spectrum, the KA is said to be generally accurate for kL >6, that is when L λ > 1. We have seen in Section 3.2.2 the recent definition of a more precise criterion for the KA validity. In the following of this section the SSA and KA are thus considered for several examples with kL > 6.
In Thorsos' previous Figure 3 [112], KA is valid away from grazing angles, whereas PT (not shown) is not. The SSA(2) result reduces to the KA result in the specular region; the third and fourth-order results do as well. KA behaves better than SSA for backscattering, but it is the contrary in the forward direction at grazing angles.
Finally, several cases have been examined beyond the region of validity of the fourth-order Perturbation Theory PT(4) and of KA. For example, in their Figure 6 [112] since k σ = 1, PT(4) breaks down whereas KA is globally correct and SSA(4) the more realistic (Figure 11).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 27 order results do as well. KA behaves better than SSA for backscattering, but it is the contrary in the forward direction at grazing angles. Finally, several cases have been examined beyond the region of validity of the fourth-order Perturbation Theory PT(4) and of KA. For example, in their Figure 6 [112] since k = 1, PT(4) breaks down whereas KA is globally correct and SSA(4) the more realistic ( Figure 11). Other validations showed that the fourth-order SSA (SSA(4)) mimics some shadowing phenomena and the SSA results are globally accurate for all observation angles for incident angles up to 80°, but the validity decreases with the incident angle.
Another study [113] of the Bragg scattering from surface waves by a rough interface recently provided comparisons between the Kirchhoff approximation, the perturbation theory and the integral equation method; Kirchhoff leads to results relatively close to the integral equation method in the studied configurations (in forward and backscattering).
Voronovich (see Figure 1.1 in [20]) argued that SSA is valid for σ/L less than a constant value. According to him, KA and PT also require this condition but PT needs in addition a small kσ and KA a great kL.
SSA has also been developed for rough solid elastic interfaces [114]. Berman [115] showed that it provides an effective solution compared to Monte-Carlo integral equations simulations for grazing incident angles, whereas Kirchhoff leads to an overall important underestimation. SSA and firstorder PT give rather good predictions and exhibit minima or peaks at the Rayleigh angle that can be exaggerated or discontinuities at critical angles which are not physical. The authors noted that Kirchhoff is rather smooth or does not have the valid behaviour for these particular angles: it is logical because the associated phenomena (Rayleigh or head waves) are not accounted on the simplified Kirchhoff theory (see also the perspectives section). We nevertheless notice that Kirchhoff yields an overall rather good prediction versus scattering angle despite some quantitative discrepancies. We also emphasise that one comparison shows an overestimation by first-order PT of the specular Other validations showed that the fourth-order SSA (SSA(4)) mimics some shadowing phenomena and the SSA results are globally accurate for all observation angles for incident angles up to 80 • , but the validity decreases with the incident angle.
Another study [113] of the Bragg scattering from surface waves by a rough interface recently provided comparisons between the Kirchhoff approximation, the perturbation theory and the integral equation method; Kirchhoff leads to results relatively close to the integral equation method in the studied configurations (in forward and backscattering).
Voronovich (see Figure 1.1 in [20]) argued that SSA is valid for σ/L less than a constant value. According to him, KA and PT also require this condition but PT needs in addition a small kσ and KA a great kL.
SSA has also been developed for rough solid elastic interfaces [114]. Berman [115] showed that it provides an effective solution compared to Monte-Carlo integral equations simulations for grazing incident angles, whereas Kirchhoff leads to an overall important underestimation. SSA and first-order PT give rather good predictions and exhibit minima or peaks at the Rayleigh angle that can be exaggerated or discontinuities at critical angles which are not physical. The authors noted that Kirchhoff is rather smooth or does not have the valid behaviour for these particular angles: it is logical because the associated phenomena (Rayleigh or head waves) are not accounted on the simplified Kirchhoff theory (see also the perspectives section). We nevertheless notice that Kirchhoff yields an overall rather good prediction versus scattering angle despite some quantitative discrepancies. We also emphasise that one comparison shows an overestimation by first-order PT of the specular amplitude, which is a strong inconvenient for numerous applications as NDE.
SSA for layered fluid seafloors has only been studied recently. A study published this year [116] reviewed three different approximations. Jackson [117] has proposed an acoustic SSA method for rough water-sediment interfaces with a fluid layer of finite thickness overlying a semi-infinite fluid layer. A second method consists in using in acoustics the SSA developed for the general layered case in electromagnetics [118][119][120]. A third approach is developed in [116] by translating usual small-slope ansatz from k-space to coordinate space. It is then shown in [116] using different examples that the different approximations give different results in the case of strong internal reflections and also differ from Kirchhoff in specular direction. In an example with a mud layer overlying a harder basement, the first approximation differs strongly from the two others and consequently seems erroneous. Unfortunately, no numerical validation is provided to define the best ansatz.

Parabolic Equation (PE)
Acoustic losses at the rough ocean surface can be simulated utilizing a propagation code which can use as input an explicitly defined roughness profile. The PE (Parabolic equation [121][122][123]) model Ramsurf [124] is often used and enables to model shadowing from the surface roughness. Ramsurf is based on the split-step Padé method [121] (a pseudo-spectral numerical method for solving nonlinear partial differential equations).
Firstly, based on the rough ocean surface of Gaussian spectrum, the validity of Ramsurf acoustic propagation model has been studied [125] by comparison with measurements after an explosion at receivers located on the ocean bottom around 1000 m depth. The height for the sea surface is from 1.5 m to 2 m. Transmission loss curve under a smooth sea surface and that under the rough sea surface are simulated respectively by Ramsurf. The sea surface reflection loss caused by the surface roughness is about 4 to 6 dB. The transmission loss simulated by Ramsurf model is very close to the experimental one. Kirchhoff approximation (KA) and small slope approximation (SSA) abilities to simulate the reflection loss are also compared: SSA and Ramsurf match well for very small grazing angles (<6 • ) for which KA breaks down. We nevertheless notice that the gap between KA and Ramsurf becomes less than 3 dB above 6 • : KA is thus inaccurate only in a very small region.
Another study [126] focuses on the modeling of underwater acoustical reflection from a wind-roughened ocean surface and compares the surface loss between KA, a small-slope model SSA(2) and a stochastic modeling using Ramsurf. SSA(2) matches well with Ramsurf for the entire range of frequency and wind speed combinations studied (1.5-9 kHz, 5-12.5 m/s) and for sound incidence grazing angles from 1 to 11 • . KA can be used effectively from an angle about 2.5 to 7 • .
An extension of the parabolic equation method to deal with forward scattering has recently been proposed for both Dirichlet and Neumann boundary conditions [127] but future works are needed to make the connection between forward and backward 3D scattering and treat other boundary conditions.

Summary, Discussion and Perspectives
A review of acoustic wave scattering from rough surfaces has been done. Several analytical models proposed in the literature have been described in more detail as well as their validity. We are also focused in specific acoustic applications as scattering from ocean, non-destruction evaluation, medical sciences.
Two classical approximations have been early devised: the perturbation theory and the Kirchhoff approximation.
The perturbation theory (PT) consists in expressing the scattered field in an expansion series using as small parameter the roughness height. This method leads to a major inconvenient: it is valid for low roughness when compared to acoustic wavelength (at first order for small kσ < 2 to 4 depending on kL value, σ, L and k being respectively the rms roughness height, the correlation height and the wave number). Existing validations also showed that PT is rather accurate far from the specular direction. Consequently, it appears possibly more restrictive than the Kirchhoff method for numerous applications notably non-destruction evaluation (NDE) which usually employs inspection in specular direction.
The Kirchhoff approximation (KA) leads to a simple analytical expression; it is a relatively powerful and inexpensive approximation in terms of calculations and easily extended to impedance rough surfaces. One important fact is that KA is the most used method in NDE for rough surfaces: its implementation of the rough structure echo consists in simply modifying that of the smooth surface; the "rough" KA algorithm can thus simply mesh the smooth mean surface (surface without roughness) in small surface elements and then modify the contribution of each mesh point by including some local geometrical parameters of the rough surface.
KA has the advantage of handling surfaces of roughness smaller than a constant depending on the wavelength and on the correlation length with a relatively satisfying validity except at the vicinity of grazing observation and critical angles. Indeed we have connected and merged recent studies of the KA validity and conclude that for many applications using rather near specular observation and nongrazing incident angles, two criteria can be used for KA validity: σ 2 /L ≤ Cλp (C and λp being a constant value and the wavelength) and kL >1 or 2. If we do not want to limit the KA validity to near specular observation, the previous criteria are likely to be a little more restrictive to reach the KA validity for incidence and scattering angles over the range from −80 • to 80 • (kL > π is for instance needed) It is well admitted that KA validity for application to rough surfaces decreases for grazing angles and the region of inaccuracy in terms of scattering angles is much more important for a fluid/solid interface. One has to pay attention at the KA validity in the vicinity of particular angles (Rayleigh and critical ones). In NDE inspections, these angles are often avoided by the inspector; nevertheless it may be possible to meet them in some geometrical configurations. Improving KA or analytical approximations to better take into account head waves occurring at critical angles is an important challenge already for planar surfaces [128]. Modeling head waves on irregular surfaces is even more complicated and seems intractable using analytical methods [90,129]. Finite elements methods [130] can be useful to model both the head waves generated at critical angles, the bulk waves radiated near the surface and their possible interferences [90,129].
For planar surfaces, some smoothing modifications of KA can be proposed [69] in the vicinity of critical angles for better validity, but they do not simulate head waves. A modification of KA using an alternative integral formula for the Kirchhoff integral has just been developed [89], and it could be extended to scattering from a rough surface. Nevertheless, this perspective has to be tested and validated, since such a model will generate head waves as if they have been created on successive local non-finite tangential plane.
PT and KA are usually seen as complementary, PT being valid for small roughness scales and Kirchhoff for large scales. Several other approximations have been developed in the view of handling the inconveniences of the classical methods (PT and KA) and enlarge the simulation validity. Some of them (phase perturbation theory, two-scales models) found some applications in electromagnetism. A recent study [131] has just began in order to try to improve the two-scale models. We are not convinced that improving these models is the best perspective because they inherently own a certain complexity and there are probably better candidates (evoked in the following) in simulation methods.
The small slope approximation (SSA) is a model based on an integral field equation; it is more efficient than Kirchhoff for grazing angles but is only valid for small surface slopes. Indeed, the SSA is of great importance in ocean applications for modeling scattering from ocean surfaces on which the incidence is generally near from grazing: in such a configuration the SSA accuracy is well established. One main future challenge is to extend this approximation to layered fluid seafloors. In a study published this year [116], different approximations are proposed but they give different results in the case of strong internal reflections and also differ from Kirchhoff in specular direction. Numerical validations are clearly needed to try to find the best approximation among them; another problem is that only one extended from electromagnetism [118] can deal with multiple layers.
The parabolic equation describes all physical effects including shadowing, but this method needs a complex and heavy numerical technique to solve the parabolic equation. The parabolic equation method has recently been extended to forward scattering, but numerous efforts should be necessary in the future to build a more robust model applicable for both forward and backward 3D scattering. Such a strategy can be seen complex when compared to other methods as fully numerical ones which are generic for all scattering directions.
There is also recent progress in methods based on the integral equation; for instance the new development described in [132] leads to impressive results, notably for grazing angles; but the calculation time is still for the moment a limitation by comparison with simpler methods, if we are rather interested in near normal incidence.
Numerical methods are obviously the main rivals of the existing analytical models. For numerical techniques, obtaining 3D simulations with a reasonable computation time remains a challenge. In a very recent study, the boundary element method has been adapted [133] to reduce the 3D calculation time; the proposed adaptation is validated by using approximated models (KA and PT) showing that the classical approximations are also useful for validation purposes.
In the field of NDE, one major future challenge is to model scattering from roughness in components fabricated by additive manufacturing (AM). The fatigue life of AM specimens is very influenced by their surface topography [134]. AM processes for metallic parts are notably the selective laser melting (SLM), selective electron beam melting (SEBM), laser metal deposition (LMD) and also wire arc additive manufacturing (WAAM [135]). The corresponding roughness value of each process is discussed in [136,137]: it is 4-11 µm Ra for SLM, 15-35 µm Ra for SEBM, 10 µm Ra and 200 µm Ra for LMD, and can be very important for WAAM. Note that models are searched for the prediction of surface roughness (see a list of studies in the introduction of [138]), the optimization of process parameters as well as the roughness minimization (as in [138]) for manufacturing using WAAM and milling.
In AM processes there a lot of phenomena causing roughness [139,140], notably the sticking of nonmelted or partially melted particles on the surfaces and the formation of menisci. The latter phenomenon is more pronounced for the WAAM process (see Figure 12, below).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 19 of 27 Numerical methods are obviously the main rivals of the existing analytical models. For numerical techniques, obtaining 3D simulations with a reasonable computation time remains a challenge. In a very recent study, the boundary element method has been adapted [133] to reduce the 3D calculation time; the proposed adaptation is validated by using approximated models (KA and PT) showing that the classical approximations are also useful for validation purposes.
In the field of NDE, one major future challenge is to model scattering from roughness in components fabricated by additive manufacturing (AM). The fatigue life of AM specimens is very influenced by their surface topography [134]. AM processes for metallic parts are notably the selective laser melting (SLM), selective electron beam melting (SEBM), laser metal deposition (LMD) and also wire arc additive manufacturing (WAAM [135]). The corresponding roughness value of each process is discussed in [136,137]: it is 4-11 µ m Ra for SLM, 15-35 µ m Ra for SEBM, 10 µ m Ra and 200 µ m Ra for LMD, and can be very important for WAAM. Note that models are searched for the prediction of surface roughness (see a list of studies in the introduction of [138]), the optimization of process parameters as well as the roughness minimization (as in [138]) for manufacturing using WAAM and milling.
In AM processes there a lot of phenomena causing roughness [139,140], notably the sticking of nonmelted or partially melted particles on the surfaces and the formation of menisci. The latter phenomenon is more pronounced for the WAAM process (see Figure 12, below). The balling effect ( Figure 13) is a detrimental phenomenon that may result in breaking up the liquid scan track during SLM and produce particles in spherical shapes. The balling effect ( Figure 13) is a detrimental phenomenon that may result in breaking up the liquid scan track during SLM and produce particles in spherical shapes. WAAM). Reproduced from [141] (Open access, no copyright). he balling effect ( Figure 13) is a detrimental phenomenon that may result in breaking up scan track during SLM and produce particles in spherical shapes. igure 13. Surface topography of a specimen manufactured by selective laser melting (SLM). eproduced from [137] (Open access, no copyright).
he rough surface of AM parts can thus exhibit a kind of periodicity and elementary egg-sh larities [142]. Finite element models [143] are likely to be the best solution to both mode macrostructure of such components and their surface roughness-such models have tly used to model scattering in virtual simulated microstructures of polycrystalline mate with equiaxed grains [144][145][146] or transverse isotropic as in austenitic welds [147]. sing power brought by the spectral finite elements (SEM) method is noticeable [148-15 Figure 13. Surface topography of a specimen manufactured by selective laser melting (SLM). Reproduced from [137] (Open access, no copyright).
The rough surface of AM parts can thus exhibit a kind of periodicity and elementary egg-shaped irregularities [142]. Finite element models [143] are likely to be the best solution to both model the bulk macrostructure of such components and their surface roughness-such models have been recently used to model scattering in virtual simulated microstructures of polycrystalline materials either with equiaxed grains [144][145][146] or transverse isotropic as in austenitic welds [147]. The increasing power brought by the spectral finite elements (SEM) method is noticeable [148][149][150] in terms of calculation time and could be very useful for roughness with complex shape irregularities; SEM has recently been used for modeling ultrasonic propagation in granular materials, such as concrete [151].
In another domain, a FEM model has just been employed to model the ultrasonic response of bone-dental implant interfaces [5,152] whose irregularities can generate multiple scattering. Such a numerical method also presents the advantage of being easily adaptable in the future, for instance to treat non-linear effects at the contact interface. All the previous examples show the increase of applications for numerical methods.