## Abstract

The Kirchhoff approximation (KA) for elastic wave scattering from two-dimensional (2D) and three-dimensional (3D) rough surfaces is critically examined using finite-element (FE) simulations capable of extracting highly accurate data while retaining a fine-scale rough surface. The FE approach efficiently couples a time domain FE solver with a boundary integration method to compute the scattered signals from specific realizations of rough surfaces. Multiple random rough surfaces whose profiles have Gaussian statistics are studied by both Kirchhoff and FE models and the results are compared; Monte Carlo simulations are used to assess the comparison statistically. The comparison focuses on the averaged peak amplitude of the scattered signals, as it is an important characteristic measured in experiments. Comparisons, in both two dimensions and three dimensions, determine the accuracy of Kirchhoff theory in terms of an empirically estimated parameter *σ*^{2}/λ_{0} (*σ* is the RMS value, and λ_{0} is the correlation length, of the roughness), being considered accurate when this is less than some upper bound *c*, (*σ*^{2}/λ_{0}<*c*). The incidence and scattering angles also play important roles in the validity of the Kirchhoff theory and it is found that for modest incidence angles of less than 30°, the accuracy of the KA is improved even when *σ*^{2}/λ_{0}>*c*. In addition, the evaluation results are compared using 3D isotropic rough surfaces and 2D surfaces with the same surface parameters.

## 1. Introduction

The Kirchhoff approximation (KA) is the classical methodology often used to calculate the scattering behaviour from rough surfaces. It is widely used in many fields, such as underwater acoustics [1], radar detection [2] and seismology [3,4], to produce fast predictions of statistics for the scattering of waves via a Monte Carlo approach. For elastic media, there is strong interest in non-destructive evaluation (NDE) and ultrasonics in applying the Kirchhoff model to understand the scattering behaviour from rough defects as this is of critical importance in assessing the structural integrity of, for instance, components in the nuclear industry. In that area early work by Ogilvy [5] studied the effects of roughness on scattered wave amplitude-based and time-based sizing techniques, and concluded that careful interpretation of measured data is required when roughness is present. More recently, Zhang *et al.* [6] have investigated the effect of roughness and inclined angles on sizing the defect using the ultrasonic array imaging with the total focusing method. Given the wide ranging use of Kirchhoff theory in NDE, it is natural to want to determine, or clarify, the range of validity of KA.

It is well known, as pointed out by Ogilvy [7], that Kirchhoff theory is a high-frequency approximation and its accuracy depends on many factors (surface roughness, grazing incidence/scattering angles, etc.). It is therefore important to find the range of validity, and when the use of KA should be expected to be reliable, particularly as it is used for safety critical inspection. Historically, there have been several attempts to evaluate the accuracy of KA, and Ogilvy provided a detailed summary in her book [7]. For instance, in [8], Ogilvy implemented a variational principle to quantify the accuracy of the scalar KA. Another evaluation was provided by Thorsos [1] comparing the boundary integration method and acoustic Kirchhoff theory using random Gaussian surfaces. The conclusion was that the correlation length λ_{0} is the most critical parameter and the RMS *σ* also affects the accuracy of the Kirchhoff model. A similar study was extended to a Pierson–Moskowitz sea surface [9] and the KA was shown to accurately predict the incoherent scattering near the specular direction for incident grazing angles as low as 10°; this is typical of the impressive accuracy and utility that has led to the widespread use of Kirchhoff theory.

The above examples are all for scalar problems, and the validity of the vector wave KA (occurring in electromagnetics and elastodynamics) differs from the scalar case. The difference is because restrictions on surface properties that arise from effects such as shadowing and multiple scattering can be made more severe by mode conversion, into bulk waves and surface, or interfacial, waves. Chen & Fung [10] investigated electromagnetic wave scattering from rough surfaces by the KA and the moment method finding that electromagnetic Kirchhoff theory provides good agreement over small angles of incidence when *kσ* is less than 0.2 and *kλ*_{0} is less than or equal to 2.0. Here, *k* is the incident wavenumber. Robertsson *et al.* [3] compared the elastic KA with the finite difference method and the spectral element method with an application in reflection seismology. In this work, some discrepancies between the KA and the other two methods were shown particularly in terms of amplitudes. Elastic Kirchhoff theory was also investigated in connection with the application in NDE by Roberts [11], who compared it with the boundary element method. More recently, Zhang *et al.* [12] calculated the scattering matrix from Gaussian rough defects in comparison with the finite-element (FE) method and concluded that the elastic KA is valid with the roughness *σ*≤0.3λ and λ_{0}≥0.5λ at incidence and scattering angles over the range from −80° to 80°.

