Rayleigh waves are well known to attenuate due to scattering when they propagate over a rough surface. Theoretical investigations have derived analytical expressions linking the attenuation coefficient to statistical surface roughness parameters, namely, the surface's root mean squared height and correlation length and the Rayleigh wave's wavenumber. In the literature, three scattering regimes have been identified—the geometric (short wavelength), stochastic (short to medium wavelength), and Rayleigh (long wavelength) regimes. This study uses a highfidelity twodimensional finite element (FE) modelling scheme to validate existing predictions and provide a unified approach to studying the problem of Rayleigh wave scattering from rough surfaces as the same model can be used to obtain attenuation values regardless of the scattering regime. In the Rayleigh and stochastic regimes, very good agreement is found between the theory and FE results both in terms of the absolute attenuation values and for asymptotic power relationships. In the geometric regime, power relationships are obtained through a combination of dimensional analysis and FE simulations. The results here also provide useful insight into verifying the threedimensional theory because the method used for its derivation is analogous.
I. INTRODUCTION
Elastic waves guided on the surface of a solid are well known to attenuate if the surface is not perfectly flat as a result of scattering.^{1} The attenuation of Rayleigh waves from rough surfaces has been described analytically, but experimental validations have proved difficult to achieve over a wide range of parameters of the roughness. One reason is that mathematical models are valid in regions that are difficult to replicate in practice because they either require very low roughness or predict high attenuation values, which would render the waves undetectable.
Early work to describe the attenuation of Rayleigh waves resulted in the derivation of expressions for waves on flat surfaces of solids whose material properties induced attenuation. For instance, Maris^{2} derived expressions for the attenuation of Rayleigh waves on a dielectric crystal with arbitrary crystallographic orientation, where temperature and viscosity were also considered. Their work followed the experimental results of Salzmann et al.,^{3} who used lasers to measure the effect of the temperature and frequency on the attenuation of Rayleigh waves propagating along quartz crystals. These early studies considered specific attenuation cases related to material properties but did not take into account the roughness that unavoidably exists on all surfaces and causes attenuation even when the waves propagate on perfect lossless elastic materials; the attenuation occurs by partial scattering of the waves from the geometric features of the roughness.^{1}
Some of the first analytical expressions for the attenuation of Rayleigh waves from rough generalised surfaces were derived by Maradudin and Mills^{4} and Urazakov and Fal'kovskii.^{5} The rough surface was described using two statistical parameters—the root mean squared (RMS) height, δ, which describes the amplitude of the peaks and troughs of the roughness, and the correlation length, Λ, which is a measure of the spacing of those peaks. Urazakov and Fal'kovskii used the Rayleigh method for their studies. This method solves the Rayleigh wave equation, which predicts the creation of Rayleigh waves on a stressfree boundary^{6} but was extended by the authors to apply that stressfree surface condition across the rough surface.
Maradudin and Mills^{4} used a Green's function approach to solve the relevant equations. Both approaches predict that for the threedimensional (3D) case, when the roughness is low (for instance, $ \delta / \Lambda < 0.3$), the attenuation coefficient is proportional to the fifth power of the wave's frequency, f, in the region where the Rayleigh wavelength, λ_{R}, is much greater than Λ. Maradudin and Mills also demonstrated that the attenuation coefficient is proportional to $ \delta 2$. Subsequently, Eguiluz and Maradudin^{7} published an updated version of their derivations in which the additional scattering of the Rayleigh waves to bulk waves was considered. The new results also demonstrated a proportionality of the attenuation coefficient to f ^{5}, as well as a new result of a proportionality to the ratio $ \delta 2 / \Lambda 2$.
Following the work by Eguiluz and Maradudin, de Billy et al.^{8} completed some experimental work to verify the theory derived in Ref. 7. In their work, attenuation measurements were taken for rough duraluminum and titanium samples, and the effect of varying λ_{R} on the attenuation measurements was investigated. The authors found good agreement between their experimental results and the f ^{5} relationship predicted by the theory at the long wavelength limit; however, they observed that this fifthorder proportionality did not hold true at smaller λ_{R} values. More recently, Kosachev and Gandurin^{9} studied the dispersion attenuation of Rayleigh waves on statistically rough hexagonal crystals to expand the work in Ref. 7. In their work, the scattering from the rough surface of an anisotropic hexagonal crystal was theoretically studied with the authors deriving an expression for the attenuation coefficient when the scattering occurs at a generalised crystal orientation. Their expression reduces to the same f ^{5} relationship derived in Ref. 7. for the isotropic case in the long wavelength region. Finally, Chukov^{10} used the RayleighBorn approximation to derive similar relationships to Ref. 7 for the isotropic case, arriving at the same power relationships. In addition, Chukov derived expressions for 2D roughness, which are discussed in more detail below.
However, the theory is restricted to specific $ \delta / \Lambda $ and λ_{R} combinations. In particular, a large proportion of the theory has been derived in the Rayleigh regime—this is a region where $ \lambda R \u226b \Lambda $, which is a very low frequency regime. Therefore, the abovementioned considerations motivate two research problems, which this paper attempts to solve—first, to create a finite element (FE) model to validate the existing analytical expressions in the Rayleigh region and, second, to extend the results by FE modelling to other more practically relevant regimes.
To obtain a representative attenuation value for a combination of δ, Λ, and f from numerical modelling, it is necessary to average over a sufficiently large number of attenuation values, which were obtained from individual surfaces characterised by those specific statistical parameters. It is also necessary to have a model sufficiently large to accommodate a representative scattering distance. The implementation of a 3D model would incur a significant computational burden and runtimes for such a wide range of parameter values and, therefore, this paper conducts a comprehensive 2D analysis. The validation in two dimensions also infers validation in three dimensions as the same theoretical approach has been implemented for the derivations in both cases.
In addition to the 3D case discussed above, Chukov demonstrated that for the 2D case, the attenuation coefficient is proportional to f ^{4}, $ \delta 2$, and Λ in the Rayleigh region in comparison with the $ f 5 \delta 2 \Lambda 2$, which has been derived for 3D roughness. Identical power relationships were also derived by Huang and Maradudin,^{11} using the same small perturbation method implemented in the 3D analysis,^{7}and the f ^{4} relationship was observed experimentally in Ref. 8. These are the proportionality relationships that this paper validates. In addition to this, Maradudin and Eguiluz^{7} and Huang and Maradudin^{11} have also derived an expression that gives quantitative values for the attenuation coefficient—the results from the FE simulations in this study are also compared against this expression.
Following this validation, this paper looks into the attenuation coefficient's behaviour in regimes outside of the Rayleigh region. More specifically, both the stochastic ( $ \lambda R < \Lambda $) and geometric ( $ \lambda R \u226a \Lambda $) regimes are investigated. It is worth noting that in the literature, these scattering regimes are studied separately or under different assumptions. In our study, the same approach was used across all three scattering regimes, creating a unified method for studying the scattering of Rayleigh waves from rough surfaces regardless of the scattering regime.
Analytical expressions for the attenuation of Rayleigh waves from statistically rough surfaces are, to our knowledge, very limited in both the geometric and stochastic regions. However, a useful analogous problem has been studied in detail by Van Pamel et al.^{12} In their study, the authors investigated wave scattering within heterogeneous media—where the wave's propagation was impeded by the presence of grain boundaries within a material. In the Rayleigh regime, the authors found a reduction of the attenuation coefficient dependence by one power of frequency between the 2D and 3D scattering—this is consistent with the reduction of the proportionality from f ^{5} to f ^{4} in two dimensions for roughness scattering, derived by Refs. 10 and 11. The same fourth power proportionality between the attenuation coefficient and frequency has also been derived by Kaganova and Maradudin for the scattering of surface waves in a polycrystalline material.^{13} Additionally, it was found that at large λ_{R} values, belonging to the stochastic region, the attenuation coefficient is proportional to the same powers of δ, Λ, and f regardless of the number of dimensions. Therefore, for the problem studied here, the same power relationships derived for the 3D case by Kosachev et al.^{14} can be suggested to hold true in two dimensions. Regarding the geometric region, the authors in Refs. 12 and 15 state that the attenuation coefficient is independent of f.
II. THEORY
This section presents the theory related to both the generation of rough surfaces and power relationships between the attenuation coefficient and different parameters that characterise the incident wave and rough surface. The analytical expressions for calculating the attenuation coefficient, α, along with a brief derivation are also presented and discussed.
A. Rough surfaces
B. Power relationships
An analogous analysis in Refs. 11 and 20 showed that for 2D roughness, the ω_{2} function is proportional to the third power of ω and Λ in the Rayleigh regime. The derivations were performed using similar methods to those in Ref. 7, which allow ω_{2} to be directly substituted into the righthand side of Eq. (7). Therefore, the theory predicts that $ l 2 D R \u2212 1 \u221d f 4 \delta 2 \Lambda $, where the subscript “2D” denotes the case of twodimensional (2D) roughness.
For the stochastic region, in an inhomogeneous medium in three dimensions, where $ \lambda R < \Lambda $, it has been shown analytically in Ref. 12 that there is no difference in the power relationships in three dimensions and two dimensions between $ l \u2212 1$ and the relevant parameters—this observation was subsequently verified numerically by the authors. In this study,^{12} the authors derived this theory for the attenuation arising from inhomogeneous media—however, their analysis used similar metrics to ours. The size of the “obstacle” causing the attenuation was characterised by its correlation length, and the stochastic regime was defined as the region where $ q \Lambda > 1$. We can, therefore, suggest the independence of the attenuation length with respect to dimensionality to hold true in our case as well, i.e., the same power relationship exists between $ l \u2212 1$ and δ, Λ and f in both 3D and 2D roughness. Kosachev et al.^{14} have derived an analytical expression for the stochastic region for 3D roughness—therefore, based on their derivation, we expect that $ l 2 D S \u2212 1 \u221d f 2 \delta 2 \Lambda \u2212 1$, where the subscript “S” denotes the stochastic region.
To present power relationships relating to the geometric regime, it is first necessary to introduce the dimensional analysis associated with the asymptotic study of the attenuation coefficient. The asymptotic study is based on the principle of similitude, which stipulates that the study of different phenomena can be treated using equivalent equations if they can be described by the same dimensionless variables.^{6} For studying attenuation phenomena, an equivalent mathematical analysis can be implemented if roughness parameters (δ and Λ) and the loss in energy (α) are normalised by the Rayleigh wavelength.
Equation (11) holds true regardless of the number of dimensions and the scattering regime—this can also be confirmed by observing both the 3D and 2D power relationships demonstrated in the previous paragraphs, which all obey Eq. (11). In the geometric regime, the scattering is independent of the frequency.^{15} Therefore, from Eq. (11), $ m \delta G = \u2212 1 \u2212 m \Lambda G$ and $ l 2 D G \u2212 1 \u221d \delta m \delta G \Lambda \u2212 1 \u2212 m \delta G$, where the subscript “G” denotes the geometric regime. Here, the power coefficients relating to δ and Λ will be treated as unknowns to be found, and the independence of the attenuation coefficient on f will be validated by the FE model. Finally, it is worth noting that for all of the theoretical derivations, it was assumed that $ \delta < \Lambda $, which follows numerous studies on rough surface scattering in the literature such as Refs. 18 and 21–23.
A summary of the expected power relationships between α and δ and between Λ and f is shown in Table I. For ease and uniformity of presentation, we have introduced the dimensionless notation δ_{n}, Λ_{n}, α_{n}, and β, where $ \alpha n = \alpha \lambda R$ and $ \beta = \alpha n \Lambda / \delta $. The variable β is defined to later allow us to generate a master curve in which we plot the numerical results against a single variable. If just the conventional normalised attenuation coefficient (α_{n}) is used, this yields results which are functions of both δ_{n} and Λ_{n}, making it impossible to plot all of the results in a single graph as they are not functions of a single variable. The variable β is defined such that all of the numerical results become a function of a sole variable (δ_{n}), assuming that $ m \delta G$ is zero.
Regime .  Rayleigh .  Stochastic .  Geometric . 

Limits  $ q \delta < q \Lambda < 1$  $ q \delta < 1 < q \Lambda $  1 $ < \u2009 q \delta \u226a q \Lambda $ 
$ \alpha ( \delta , \Lambda , f )$  $ \delta 2 \Lambda f 4$  $ \delta 2 \Lambda \u2212 1 f 2$  $ \delta m \delta G \Lambda \u2212 1 \u2212 m \delta G$ 
$ \beta ( \delta n , \Lambda n )$  $ \delta n \Lambda n 2$  δ_{n}  $ \delta n m \delta G \u2212 1 \Lambda n \u2212 m \delta G$ 
Regime .  Rayleigh .  Stochastic .  Geometric . 

Limits  $ q \delta < q \Lambda < 1$  $ q \delta < 1 < q \Lambda $  1 $ < \u2009 q \delta \u226a q \Lambda $ 
$ \alpha ( \delta , \Lambda , f )$  $ \delta 2 \Lambda f 4$  $ \delta 2 \Lambda \u2212 1 f 2$  $ \delta m \delta G \Lambda \u2212 1 \u2212 m \delta G$ 
$ \beta ( \delta n , \Lambda n )$  $ \delta n \Lambda n 2$  δ_{n}  $ \delta n m \delta G \u2212 1 \Lambda n \u2212 m \delta G$ 
III. SETTING UP THE FE MODEL
This section presents the method used to create the FE models, which were used in all subsequent simulations. The purpose of the FE models is to allow us to study the phenomenon of surface wave scattering from a large range of roughness parameters meaningfully and efficiently. Therefore, the necessary steps were taken to minimise the model's size and ensure that the surface wave was of good quality with minimal noise. The process to achieve these properties, as well as the method used to calculate the attenuation coefficient, and the computational resources used are presented below.
An example of the lower portion of the FE domain after the rough surface was attached to its lower boundary is shown in Fig. 2. There are two interesting features in Fig. 2. First, the ability of Pogo to mesh efficiently can be observed as the irregular mesh only lies close to where the rough surface is located, whereas the mesh efficiently reverts to a regular form in the main bulk of the material as the distance from the rough surface increases. Second, the smooth joining of the rough surface to the FE domain can also be seen. Here, it is worth acknowledging that rough surfaces can be described by fractals^{32} with past roughness studies using the Weierstrass function,^{33} which exhibits selfsimilarity, to model them. Meshing unavoidably truncates this fractal feature. However, it is expected that the absence of this fractal nature will not affect the result as geometrical features significantly smaller than the wavelength cannot be resolved by the wave.^{22,34}
In the FE model, the input signal used was a fivecycle Hann windowed tone burst. The absolute value of signal's amplitude was arbitrarily selected because the simulation is linear and we are concerned about ratios of measured results and not their absolute values. A source line, comprised of multiple source nodes, was located to the left of the rough surface. For each simulation, the size of the source was set to be equal to three Rayleigh wavelengths, calculated from the simulation's central frequency. To obtain a Rayleigh wave travelling toward the rough surface, a phase delay was applied to each node such that constructive interference from the signal from each node occurred in the desired direction. This method was implemented as follows:

Two sinusoidal time domain signals were created with a 90° phase shift between them.

Pogo provides the ability to assign each source node a unique amplitude, which scales the time domain signal assigned to that node accordingly. Therefore, the amplitude at each node was selected in a way such that a clean Rayleigh wave with the correct amplitude and phase was created by the interference of the signal from all of the nodes.
 The complex amplitude assigned the ith source node, a_{i}, located at the x_{i} position in the model is given by$ a i = 1 2 [ 1 \u2212 cos \u2009 ( 2 \pi x i \u2212 x min x max \u2212 x min ) ] e j q x i ,$
where $ x min$ is the position of the leftmost source node, $ x max$ is the position of the rightmost source node, and j is the imaginary unit. The collective use of a unique amplitude at each node, according to Eq. (17), results in the constructive interference of the signals from all of the source nodes in the correct direction, which generates a clean Rayleigh wave.

The first signal was applied to the ith node with a weighting of [Im(a_{i}),Re(a_{i})] and the second signal was applied with a weighting of [Re(a_{i}),Im(a_{i})], where the two entries in the previous vectors denote the x and z direction, respectively, and the notations Re( ) and Im( ) denote the real and imaginary parts of their argument, respectively. The amplitudes of the x and z components of the Rayleigh wave are arbitrary because one is interested in the ratio of the amplitudes before and after the rough surface rather than the absolute values.
An example of a Rayleigh wave created using the method described above is shown in Fig. 3. As shown in Fig. 3, a pure Rayleigh wave is created by implementing this method. The minimal secondary waves that exist in Fig. 3 just above the righthand end of the Rayleigh wave lie away from the rough surface and, therefore, do not interfere with the attenuation measurements.
Finally, to avoid unwanted noise from the source or scattered waves interfering with the attenuation measurements, absorbing layer regions were applied to the top, left, and right sides of the FE domain. These regions are defined by material parameters whose damping gradually increases, leading to attenuation rather than reflection of the waves.^{35} The addition of absorbing layers in similar setups has proven to be beneficial in related rough surface studies.^{36}
For each Monte Carlo simulation, 100 unique surfaces and, hence FE domains, were generated. Using 1 Nvidia GTX 1080Ti (Santa Clara, CA) with 11 GB of memory, each set of 100 models takes about 1.5 h to complete. This efficiency allowed one to run several statistically stable sets of simulations, covering a wide range of roughness, which better validates the theoretical models. A 3D approach would not have been able to achieve such wide variety of simulation parameter values due to computational limitations. The 2D study presented here provides useful and meaningful insight in the verification of the existing 3D theory as the 2D^{11} and 3D^{7} theories have been derived under the same assumptions using similar methods. Additionally, the 2D work here identifies the most important regimes in which 3D investigations may be conducted. It is also worth noting that the same method described here was used to generate all of the FE models required for our study, regardless of the scattering regime. This strengthens the universality of our findings as it has eliminated the need to study each regime separately, which is often the case in the literature.
IV. RESULTS AND DISCUSSION
This section presents the results from our FE simulations and is split into three parts. In the first part, our FE results are compared quantitatively with Eq. (7) to investigate the agreement between the theory and FE model. After this agreement is established, the second part presents the results relating to the power relationships presented in Table I.
The initial roughness statistical parameters used for the Rayleigh regime were selected to reflect values of roughness that can be found on metal parts created by additive manufacturing. More specifically, values of δ in the range of 10–25 μm^{37–39} were chosen. The correlation length values were then adopted such that they fulfilled the $ \delta < \Lambda $ condition as required by the limiting conditions in Table I. When the investigation relating to the Rayleigh regime was completed, we explored the stochastic and geometric regimes by expanding our initial selection of δ and Λ values.
Each surface realisation, despite being characterised by its δ and Λ values, has a unique profile due to the inherent randomness of roughness. It is, therefore, necessary to average over a sufficient number of realisations to get results that are statistically meaningful. Each datapoint in Figs. 5–9 is an ensemble average α value, obtained from the Monte Carlo simulation of 100 realisations. This number of realisations lies within the range of 50–200, which has been used in similar studies^{22,26,31} but has also been verified for its statistical stability for our specific study by conducting a convergence analysis similar to that of Ref. 40. An example of the convergence plots for three roughness cases is shown in Fig. 4.
As shown in Fig. 4, the number of realisations required to obtain a converged α value is a function of the surface roughness. For the case with the longest correlation length in Fig. 4, the α value converges after approximately 15 realisations, whereas for the other 2 cases, where the correlation length is shorter and, therefore, the peaks and troughs of the surface are closer together, a higher number of realisations is required for convergence. However, all cases converged before the 100 realisation limit was reached, further supporting our choice.
Finally, following the verification of the quantitative and asymptotic results, the last part of this section presents a “master” attenuation curve on which attenuation values from various δ, Λ, and f combinations are plotted as a method to further verify the agreement of the FE and theory over a wider range of parameters. A summary of all of the power relationships obtained by FE modelling is also given.
A. Quantitative results
A comparison between the theoretical α and values predicted by combining Eqs. (7) and (12) is shown in Fig. 5. The statistical parameters used to obtain Fig. 5 were $ \delta = 200 \u2009 \u2009 \mu m$, $ \Lambda = 800 \u2009 \u2009 \mu m$, and the frequency was varied from 0.4 to 1.75 MHz.
It is clear that the FE results follow the curve predicted by the theory. The agreement is also evident in more than one of the scattering regimes described in Sec. II B because Fig. 5 contains results for $ q \Lambda $ values both greater and smaller than one, which is defined as the value at which the scattering behaviour changes from the Rayleigh to the stochastic regime. Using the definitions in Table I, the results of Fig. 5 cover the Rayleigh and stochastic regions with some of the points at higher $ q \Lambda $ values lying in the transition region between the stochastic and geometric regimes. It is worth noting that there are some very small discrepancies between the theoretical attenuation coefficient and FE results. However, we are looking at differences between two approaches here, so it would be inappropriate to view this as errors in the FE simulation not matching the theory. We believe there are three possible sources of discrepancies—approximations in the theory, errors in the FE simulations, and insufficient convergence of the attenuation coefficient at 100 realisations. The outcome of the recent studies using highfidelity FE analysis is that the results from the FE modelling are highly accurate. The authors in Refs. 26–28 discuss in detail the high degree of accuracy achieved through FE, hence, the error associated with the FE approach here is expected to be very small. Regarding the insufficient convergence issue, the results in Fig. 4 show good convergence, hence, we expect this error to also be small. Therefore, it is possible, indeed likely, that the approximations in the theory are a bigger contributor to the differences between the theory and FE results. This confirms the validity of our FE model and allows us to proceed with a deeper analysis of each scattering regime.
B. Power relationships
In this section, we present the results relating to the power relationships in Table I. To calculate the power relationship between α and the variable of interest, the sought power coefficient was determined numerically by a least squares regression analysis, and the results are plotted on a loglog scale in which the power coefficient of the best fitting power function is represented by the slope of the regression line. Similar to Fig. 5, SE bars have been added to all of the plots relating to the power relationships.
1. Rayleigh regime
The results relating to the Rayleigh region are shown in Fig. 6. We are plotting against the characteristic parameters f, δ_{n}, and Λ_{n}, which are the frequency and normalised RMS height and correlation length, respectively, as defined in Sec. II B. According to the theory presented in Sec. II B, we are expecting $ \alpha \u221d f 4 \delta 2 \Lambda $ in the Rayleigh regime. Figure 6(a) shows the simulation data relating to the f ^{4} relationship. To generate the datapoints in Fig. 6(a), the rough surfaces were defined to have δ = 10 μm and $ \Lambda = 20 \u2009 \u2009 \mu m$. The frequency was varied from 1.75 to 2.25 MHz, and at each frequency point, a Monte Carlo simulation of 100 realisations was completed. The gradient of the best fit line in Fig. 6(a) is 3.77, which is close to the expected value of four.
Figure 6(b) shows the simulation results relating to the $ \delta 2$ relationship. To produce Fig. 6(b), the frequency of the simulations was set to 0.5 MHz ( $ \lambda R = 5800 \u2009 \u2009 \mu m$), and Λ was set to 80 μm. Then, δ was varied from 30 to 80 μm. The gradient of the best fit line in Fig. 6(b) is 1.77, which is fairly close to two. It is worth noting that the power relationship here is calculated for $ \alpha n ( \delta n )$, however, it holds true against δ as well—this is because the simulations were completed at a fixed frequency (hence, wavelength) and, therefore, the values in Fig. 6(b) have been normalised by the same scalar.
Figure 6(c) shows the simulation results relating to the Λ relationship. To produce Fig. 6(c), the frequency of the simulations was set to 0.5 MHz ( $ \lambda R = 5800 \u2009 \u2009 \mu m$) and δ was set to 20 μm. Then, Λ was varied from 50 to 90 μm. The gradient of the best fit line in Fig. 6(c) is 0.94. Again, although the power relationship has been calculated for $ \alpha n ( \Lambda n )$, it remains the same for Λ for the same reason explained in the previous paragraph.
It appears that our FE results match the expected power relationships presented in Table I very closely. The power relationships can be explained as follows. The analysis in Maradudin and Huang^{11} demonstrated that the ω_{2} function is proportional to $ ( \omega \Lambda ) 3$ as $ q \Lambda $ tends to zero. When this is substituted into Eq. (7), it yields the $ \alpha \u221d f 4 \delta 2 \Lambda $ relationship, demonstrated by our FE results here. Physically, the attenuation coefficient tends to zero at low frequencies because the Rayleigh wavelength becomes so long relative to the statistical parameters characterising the surface that the surface appears flat to the wave.
2. Stochastic regime
The results relating to the stochastic region are shown in Fig. 7. Based on the theory presented in Sec. II B, a proportionality of α to $ f 2 \delta 2 \Lambda \u2212 1$ is expected in the stochastic regime.
Figure 7(a) shows the FE results relating to the f^{2} relationship. For this set of simulations, δ was set to 200 μm and Λ was set to 800 μm. The gradient of the best fit line is 1.94. To produce Fig. 7(b), the frequency of all simulations was set to 1 MHz ( $ \lambda R = 2900 \u2009 \u2009 \mu m$) and Λ was set to 800 μm. Then, δ was varied from 100 to 400 μm, satisfying the stochastic region's condition. The gradient of the best fit line was found to be 2.15. Finally, Fig. 7(c) was generated by setting the frequency again to 1 MHz. Then, δ was fixed to 200 μm and Λ was varied from 400 to 1600 μm. The gradient of the best fit line was found to be −1.25.
Overall, it appears that our FE model follows the already established theory in the stochastic region well. The power relationships here can be again explained by the behaviour of the ω_{2} function in the stochastic region. As has been demonstrated in Ref. 11, in the stochastic region, $ \omega 2 \u221d \omega \Lambda $. When this relationship is substituted in Eq. (7), it yields the $ f 2 \delta 2 \Lambda \u2212 1$ proportionality predicted by the theory and supported by our FE results.
3. Geometric regime
The results relating to the geometric region are shown in Fig. 8. The geometric region is a region in which the RMS height is greater than $ \lambda R / 2 \pi $ as per the definition made in Table I. Therefore, the frequency in this set of simulations was set to 5 MHz ( $ \lambda R = 580 \u2009 \u2009 \mu m$). Then, δ was set to 200 μm and Λ was varied from 800 to 1600 μm. The gradient of the best fit line in Fig. 8 is −0.83.
The attenuation coefficient can be seen to be decreasing with an increase in the frequency in Fig. 8. Physically, in this region, the wavelength has become so small compared with the correlation length that the roughness does not impede its motion—the wave travels along the peaks and troughs without being scattered. Additionally, from the dimensional analysis presented in Sec. II B, the fact that $ m \Lambda G$ is approximately equal to one implies that $ m \delta G$ is approximately equal to zero in the geometric regime. This is further investigated and validated in Sec. IV C.
C. Summary of results
In Table I, we have also introduced the generalised attenuation coefficient, β. Using this allows us to further verify the validity of the theory in the stochastic and geometric regions, by plotting a wider range of Monte Carlo results against only the variable δ_{n}. As shown in Fig. 9, the results follow the asymptotic approximation lines independently of which parameter was the variable in each FE Monte Carlo set. Additionally, the transition between the stochastic and geometric region can clearly be seen at $ \delta n = 1 / 2 \pi $. This is the expected location of the transition and can be derived by identifying that in Table I, the transition point between the stochastic and geometric regimes is defined to be where $ q \delta = 1$.
It is now worth noting that the gradient of the asymptote on the righthand side of the plot, which corresponds to the geometric regime, is equal to −1. Therefore, from our dimensional analysis and the values in the last row of Table I, $ m \delta G$ is again found to be equal to zero because $ m \delta G$ − 1 = −1. This also confirms that $ m \Lambda G = \u2212 1$, which is demonstrated by our FE results, indicates that the attenuation coefficient is independent of the frequency in the geometric regime, because $ m f G$ must be approximately equal to zero from Eq. (11). A summary of the power relationships obtained from the FE results in all scattering regimes is shown in Table II.
Regime .  Rayleigh .  Stochastic .  Geometric . 