However, to the best knowledge of the authors, all of these previous research papers regarding the range of validity of the Kirchhoff theory, especially in the elastodynamic regime, are only based on two-dimensional (2D) rough surfaces or at most 2.5D (corrugated) surfaces due to limitations in the simulation methods used. We are now in a position to sidestep these limitations and perform highly accurate numerical simulations in three dimensions due to recent advances such as hybrid methods [13,14]. We have been motivated to explore the validity of Kirchhoff theory in three dimensions as researchers developing modelling tools using three-dimensional (3D) KA for NDE in the power industry have observed significant differences of several dB between experimentally measured reflections from rough defects and KA simulations of the same cases [15]. These discrepancies arise because surface scattering in a real inspection is an inherently complex process, involving phenomena such as mode conversions, surface waves and shadowing, which are not at all included in the KA. These phenomena are known to cause the greatest errors when the angle of incidence and the height of the roughness are large, and are also expected to be different for two dimensions and three dimensions. Therefore, to gain a comprehensive understanding of the performance of the 3D Kirchhoff theory on random rough surfaces, with height variations in both *x*- and *y*-directions, it is necessary to actually simulate the full 3D problems. Motivated by the experimental observations, and the apparent discrepancies with the KA in elasticity, we carefully examine the range of validity of both the 3D and 2D elastic KA by comparison to a reference numerical method. Using an empirical procedure, the ratio *σ*^{2}/λ_{0} is found to be important in determining the region of validity of the KA, and this ratio has different upper bounds in two dimensions and three dimensions. In addition, the incident angle also affects the performance of the KA.

This paper is organized as follows: §2 introduces the numerical models, including the rough surface profile, the Kirchhoff model and the FE boundary integration method used as a reference for comparison. Section 3 describes details of the simulations in both two dimensions and three dimensions. Error analysis in two dimensions and three dimensions is discussed in §4, and the concluding remarks are made in §5.

## 2. Numerical models

We intend to simulate, as closely as possible, the scattering of a plane wave from an infinite rough surface and the scattering geometry (in 2D) is depicted in figure 1. The incident and scattering angles are denoted as *θ*_{i} and *θ*_{s}, respectively, and a stress-free boundary condition is assumed along the rough surface.

### (a) Gaussian random rough surface

The surface of the defect is represented as ‘rough’ in the sense that the surface height data obey a certain statistical distribution, and these realistic surface profiles are generated using some statistical parameters. Gaussian random surfaces were commonly used historically [1,16], and also some fatigue cracks are found experimentally to have a Gaussian spectrum [12]. Therefore, we use Gaussian surfaces herein to find the region of validity of KA. Figure 2 shows a typical 3D isotropic (the roughness varies in the same manner in both *x*- and *y*-directions) Gaussian rough surface profile, with the mean value of height data *h*(*x*,*y*) being zero. The probability density function *p*(*h*) is used to describe the statistical distribution of the height data, which is a Gaussian function in this study
*σ* is the RMS of height, and the empirical formula to estimate *σ* is then given by
_{x} and λ_{y} are called the correlation lengths in the *x*- and *y*-directions, as distances over which the correlation function falls by 1/*e*; we assume an isotropic rough surface hereafter, implying λ_{x}=λ_{y}=λ_{0}. In this isotropic case, the correlation function simplifies to
*R*^{2}=*x*^{2}+*y*^{2}.

The commonly used moving average method [16] is used to generate the rough surface by controlling the RMS value *σ* and the correlation length λ_{0}. They are in units relative to the compressional wavelength λ_{p}. This method produces a series of uncorrelated random numbers that are convolved with a set of weighting functions to calculate the required surface height data *h*(*x*,*y*). Figure 2*a*,*b* shows one realization of a 3D Gaussian rough surface and the corresponding height distribution function. The effect of changing *σ* and λ_{0} on the roughness is clearly seen in figure 2*c*,*d*, compared with the original surface profile in figure 2*a*. It should be noted that each generated realization of the surface profile is not identical given the random nature of the surface. However, for a specific *σ* and λ_{0}, all realizations should follow the same statistical parameters and can be classified into the same category. In later sections, Monte Carlo simulations are run for different categories of rough surfaces, and a general rule with respect to the statistical parameters regarding the validity of the KA is obtained.

### (b) Finite-element boundary integration approach

The finite-element method is a very general numerical tool that is highly effective for elastic wave problems, and we use it as a reference against which to evaluate the accuracy of the Kirchhoff theory. However, full 3D FE modelling is very computationally expensive and even the 3D free meshing of the irregular surface depends on the irregularity of the defect shape. To overcome these issues, that have limited previous studies to 2D or 2.5D simulations, we use a state-of-the-art FE boundary integration method implemented in a GPU-driven software package Pogo [17], with a mixed meshing algorithm in three dimensions. This package is sophisticated in its use of computer hardware, algorithms and in using hybrid methods. Efficient and accurate Monte Carlo simulations using the FE method in three dimensions are performed; for brevity, the full details of the methodology are not included here as they are given in [13].

The FE model to calculate the scattering from a rough surface when internal waves are incident, is shown in figure 3. A rough backwall is surrounded by smooth surfaces extended into an absorbing region [18]. The smooth transition between the rough part and the smooth part is achieved by multiplying a spatial Hanning window to the whole bottom surface. A plane P wave is excited by forcing along a line located just above the rough backwall. Specific time traces of forces, which are calculated from a finite difference time domain method using a modified forcing approach [19], are fed into each excitation node. This approach was initially designed to calculate forcing signals to excite the correct incident waves inside a rectangular box. Here, we simply use the code to obtain the forces on one side of the box, denoted as the excitation line. It should be noted here that the two tips of the excitation line need to be buried into the absorbing region to prevent any unwanted circular-crested waves, which would otherwise be generated from the two ends.

After executing the FE computation in the local box, the boundary displacements recorded at the rough surface are used in the Helmholtz integration with a stress-free boundary condition [20] to calculate the scattered waves:
*Σ*_{ijk} is the stress Green's tensor, *u*_{i} is the *i*th component of the total displacement, *n*_{j} is the *j*th component of the unit normal vector surface pointing towards the observation point at ** R** and

**represents the position of the surface point.**

*r*With a far-field approximation, equation (2.5) can be simplified as
** D**=

**−**

*R***. Before using the boundary displacements in equation (2.6), the Hanning window function is applied to them, to reduce the ‘edge effects’ in the approximation of an infinite surface. In two dimensions, a popular way to approximate plane wave scattering from a truncation to an infinite surface is to use a tapered plane wave [21]. In this approximation, the half beam width must be much smaller than the surface length, but roughly larger than 3λ to eliminate the tip effects. Unfortunately in three dimensions, it is difficult to satisfy both criteria simultaneously due to the limited size of the 3D surface that can be simulated. Our approach is to apply a Hanning window to smooth the boundary displacement and, although this does not completely cancel the circular-crested waves from the edges, it does mitigate the effect substantially.**

*r*In practice, equation (2.6) is implemented in the frequency domain and the results are synthesized back to the time domain to obtain the scattering signals using the inverse fast Fourier transform. The model is extended to three dimensions by replacing the FE rectangular box with a cubic box, the 2D surface with a 3D surface and a source line with a source plane, as shown in figure 3*b*. Note that, for brevity, the absorbing region in three dimensions is not drawn here.

### (c) Elastic Kirchhoff model

The Kirchhoff model with a rough surface in the *x*–*y* plane is illustrated in figure 4. The KA assumes that the motion of a point on the surface is the same as if the point were part of an infinite plane insonified by the incident wave, as depicted in the zoom-in of figure 4. The total displacement at this point is calculated as a summation of the incident wave and the reflected compressional(P) and shear(S) waves:

Here, *u*_{t} represents the total displacement, *u*_{0} is the amplitude of the incident P wave, *A*_{pp} and *A*_{ps} are reflection coefficients of P and S waves, respectively, and *d*_{0}, *d*_{p} and *d*_{s} are the displacement polarization vectors for the incident P and reflected P/S waves. The same Hanning window as implemented in the FE method is also applied to smooth the total displacement from KA to guarantee a fair comparison. By substituting the total displacement into equation (2.6), the scattering displacement is computed explicitly in the far field.

## 3. Simulations in two- and three dimensions

In this paper, although we mainly investigate the compressional wave incidence/scattering case, the method can also equally be applied to study shear wave scattering. For the simulations shown, we use the material in the bulk medium to be aluminium with Young's modulus of 70 GPa, density of 2700 *kg* *m*^{−3} and Poisson's ratio of 0.33. Thus, the compressional wave speed is 6198 *m* *s*^{−1}, and the shear wave speed is 3122 *m* *s*^{−1}. The surfaces used follow the Gaussian distribution, which is generated by the moving average method discussed in §2. In two dimensions, 15 different rough surface profiles are tested with RMS *σ*=λ_{p}/8, λ_{p}/6, λ_{p}/5, λ_{p}/4 and λ_{p}/3, and correlation length λ_{0}=λ_{p}/2, λ_{p}/3 and λ_{p}/4. While in three dimensions, 10 rough surface profiles are used with the same RMS values as those of the 2D cases, and with correlation lengths λ_{0}=λ_{p}/2 and λ_{p}/3. In the 3D case, the example with λ_{0}=λ_{p}/4 is not used because this is too demanding for the 3D free-meshing algorithm which tends to fail when such a small correlation length is tested, and it should be noted that despite the efficiencies of the methodology, these calculations remain extremely demanding in terms of computer resources. Numerous simulations are run with incidence angles *θ*_{i} ranging from 0° to 30° to calculate scattering signals with *θ*_{s} from −90° to 90°. The distance from the surface to the far-field projection point is 50 mm (≈32λ_{p}). For each roughness, 50 realizations of surfaces are generated and Monte Carlo simulations are performed to obtain a statistically meaningful result for comparison between the KA and the FE method. The number of realizations is chosen by considering the conflicting requirements of both the statistical stability and the computational cost, especially in three dimensions. In NDT inspection scenarios, we use the amplitude of the scattering signal because the amplitude directly determines the probability of detection [22]. The statistical parameter used for comparison is therefore the ensemble average of the peak of the Hilbert transformed scattered signals.

### (a) Simulations using two-dimensional rough surfaces

Simulations are first run using 2D surfaces and the corresponding FE model is meshed by linear triangular elements (equivalent to CPE3 in Abaqus (Dassault Systemes Simulia Corp., Providence, RI)) with an element size of λ_{p}/30. From earlier studies [23], this element size is sufficient for the convergence requirement in this study. With a proper partition, only the region surrounding the rough backwall is free meshed and the remaining region can be regularly meshed. A five-cycle Hanning windowed compressional wave with centre frequency of 4 MHz is used and the length of the 2D rough surface is about 5.2λ_{p}. The dimension of the FE domain including the absorbing region is 14.8×9.4 *mm*^{2} (≈9.6λ_{p}×6.1λ_{p}), and the thickness of the absorbing layers is 2.4 mm (≈1.5λ_{p}).

For each scattering angle, the compressional displacement is decomposed from the total displacement by using *θ*_{i}=*θ*_{s}=0°, with increased RMS values and a fixed correlation length. The amplitude is shown on a linear scale, and it is divided by the peak of the normal pulse echo response from a smooth surface of the same size. From figure 5, with a positive increment of the roughness, the amplitude of the scattering waves decreases significantly. In addition, the scattering waveforms become more complicated and a clear second arrival can be seen around 17 *μs*. This second arrival is the mode converted shear wave because it matches the correct arrival time, which equals to 4 *mm*/(6198 *m* *s*^{−1})+50 *mm*/(3122 *m* *s*^{−1})=16.7 *μs*. As shown in figure 5*a*–*d*, the Hilbert transformed scattering signals calculated from KA agree well with those from the FE method, when *σ* is smaller than λ_{p}/4. Such good agreement can no longer be seen in figure 5*e*, as the roughness increases to *σ*=λ_{p}/3, indicating that KA then becomes inaccurate.