Limits  $ q \delta < q \Lambda < 1$  $ q \delta < 1 < q \Lambda $  $ 1 < q \delta \u226a q \Lambda $ 
$ \alpha ( \delta , \Lambda , f )$  $ \delta 1.77 \Lambda 0.94 f 3.77$  $ \delta 2.15 \Lambda \u2212 1.25 f 1.94$  $ \Lambda \u2212 0.83$ 
Regime .  Rayleigh .  Stochastic .  Geometric . 

Limits  $ q \delta < q \Lambda < 1$  $ q \delta < 1 < q \Lambda $  $ 1 < q \delta \u226a q \Lambda $ 
$ \alpha ( \delta , \Lambda , f )$  $ \delta 1.77 \Lambda 0.94 f 3.77$  $ \delta 2.15 \Lambda \u2212 1.25 f 1.94$  $ \Lambda \u2212 0.83$ 
Comparing Tables I and II, it is clear that there is good agreement in the asymptotic power relationship coefficient across all scattering regimes. This has two implications. First, we have managed to verify the wellestablished theory regarding scattering in the Rayleigh regime both quantitatively and asymptotically. Second, our FE model was able to also verify the asymptotic in the stochastic and geometric regimes, confirming the applicability of assumptions from scattering^{12} to our study.
V. CONCLUSION
A comprehensive study of the attenuation of Rayleigh waves from 2D statistically rough surfaces using FE modelling has been presented. Three distinct scattering regimes have been identified from the literature—the Rayleigh (low frequency), stochastic (low to medium frequency), and geometric (high frequency) regimes. Analytical formulas, predicting attenuation values, have been derived in the past,^{11} as well as asymptotic power relationships^{10,11} between α and δ and between Λ and f.
Here, we attempted to validate the existing theory using the FE analysis and extend the results to regions in which the theory is less established or obtain results with a wider combination of δ, Λ, and f values. We have found good agreement between the theory and FE results in all three regimes—in the Rayleigh and stochastic regimes, good agreement was found both quantitatively and asymptotically, and the f ^{4} relationship between α and the frequency in the Rayleigh regime was also observed. For the geometric regime, power relationships were derived by a combination of FE modelling and dimensional analysis.
The FE model's ability to follow the theory creates a plethora of useful implications. The theoretical formulas rely heavily on the ω_{2} function, which is a complicated function comprising multiple subfunctions, many of which have a different form, depending on the region of interest, meaning that calculating ω_{2} is far from straightforward. The ω_{2} function's validity is also limited in terms of the roughness parameters for which it can produce results ( $ \delta / \Lambda < 0.3$), and its behaviour has not been studied extensively in the literature in the geometric regime. FE modelling removes the necessity to obtain this function and allows for direct calculation of the attenuation coefficient. Additionally, the FE models can potentially be extended to regimes where the literature is more limited, such as the geometric regime, and can also simulate δ and Λ parameters outside the theory's region of validity. Finally, the use of FE modelling has provided a more unified approach to the study of rough surface scattering. In the literature, each scattering regime is largely studied in depth on its own, whereas the FE approach here has been able to verify the theory in all three regimes by always implementing the same method.
Finally, it is worth noting that despite the results here being obtained from 2D simulations, they are still relevant for the 3D analytical formulas. The mathematical approach used to derive the 3D^{7} and 2D^{11} theories are analogous—therefore, the FE validation of the 2D theory in our study provides important insight for the 3D theory, indicating that it will also hold true for FE simulations and experimental scenarios.
ACKNOWLEDGMENTS
G.S. is funded by the United Kingdom (UK) Research Centre in nondestructive evaluation (NDE), iCASE No. 17000191, with contributions from RollsRoyce Holdings plc and Jacobs Engineering Group Inc. During this work. P.H. was partially funded by the UK Engineering and Physical Sciences Research Council (EPSRC) Fellowship No. EP/M020207/1. M.L. is partially sponsored by the EPSRC.