Taking the mean value of the peak of these time traces from all realizations, the error between the KA and the FE across all scattering angles can be more accurately demonstrated by comparing the scattering patterns as in figure 6. For comparison, the scattering pattern from a smooth surface with the same dimension is also shown. Only one curve is seen because the KA and the FE predict the same result for a smooth surface. Note that the angular spread is caused by both the finite length of the crack and the presence of the roughness. The amplitude is shown on the dB scale, and the reference is again the peak of the normal pulse echo response from a smooth surface. If more realizations are run, we would expect the graph to become symmetric about *θ*_{s}=0°. If the tolerance of error is set to be 1 dB, as commonly used in NDE, the range of the scattering angle *θ*_{s} for which KA is acceptable is from −70° to 70° when *σ*≤λ_{p}/4. However, if the RMS *σ* increases to λ_{p}/3 then, even at the specular direction the KA no longer matches with the FE. In addition, at near specular directions the KA underestimates the amplitude but at near grazing angles it overestimates the amplitude compared with the FE results; this is because multiple scattering phenomena are neglected in the KA. The consequences of neglecting the multiple scattering is that less energy is reflected back to the normal, or near normal, scattering angles and more energy than expected is distributed at near grazing angles. Figure 6*f* shows that both the mean absolute value and the standard deviation of the error at the normal backscattering direction increase as *σ* increases. The mean absolute error and the standard deviation are calculated using the following equations:
*E*_{n} is the error for one realization and *N* is the number of total realizations. *A*^{KA}_{n} and *A*^{FE}_{n} are the scattering amplitudes calculated from KA and FE for one realization, respectively. *std*(*E*) is the standard deviation of the error.

### (b) Simulations using three-dimensional rough surfaces

The 2D simulation results shown above are equivalent to 3D models using corrugated surfaces which have a height variation only in the *x*-direction. Such simplifications have been performed in the past to ease computational demands [1,8,12], but in reality all defects have roughness in two dimensions. Although the range of validity for the 3D KA might, intuitively, not differ from the 2D KA much, some discrepancy between the two is expected given we now have roughness in one more direction.

The FE model as discussed and implemented with a 2D surface is extended to three dimensions with minor modifications. The 3D CAD software Rhino (Robert McNeel & Associates, Seattle, WA, USA) was implemented to build the rough surfaces as shown for example in figure 7. The dimension of the rough surface is 5×5 *mm*^{2}, that is approximately 3.3λ_{p}×3.3λ_{p}. The surface is created by spline interpolating the surface point cloud generated using the moving average method. Restrictions are included in the interpolation algorithm so that the CAD surface is exactly on the positions of the generated points. This is crucial to guarantee that the produced surface does not deviate from the required point position, and consequently to avoid any smoothing effect caused by general surface interpolation algorithms.

Free meshing algorithms using linear tetrahedral elements (C3D4 in Abaqus) are commonly used for the 3D FE modelling of scattering from obstacles with irregular geometries [24,25]. As is well known, free meshing tends to randomize the distribution and the shape of each element, and thus may introduce unwanted mesh-scattering and dispersion [26]. We therefore combine two different meshing algorithms in three dimensions with a code developed using Matlab (MathWorks, Natick, MA, USA). The region just above the 3D rough surface is meshed via a free meshing algorithm to produce elements around the surfaces. This very local meshing profile is then used as an input to the code to grow a regular mesh to fill up the remaining region of the 3D model. The regular meshed region is meshed with many hexahedral cells and each cell is composed of six linear tetrahedral elements as shown in figure 8*a*. The dimension of the 3D FE domain including the absorbing region is 12.8×12.8×4.4 *mm*^{3} (≈8.3λ_{p}×8.3λ_{p}×2.9λ_{p}), and the thickness of the absorbing layers is 2.4 mm (≈1.5λ).

Figure 8*b* shows a local view of the mesh profile of one 3D FE model. It should be noted that at the boundary of the free meshing and the regular meshing regions, the two neighbouring elements need to have the same hypotenuse to prevent any spurious reflection. In this manner, the mesh minimizes the complexity caused by the free meshing algorithm while still capturing the exact shape of the complex rough surface.

In a similar manner to that deployed in the 2D studies, a Monte Carlo simulation is performed, with various roughness and incidence/scattering angles, and the averaged peak amplitudes from 50 realizations are compared between the KA and the FE model. Figure 9 shows the comparison for the case of a normally incident wave. Good agreement can be seen from −70° to 70° when *σ*=λ_{p}/8, λ_{p}/6, λ_{p}/5 and λ_{p}/4. The scattering pattern from a smooth surface with the same dimension is also shown for comparison. However, when *σ* increases to λ_{p}/3, the agreement no longer exists even in the normal backscattering direction. Note that the errors between the KA and the FE are relatively higher than the corresponding 2D cases in figure 6. Again the mean absolute error and the corresponding standard deviation for each *σ* are plotted in figure 9*f*.

## 4. Error analysis in two- and three dimensions

### (a) The effect of surface roughness

As mentioned by Ogilvy [7], historically the most frequently cited validity criterion from a geometrical point of view for the KA is *r*_{c} is the local radius of curvature and *θ*_{i} is the global angle of incidence. This equation explains the fundamental restriction of the rate of change of the surface height gradient, which determines the validity of the tangential plane assumption used as a starting point for the KA. However, other related works [1,10] indicate that this commonly used criterion has not proved to be practically useful. The errors from the KA are not only from the ‘local’ geometrical tangential plane assumption, but also arise from ‘global’ effects, such as multiple scattering and shadowing effects due to the overall shape of the rough surface. These ‘global’ errors become larger as the roughness or the frequency increases, and the corresponding ‘global’ effects are not easy to quantify directly. Thorsos [1] suggests the use of λ_{0} as a crucial surface parameter including both the local and the global effects. In this section, we propose a different surface parameter for elastic materials, a function of both λ_{0} and *σ* to give a rough idea when the use of the KA is valid for practical applications.

Figure 10 shows the errors of the averaged peak amplitude between the 2D KA and the FE with respect to the RMS *σ*, with a normal incident/backscattering configuration. Different curves represent the numerical errors using surfaces with different correlation lengths. The tolerance of error is set to the value commonly used in NDE applications of 1 dB, as shown by the flat dashed line, indicating that the roughness is acceptable for the use of KA once the corresponding error is below this threshold. The maximum values of *σ* for a given correlation length λ_{0} and error tolerance can hence be found by observing the trend of the error curves. These values are summarized in table 1.

Clearly, as can be seen from figure 10, increasing the RMS *σ* and decreasing the correlation length λ_{0} result in a larger value of the error. This increase in error is because the surface effectively becomes more rough when changing the two surface parameters in this manner. Therefore, we propose a simple criterion as a function of both *σ* and λ_{0}, based on the inversely proportional relationship
*a* is the weighting factor for *σ* and *c* is an unknown constant representing an upper bound for this inequality. Note that this form is not unique, for example one may also assume a weighting factor for λ_{0} instead of *σ*. However, in this paper, we propose equation (4.1) and it will be shown later that this formulation can approximately estimate the region of validity of the KA in both two dimensions and three dimensions.

The acceptable values of *σ* and λ_{0} can be explicitly calculated from this formula once *a* and *c* are known. An alternative form converted from equation (4.1) is used here to estimate these two unknown coefficients.

Substituting the three values of *σ*_{max2D} from table 1 into equation (4.2) and taking the log of both sides yields the following form of matrix multiplication:

By multiplying the term on the right-hand side with the pseudo inverse of the first term on the left side, the weighting factor *a* and the upper bound *c* are calculated in the sense of least squares. Recalling that we are using the tolerance of 1 dB, typical of NDE, from the calculation we see that for this error *a*≈2. The value of *a* could change if a more or less stringent tolerance of error were required. The observed criterion for the validity of KA can be expressed as
*σ*_{max} with a given λ_{0} or the λ_{0min} with a given *σ*

A numerical case is performed here to test the proposed criterion. As can be seen from figure 10, when *σ*=λ_{p}/3 the errors for all curves are above the 1 dB threshold, which implies that the acceptable λ_{0} needs to be larger than λ_{p}/2. By substituting *σ*=λ_{p}/3 into equation (4.5), λ_{0min} is estimated to be roughly 0.6λ_{p}. A Monte Carlo simulation with 50 surfaces is run for this specific *σ* and λ_{0}, as was performed with the other roughness tested. As can be seen from figure 11, the error between the KA and the FE method when *θ*_{i}=*θ*_{s}=0° is 0.86 dB. This error is just below the 1dB threshold, implying that when *σ*=λ_{p}/3, λ_{0min}=0.6λ_{p}, which is the same as the value predicted from equation (4.5). The simple empirical criterion proposed from observation can therefore be applied to roughly estimate the region of validity of KA.

### (b) The effect of scattering/incidence angle

The Kirchhoff theory is known to be insufficient to model the scattering of waves at large scattering angles [7]. For example, as seen from figure 6*a*–*d*, the scattering amplitude calculated by the KA is several dB higher than the FE when the scattering angle *θ*_{s}>70°. The explanation given in §3 is that the KA cannot account for multiple scattering effects.

In the same way, large incidence angles also result in considerable errors of the Kirchhoff model. In this study, we alternatively focus on the effect of modest incidence angles on the accuracy of the KA. Specifically, we compared the averaged scattering amplitude of the KA and the FE at incidence angles from 0° to 30° with an interval of 5°. Figure 12*a* shows the comparison results of the amplitude with a slightly oblique incidence angle (*θ*_{i}=30°) when *σ*=λ_{p}/3, λ_{0}=λ_{p}/2. Excellent agreement is found with the scattering angle *θ*_{s} ranging from −65° to 65°. It is different from the case with the normal incidence angle as shown in figure 6*e*, where the error is above 1 dB. It indicates that even with a very rough crack, an acceptable inspection result may still be achieved using the Kirchhoff model with a small oblique incidence angle. Figure 12*b* further illustrates this point as the error of the scattering signal at the specular direction is plotted as a function of the incidence angle. As can be clearly seen, the error decays as the incidence angle increases for the curve representing *σ*=λ_{p}/3; for the other curves when *σ*<λ_{p}/3, the error roughly remains the same, well within 1 dB. The results for surfaces with a correlation length of λ_{p}/2 are shown here but the trend is typical for all other correlation lengths tested in this paper.

The ‘Rayleigh parameter’ [7] which has been used to classify the rough surface and the smooth surface can be quoted here to give a qualitative physical explanation. The Rayleigh parameter is
*R*_{a}, the more destructively the two reflected waves interfere with each other. Therefore, it is one feature which is proportional to the roughness of the surface. According to equation (4.6), if the Rayleigh parameter *R*_{a} is fixed then the RMS *σ* is inversely a measure of *R*_{a}, with a modest incidence angle the maximum value of *σ* can be relatively larger than that when *θ*_{i}=0°. It should be noted that the effect of a modest incident angle is qualitatively analysed here by the use of the Rayleigh parameter, but it cannot be applied in the same manner when the incidence angle is large.

### (c) The effect of dimension (two or three dimensions)

Similar criteria for the 3D KA can be proposed in the same manner as was deployed for two dimensions. Figure 13 shows the error between the KA and the FE method with respect to *σ* for two different values of λ_{0}. Similar trends can be observed in the two curves but the errors are relatively higher, compared with those in figure 10 for the 2D KA. The same strategy is implemented to estimate the unknown parameters in equation (4.2), but using *σ*_{max3D} from table 1. By using the least-squares method, the best-fitted weighting factor remains at *a*≈2, and the upper bound reduces to *c*=0.14λ_{p}. The criterion in three dimensions can therefore be expressed as
_{p} is smaller than the corresponding value 0.2λ_{p} in two dimensions, which highlights one important discrepancy of the validity between the KA in two dimensions and three dimensions. Specifically, as we can see, the values of *σ*_{max} for a given λ_{0} in three dimensions are smaller than those in two dimensions as shown in table 1; this implies that the criterion for the 3D KA is stricter than that in two dimensions. We can consider this by examining the surface locally, where we know that the magnitude of the gradient of one surface point in two dimensions and three dimensions can be expressed as

Obviously from the above equation, |∇*h*(*x*,*y*)|_{3D} is larger than |∇*h*(*x*)|_{2D} if it is assumed that |∂*h*/∂*x*|_{2D}=|∂*h*/∂*x*|_{3D}, which is due to adding an additional term |∂*h*/∂*y*|_{3D}. The conclusion is that the rate of change of height is faster for the 3D surface which consequently appears rougher than the corresponding 2D surface.

Apart from this, height variations in the *y*-direction can contribute to out of plane scattering effects besides the scattering waves inside the *x*–*z* plane. Globally more shadowing effects and multiple scattering would occur, due to the roughness in the extra dimension. It would result in extra errors of the 3D KA in addition to the local errors from the tangential plane assumption. As a result, restrictions on the validity of 3D KA are made more severe than the 2D KA due to this one more dimension of roughness.

If the incidence angle is slightly oblique at 30°, as in figure 14*a*, then the FE and the KA models show agreement, but only within a very narrow angular range around the specular direction, roughly from 30° to 50°. By contrast, the acceptable region of the scattering angles in dimension with the same roughness is much larger, ranging from −65° to 65°, as shown in figure 12*a*. This rapidly reduced angular range in three dimensions when *σ*=λ_{p}/3 is another important feature between the region of validity of 3D and 2D KA. It may be because in three dimensions, there are more multiple scattering and mode conversions than two dimensions, which are not modelled using the KA. The scattering waves contributed from these phenomena would spread to all scattering angles, reducing the accuracy of the 3D KA at non-specular directions. On the other hand, the effect of a modest incidence angle on the 3D KA is demonstrated in figure 14*b*. All the curves representing different *σ* show a decay of the error when increasing the incidence angle *θ*_{i}. The errors are relatively higher than those shown in figure 12*b* for the 2D cases.

## 5. Conclusion

A numerical study to determine the range of validity of the KA as used to calculate the elastic wave scattering signals from Gaussian rough surfaces has been performed. The Kirchhoff theory is evaluated for comparison with a FE method using 2D and 3D isotropic surfaces with different roughness characterized by RMS *σ* and correlation length λ_{0}. Monte Carlo simulations of multiple realizations are run with a variety of incidence/scattering angles, and the averaged peak amplitude of the scattering signals are used for comparison. It is found that, with a normal incidence angle *θ*_{i}=0°, the KA is valid when *σ*^{2}/λ_{0}≤*c*, with −70°<*θ*_{s}<70°. In addition, a modest incidence angle within 30° can improve the accuracy of the KA when *σ*^{2}/λ_{0} exceeds the constant upper bound *c*. The above criteria are derived empirically for an estimation of the region of validity of the KA.

We have particularly examined the difference of the valid region with 2D and 3D surfaces and have found that the criterion for the 3D KA is stricter. First of all when *θ*_{i}=0°, the upper bound *c* of *σ*^{2}/λ_{0} in two dimensions is 0.20λ_{p}, and in three dimensions it reduced to 0.14λ_{p}. In other words, the acceptable range of *σ* and λ_{0} in three dimensions is smaller than those in two dimensions, which is the most important discrepancy between the 2D and 3D KA. This is caused by the increased local RMS gradient, and the global multiple reflections and shadowing effects as well. Furthermore, similar to two dimensions, a modest incidence angle of less than 30° can also improve the accuracy of the 3D KA. However, the acceptable angular range of the scattering angle is dramatically reduced compared with the same situation in two dimensions, being only around the specular direction. In fact, these differences arise from extra surface height variations in one more dimension, which makes 3D surfaces as viewed by the incident wave appear to be rougher. In practice, the criteria for the 3D KA must be taken into consideration in a real NDE inspection because naturally all real defects are rough in three dimensions.

## Ethics statement

The paper does not include any contents regarding animals or human subjects.

## Data accessibility

Data have been submitted as electronic supplementary material.

## Funding statement

The authors are grateful to the UK Research Centre in Non Destructive Evaluation (RCNDE), the Engineering and Physical Science Research Council (EPSRC), and to Amec Foster Wheeler, Rolls-Royce Marine and EDF Energy, for funding this work.

## Authors' contributions

F.S. performed all simulations, analysed data and wrote the paper. W.C., M.J.S.L., E.A.S. and R.V.C. provided valuable discussions with the first author and critically revised the paper.

## Conflict of interests

We have no competing interests.

- Received December 18, 2014.
- Accepted April 4, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.