## Abstract

The surface roughness of soil grains affects the mechanical behaviour of soils, but the characterization of real soil grain roughness is still limited in both quantity and quality. A new method is proposed, which applies the power spectral density (PSD), typically used in tribology, to optical interferometry measurements of soil grain surfaces. The method was adapted to characterize the roughness of soil grains separately from their shape, allowing the scale of the roughness to be determined in the form of a wavevector range. The surface roughness can be characterized by a roughness value and a fractal dimension, determined based on the stochastic formation process of the surface. When combined with other parameters, the fractal dimension provides additional information about the surface structure and roughness to the value of roughness alone. Three grain sizes of a quarzitic sand were tested. The parameters determined from the PSD analysis were input directly into a Weierstrass–Mandelbrot function to reconstruct successfully a fractal surface.

## 1. Introduction

The surface of soil grains is not smooth, especially when examined at increasingly smaller scales. In the field of geotechnical engineering, surface roughness has been shown to affect the packing, shear modulus, compression and shearing behaviour of granular assemblies [1–5], but typically these studies were made using analogue soil grains such as glass ballotini or steel balls which surface was altered by chemical or mechanical action. Characterizing the surface roughness of real soil particles, although pivotal to any quantitative analysis of its effect, remains rare (e.g. [5,6]; table 1).

The significant role of surface roughness in the contact behaviour between objects has led to a large body of research in the fields of tribology and industrial manufacturing. In their pioneering work, Greenwood & Williamson [7] showed that real contact surfaces with asperities have larger contact areas and smaller contact pressures than predicted by Hertz [8] classic solution for the elastic contact of smooth spheres. Advanced experimental imaging methods have been developed to determine accurate surface topographical information, including mechanical stylus profilometry, optical interferometry, scanning electron microscopy and atomic force microscopy. These methods are generally applied to usually flat, engineered materials, e.g. milled, sand-blasted or thin-coated. Materials investigated by interferometry tend to have good reflectivity.

In soil mechanics, there has been a growing effort to study and model the behaviour of granular geomaterials at the particle level (e.g. [3,9,10] and subsequent DEM studies). Advances in the theory of contact mechanics of rough curved contacts, such as that proposed by Greenwood *et al*. [11], have been used in analytical and numerical analyses of soils [4,12,13]. These require a representative value of the grain surface roughness, but unlike with engineered materials, the use of some of the advanced techniques developed in other fields can be difficult to apply to soil particles, which can have a variety of shapes and roughness resulting from their mineralogy and diagenetic geological history. Depending on their mineralogy, their reflectivity can also be very low (e.g. quartzitic grains).

Sands may have diverse origins, typically clastic or bio-clastic, which has a marked effect on their nature. Some have suggested that the roughness should be a proportion of the size of the particle (e.g. [3]), which finds rationale in the effect of grain size on the processes of transportation and deposition of the soil, but the relationship between the two is not obvious. Both roughness and shape are affected by the geological origin of the grain (e.g. biogenic, erosion of igneous, metamorphic or sedimentary rock), its mineralogy and therefore hardness, but probably in slightly different ways. Particle breakage, chipping and abrasion during particle loading or transportation, perhaps related to the grain type and the transportation type (e.g. wind, water and ice) will affect the grain morphology. Particle breakage, e.g. splitting may create more angular grains but depending on the mineralogy the created surfaces will be more or less smooth. For example, quartz breaks along conchoidal surfaces, while the crystalline structure of feldspar forces it to break along cleavage planes (e.g. [14]). On the other hand, chipping will make particles less angular. Chemical effects in the long term, which are most likely to occur after deposition, such as dissolution, either generalized or at particle contact, and modification (e.g. precipitation of iron oxide), also affect the shape and roughness in different ways. The most typically encountered materials in geotechnical engineering applications are quartz sands such as the one tested here. These have an igneous origin but which may have been through cycles and various types of weathering, erosion, transportation and deposition as a sand/sandstone so that they may have a wide range of ages.

So far characterizing soil grain surface roughness has been either by estimating the RMS (*S _{q}*) or average (

*S*

_{a}) from a cut section of the grain surface (e.g. [3,5]), or by using two-dimensional grain profiles obtained from scanning electron microscopy, but being able to visualize the whole grain has been at the cost of losing the resolution and significant detail of the surface roughness (e.g. [15–17]). Some studies have made successful use of advanced technology, such as optical interferometry, but the analysis of the measured data has been simplistic in comparison with the challenge of obtaining the data (e.g. [3,5,6,9,18,19]). The amount and quality of data on real soil particles remains small, limiting further application of the results to numerical modelling at the grain scale (e.g. [20–24]).

In order to exploit topographical measurements of soil grain surfaces better, we have adapted a method used to determine the roughness of engineered surfaces to use on particles from a natural quarzitic sand. It is found that analysing measurement data obtained from high-resolution optical interferometry as a power spectrum can lead to a more informative yet objective quantification of roughness than currently achieved. The surface is thus described by a scale-independent parameter (the fractal dimension) in addition to the RMS of the roughness, a suggestion that has been made for engineered surfaces (e.g. [25,26]).

Several methods have been proposed in fields ranging from manufacturing to medicine to assess the fractality of surfaces. Scanning electron micrographs can be used, for example, the grey scale of the images allows texture techniques to be applied, such as the ‘skyscraper’ fractal analysis (e.g. [27]) or the ‘blanket’ fractal analysis (e.g. [28]). The projected areas of the particles, obtained from SEM or other means such as image sensor analysers, can be used with the box counting method (e.g. [29]), or the area–perimeter method (e.g. [15]), but for soil grains high-resolution images are necessary to be able to capture the surface asperities. Dividers have been used to determine the fractal dimension of surfaces, such as the triangular-prism method (e.g. [30]), variograms (e.g. [31]), triangulation or cube-counting (e.g. [26,32]). Another technique is to analyse the power spectrum of the surface, for example, as a Fourier power spectrum [33], by power spectral density (PSD) [34] or as a structure function [35]. The fractality of soil grain surfaces has been suggested by researchers who have found grain contours to exhibit a self-similar or self-affine pattern down to finer scales (e.g. [36,37]). In the following, we show how natural sand grain surfaces obtained by profilometry can also be described as fractals.

## 2. Testing apparatus and tested sand particles

The roughness measurements were made with a Fogale Nanotech optical microscope (model M3D 3000). FOGALE Pilot 3D software and FOGALE Viewer 3D software were used to obtain and analyse the data. The surface topography is described by an interferogram that is a function of the sample height at discrete points. The best lateral resolution that can be achieved by this interferometer is 0.184 µm (spacing of discrete points in the *x*- and *y*-planes perpendicular to the surface height plane *h*(*x*, *y*)), while white light profilometry ensures 3 nm RMS resolution in the vertical direction. The measuring area can be up to 141.3 × 106.6 µm. A function available within the integrated software allows separating the shape from the roughness. The function filters separating the low frequencies associated with the shape (e.g. slope, curvatures) from the high frequencies associated with the roughness. The length of this spatial filter, also called motif size, is arbitrarily set as one-quarter of the size of the field of view, and of the same unit as the image unit. The roughness is deducted by subtracting the shape from the overall surface, ensuring that the sum of the shape plus roughness is always equal to the unaltered surface [38]. An illustration of this decurvature process is shown in figure 1.

The tested particles are from Leighton Buzzard sand (LBS), a silica sand consisting of strong, highly spherical particles. The shape parameters were determined by dynamic image analysis using a Qicpic image analysis sensor where soil particles are put through a vibratory feeder to disperse them before free-falling in front of pulsed light. A high-speed digital camera (450 frames s^{−1}) captures images of the particles with a resolution of 1 µm for size and shape characteristics. The sphericity, calculated as the ratio of the perimeter of the grain to that of the circle of equivalent surface area, is about 0.9, and the convexity, calculated as the ratio of the surface area of the grain to the area of the convex Hull surface, is about 1.0.

There has been no systematic study of roughness of soil particles to enable determination of whether their surface roughness is size dependent or not. The question of scaling, i.e. whether the same value of roughness can be used throughout a range of particle sizes, for example to simulate debris flows, is however being queried by discrete element modellers. Three size groups were thus selected to try to add to the limited data available, corresponding to the sieve sizes 0.6–1.18, 1.18–2 and 2–5 mm. It was made sure by visual inspection that the particles tested were of the same mineralogy (quartz). Their particle size distributions determined by dynamic image analysis are shown in figure 2. A total of 150 particles were tested, 50 for each size group. All the surface measurements were made for an area of 106.6 × 106.6 µm corresponding to 578 × 578 discrete points. The low reflectivity of the quartz meant that obtaining good measurements was laborious. A particular difficulty with soil grains is that many points of the irregular surfaces cannot be measured, which are then shown on the resulting graph as fail-to-detect points or invalid pixels. The areas measured were chosen so that invalid pixels in the observed areas were less than 1%, ensuring that removal of these points by interpolation of adjacent heights data had a negligible effect. Edge effects can also be avoided in the same way. Then the surface height data for the 578 × 578 points were exported for analysis.

## 3. Test results and surface characterization

Figure 3 shows four measured surfaces for each size group. Most of the measured surfaces have a relatively spherical local shape with some small irregular dents spread on the surfaces. This is one of the main differences between natural sand surfaces and engineered surfaces where regular wavy curves often exist. By eye, there is no obvious difference between each size group apart from the surface curvature, which is more pronounced for the small grains.

It is implied from the Greenwood & Tripp [39] solution and its variants for non-flat materials that, in order to characterize roughness, a separation procedure is needed to remove the surface curvature from the surface measurement. In optical interferometry, the motif extraction method, which was introduced by Boulanger [40], is generally used as it is integrated in the software of the testing apparatus. The concept of motif, introduced in metrology research in the 1970s for machinery tools manufacturing, refers to the filtering of a surface profile between regular (e.g. waviness) and irregular features associated with the roughness [40]. The shape of real soil grains however does not follow a regular waviness, thus measurements are very sensitive to the separation procedure [6]. For lack of better guidance, the default value of the shape motif, which is available in the software and increases with the size of the measuring area, is generally chosen [41]. Another limitation of the motif extraction method applied to soils is that there is no appreciation of the scale of the asperities, and this can diminish the application of the measured roughness since the range of roughness scales is influential in the interfacial contact behaviour (e.g. [42,43] for flat contacts).

For the measurements made here, the default filter length (i.e. motif size) was 26.7 mm for the measured area of 106.6 × 106.6 µm. This decurvature process was however only applied for comparing the values of roughness compiled with the method presented here (figure 10), as the main focus of this paper is to present a new method which addresses these drawbacks and uses the whole measured dataset is unaltered (i.e. without an artificial separation of shape and roughness). Another advantage of the method, which borrows from tribology, is that instead of a profile line we use the whole measured surface in three dimensions: *x*, *y* and *h* (height).

### (a) Surface morphology by power spectral density

Natural soil grains follow a stochastic forming process, which usually leads to an apparently random surface morphology. Nayak [44] proposed to model rough surfaces as two dimensions, isotropic, Gaussian random processes. He represented the surface morphology by using a spectrum, which can reveal periodic surface features that might otherwise appear random. The PSD is calculated by
*A*(*x*,*y*) is the auto-correlation function of surface heights *h*(*x*,*y*) and *q* is the wavevector or spatial frequency (in µm^{−1}). By using the PSD, the spatial surface heights data are transferred into the spatial frequency or wavevector domain through a discrete Fourier transform. In order to facilitate the interpretation of the surfaces, a routine angular averaging is performed where the surface is assumed to be isotropic so that the PSD(*q _{x}*,

*q*) reduces to PSD(

_{y}*q*) and is independent of the

*x*- or

*y*-direction [44]. This assumes that the PSD is the same along the

*x*- and

*y*-directions. Figure 4

*a*shows the PSDs in each in-plane direction, averaged over 289

*x*-profiles and 289

*y*-profiles, together with the angular averaged PSD. The average PSDs in both in-plane directions are almost coincident, indicating no major anisotropy as can be found in some manufactured surfaces (e.g. [25]). The angular averaged PSD is taking into account all the measurements and therefore not necessarily equal to the average of the two PSDs in

*x*- and

*y*-directions.

Figure 4*b* shows the PSDs for all the 2–5 mm particles, using the angular average PSD (referred to thereafter simply as PSD). No spike or jump is observed, indicating that there is no predominant wavevector in those surfaces. Each PSD curve contains information both about the local shape, at small wavevectors, and the roughness, at large wavevectors, thus potentially enabling the characterization of roughness and shape separately but more objectively than the routine motif method. The variation in PSD within the same size group indicates different roughnesses for different grains, and the average PSD is also shown. The zeroth, second and fourth moments of the PSD relate to physical statistical parameters of the surface [44]. For example, the zeroth moment of the PSD relates to *S _{q}* and can be expressed as [34,44]

*q*

_{0}and

*q*

_{1}denote the smallest and largest wavevectors of the measured surface, respectively. Here, we take the values suggested by Persson

*et al*. [34]:

*q*

_{0}= 2

*π*/

*L*, with

*L*= 106.6 µm (size of the measured area) and

*q*

_{1}= (

*N*/2) (2

*π*/

*L*), with

*L*/

*N*= 0.184 µm (resolution). Figure 5

*a*shows the comparison between the values of

*S*obtained from the PSD after angular averaging and from statistics using

_{q}The values are in very good agreement, with an error less than 0.1% (0.001 µm). This validates the calculation of the PSD and shows that the angular averaging which assumes the surface to be isotropic has a negligible effect on the parameters. A similar graph obtained for different sizes of field of view taken on a grain of size 0.6–1.18 mm, shown in figure 5*b*, indicates that the size of the field of view (i.e. the value of *L*) does not affect the good comparison between the values of *S _{q}* derived from PSD and from the motif method.

The value of *S _{q}* derived with the PSD takes account of the whole measured surface, unaltered, as it takes all wavevectors into account, and therefore encompasses both shape and roughness. The high values of

*S*, between 2 and 20 µm, capture mainly the shape of the grains, even more so in the smaller particles, resulting in larger roughness values for those.

_{q}In order to determine the value of the roughness alone, which we call *S _{q}*

_{,roughness}, the scales at which the roughness acts and at which the shape acts should be determined. We describe below how we simply use a cut-off wavevector,

*q*

_{c}, that relates to the largest wavelength that contributes to the surface roughness, to separate roughness and shape in the PSD. The value of

*q*

_{c}depends on the particular surface measured and therefore should vary from particle to particle.

## 4. Separation between shape and roughness surfaces

### (a) Determination of the cut-off vector *q*_{c}

As a first step the sensitivity of the value of the surface roughness *S _{q}*

_{,roughness}to the cut-off wavevector

*q*

_{c}is investigated. Figure 6 shows the three average PSD curves for the three size grain groups, with four different values of

*q*

_{c}for wavelengths between 1.3 and 5.3 µm. According to equation (3.1), by replacing

*q*

_{0}with

*q*

_{c}, i.e. only considering the data between

*q*

_{c}and

*q*

_{1}, we should obtain the value

*S*

_{q}_{,roughness}. Increasing

*q*

_{c}four times has the effect to reduce the value of

*S*

_{q}_{,roughness}by up to 48% (table 2), thus, although more suited to soil grains than the motif extraction, caution should be taken when determining

*q*

_{c}to obtain a reliable value of roughness.

We use the idea that the surface area calculated by discretizing into a grid, which will increase as the grid mesh size decreases, should show a marked increase when features associated with the roughness of the surface are captured by the grid. The surface area is estimated geometrically using the TPM [30]. The TPM being more suitable for self-similar surfaces [45], which occur rarely in nature (e.g. [46]), here we use it primarily as a means to separate shape and roughness. Figure 7 shows a visual illustration of a typical surface that is discretized with mesh sizes *δ* proportional to the resolution: [64, 48, 16, 10, 4, 1] × 0.184 µm. The smaller the *δ*, the more the discretized surface will match the measured surface. Figure 8*a* shows the estimated surface area changes with grid sizes of [64, 48, 24, 16, 10, 8, 6, 5, 4, 3, 2, 1] × 0.184 µm, for grains 2–5 mm, using a finer grid around 2 µm in order to determine the cut-off size delimiting shape from roughness and to provide more information in the small *δ* region. Using finer grids, for example, the finest [64: −1:1] (i.e. [64, 63, 62, …, 3, 2, 1] × 0.184 µm), would lead to highly fluctuating values of surface area, as shown in figure 8*b*. The estimated surface area decreases slowly down to about *δ* = 1.84 µm, accelerating as grid sizes decrease from 1.84 to 0.184 µm. This may indicate that asperities of surface roughness enter the surface area estimation for those small grid sizes. This is consistent with the suggestion that the morphology of irregular particles might be best described over two fractal subsets characterizing the texture (≈roughness) and structure (≈shape) (e.g. [15,35,36]). The grid size *δ*_{c} ≈ 1.84 µm, regarded as the maximum size of asperity or the minimum wavelength, is taken to be the cut-off grid size between shape and roughness, the smaller grid sizes defining the roughness. The parameter *δ*_{c} does not change with particle size for the sand and size range investigated and represents between 0.06 and 0.3% of the dimension of the soil grains.

The values of *S _{q}*

_{,roughness}are calculated using the PSD (equation (3.2)), with a lower cut-off value of

*q*

_{c}≈ 2

*π*/

*δ*

_{c}= 3.42 µm

^{−1}determined from the estimated surface area above. They are plotted in figure 9. For the two larger-sized grains, the majority of the roughness values is between 0.35 and 0.45 µm, while the smaller grains have values of roughness varying over a larger range. The distributions are not normal and they are more peaked for the larger particles. The average value is 0.66 (±0.29) µm, 0.46 (±0.14) µm and 0.36 (±0.08) µm for particles of sizes from 0.6 to 1.18 mm, 1.18 to 2 mm and 2 to 5 mm, respectively, i.e. while

*δ*

_{c}does not change with particle size,

*S*does. Cavarretta

_{q}*et al*. [3] reported values of

*S*

_{q}_{,roughness}for LBS sand grains of 0.7–2.2 mm diameter of 0.3 µm and Senetakis

*et al*. [9] reported 0.38 (±0.124) µm for grains of 1.18–5 mm. In both cases, the measured areas were smaller, which might affect their results [6,41], and the motif extraction method was used. Otsubo

*et al*. [6] showed that the standard deviation increases with reducing size of the measuring area. The values obtained from the PSD method and the motif extraction method using the default input value in the software (shape motif of 26.7 µm) are compared in figure 10: for the small particle sizes (0.6–2 mm) the motif extraction method underestimates

*S*

_{q}_{,roughness}, which indicates that a smaller value of

*δ*

_{c}than 1.84 µm, i.e. a smaller roughness scale, was used in the motif method, while for the larger particle sizes the motif extraction method gives a higher value of

*S*

_{q}_{,roughness}, indicating that a larger roughness scale was used. The value of the motif depends on the size of the measured area, independently of the material tested, while the cut-off wavelength has a physical relation to the surface of the grains tested. In that sense, even if sometimes the values from the two methods are comparable, the method presented here gives more control and physical meaning over the value of

*S*

_{q}_{,roughness}.

## 5. Fractal approach to characterize surface roughness

The morphology of soil grains can be characterized successfully by fractal dimensions (e.g. [36,37]), which can be used to recreate realistic numerical models of surfaces (e.g. [22]). Several researchers have proposed and compared different methods for the fractal analysis of surfaces (e.g. [47,48]). Here, we use and compare two methods: the PSD method, based on the stochastic formation process, and the TPM, based on the geometry of the surface. The PSD method is suited to natural surfaces, which typically exhibit different scales in the vertical and in-plane directions (e.g. [49]) while the TPM should be most suited to self-similar surfaces (e.g. [45]). It is therefore expected that the TPM may not be as conclusive as the PSD method to characterize the surface of natural soil grains. Unlike in previous studies of soil grains, which were carried out on the profile of the particles, here the analysis was performed on three-dimensional surfaces.

The surface area calculated by the TPM method, if fractal, is related to the grid size by [30]:
*D*_{TPM} is the fractal dimension determined by TPM. Figure 11*a* shows the averaged surface area against grid size for each size group. The slopes of the lines, equal to 2 − *D*_{TPM}, decrease with increasing grain size, indicating an increase in *D*_{TPM} with particle size and therefore an increase in the estimated real surface area for higher fractal dimensions. Average curves are shown to highlight the trends. The distribution of *D*_{TPM} values obtained from all the measurements is reported in figure 11*b*. The values of *D*_{TPM} for sizes 0.6–1.18 mm are slightly higher than those for 1.18–2 mm, while the values for sizes 2–5 mm are much smaller than the other two size groups. The distributions of *D*_{TPM} have the same trend as the distribution of *S _{q}*

_{,roughness}, which might be expected as the method is based on the graphical representation of the surface and a rougher surface is usually accompanied with a higher surface area, but here the distribution for the smaller particles is narrower. Zhai

*et al*. [26] also found on aluminium discs treated by mechanical attrition that the fractal dimension compiled by triangulation or cube-counting and the RMS roughness have similar trends, although in their case the roughness was higher for the larger particles.

Similarly, a fractal surface would imply that the PSD follows a power law when *q* > *q*_{c} ≈ 3.42 µm^{−1}, with the slope of the spectrum related to the fractal dimension. This is expressed as
*α* (<0) is the slope of the straight fitting line in the double logarithmic plane of PSD versus *q*, and *C*_{0} is related to the intercept. The fractal dimension, *D*_{PSD}, can be expressed in terms of *α*, and several relationships have been proposed: Voss [50] and Turcotte [51] adopted the relationship *D*_{PSD} = 4 + *α*/2 for a three-dimensional surface, with *D*_{PSD} the fractal dimension determined by PSD. For fractal surfaces generated using Weierstrass–Mandelbrot (WM) functions, the relationship *D*_{PSD} = 2.5 + *α*/2 was adopted by Majumdar & Bhushan [52] for two-dimensional surfaces and *D*_{PSD} = 3.5 + *α*/2 was adopted, e.g. by Liou & Lin [53] for three-dimensional surfaces.

Figure 12*a* shows the fractal parameters *D*_{PSD} calculated with *D*_{PSD} = 3.5 + *α*/2 from the averaged PSDs for each size group; the distributions of *D*_{PSD} for all particles are shown in figure 12*b*. The values of *D*_{TPM} and *D*_{PSD} for a given particle size are distributed differently. It is found that *D*_{PSD} increases with increasing particle size, which shows an opposite trend to *D*_{TPM} and *S _{q}*

_{,roughness}. The closer fit to a linear equation for the PSD data (figure 12

*a*) compared with the surface area data (figure 11

*a*) seems to indicate that the surfaces satisfy self-affinity rather than self-similarity implied by the TPM.

If the fractal rule applies, combining equations (3.2) and (5.2) leads to
*C*_{0}, the fractal dimension *D*_{PSD} as well as *q*_{c}/*q*_{1} and *q*_{c}. If we plot the values of *S _{q}*

_{,roughness}determined from the PSD (equation (3.2)) against those of

*C*

_{0}(figure 13), we can define a unique relationship between the two parameters which can be described as

*a*= 5.95 and

*b*= 0.41 the best-fitting values. There is a slight deviation from equation (5.3), for which the exponent of

*C*

_{0}is equal to 0.5. Deriving

*S*

_{q}_{,roughness}from equation (5.3) with the average values of

*D*

_{PSD}and

*C*

_{0}in figure 13 gives consistently lower values than when using equation (3.2) (table 2). The values of

*C*

_{0}associated with the average values of

*D*

_{PSD}, also in table 3, reflect the change in roughness as well. The fractal dimensions may depend on the resolution since the measured surface heights are sampled at discrete lengths, but the hierarchical structure of the surface evidenced by the straight lines in the roughness range may be independent of the instrument resolution. Persson [54] showed that the slope of the linear part of the PSD obtained from several testing methods with resolution ranging in magnitude from around 100 to 0.01 µm is rather consistent.

### (a) An example of reconstructed surfaces using the obtained parameters

Numerous studies have shown the strong dependence of interfacial mechanical properties on the PSD, the value of *D*_{PSD}, or the range of roughnesses for flat surfaces (e.g. [22,55]). The WM function and its variants can be used to reconstruct and analyse self-affine surfaces (e.g. [22]), but the identification of the input parameters can be a difficult task. A realistic reproduction of the surface relies on realistic input parameters. Following Majumdar & Bhushan [52] and Yan & Komvopoulos [56], we use the parameters determined from the PSD of the soil surface to reconstruct the surfaces using a variant of the WM function proposed by Yan & Komvopoulos [56]. The parameters determined experimentally above are input directly where possible. The WM function is expressed as
*D* is *D _{PSD}*, the length factor of the highest asperity spacing, or sample wavelength,

*L*, is taken as

*δ*

_{c}= 1.84 µm, the smallest length is

*L*

_{s }= 0.184 µm and a parameter related to the density of frequencies,

*γ*, is set to a commonly adopted value of 1.5, which leads to the maximum frequency index

*n*

_{max}= 6. The number of superposed ridges used to construct the surface,

*M*, is not a first-order parameter [56] and is set to be 20 to make the surface structure sufficiently random. The influential parameters are

*D*(=

*D*), taken as the average values of 2.22, 2.37, 2.41 for the three grain sizes (figure 11

_{PSD}*a*), and the associated coefficient

*C*

_{0}= 6.3 × 10

^{−3}, 2.4 × 10

^{−3}and 1.3 × 10

^{−3}µm

^{4}. The fractal roughness,

*G*, can be calculated from Liou & Lin [53]:

*C*=

_{p}*C*

_{0}

*q*

_{c}

^{−α}. Values of

*G*equal to 6.0 × 10

^{−3}, 8.3 × 10

^{−3}and 6.3 × 10

^{−3}µm are found for the three grain sizes, respectively, in increasing size order.

*φ*is the random phase angle and is chosen so that the equivalent fractal surfaces have roughness values

_{m,n}*S*

_{q}_{,roughness}equal to 0.70, 0.48 and 0.37 µm, in order of increasing grain size. The reconstructed surfaces are of area 9.2 × 9.2 µm, so they fit within the bounds of the cut-off grid size and should be dominated by roughness (figure 14

*a*). They represent a portion of the measured area, and although flat it would be possible to wrap it on a shape as shown by, for example, Liou

*et al*. [57] or Hanaor

*et al*. [23], who modified the WM to overlay a planar rough surface onto a sphere. The WM surfaces were reconstructed while controlling their roughness to be the same

*S*

_{q}_{,roughness}as the real sand particles tested. Their small area ensured that the shape was not affecting the comparison between the measured and reconstructed roughness surfaces (figure 14

*b*). Some differences are observed in the location and heights of the asperities, but although reconstructing the complex surface of a soil grain exactly as it is may not be possible, there is a statistical resemblance. This highlights one of the problems with determining the surface roughness of sand grains, which is that no two grain surfaces are the same. The limited size of the measuring areas, which despite being the largest possible in this paper is still much smaller than the whole grain surface, is also a drawback. It seems likely that the variation in surface roughness on a given grain may be less than the variation in shape, but one should be aware of these limitations when using the measured data.

## 6. Conclusion

A method to characterize the surface roughness of soil grains is proposed. The PSD, a powerful tool to reveal the periodic feature of a random surface which is typically used in tribology, was adapted to characterize the surface roughness of soil grains separately from their shape, a procedure less straightforward than for engineered surfaces. The scale of the roughness, information usually missing in other methods of determining soil grain roughness, was quantified in the form of a wavevector range determined from the estimation of the surface area by triangular-prism method. The surface roughness was then characterized over that range using the PSD.

For three sizes of quartzitic sand particles considered here, the surface roughness has been characterized by a roughness value, *S _{q}*

_{,roughness}, and a fractal dimension, determined from the PSD (

*D*

_{PSD}), which is more suitable for natural surfaces where different scales exist in the vertical and in-plane directions. The

*D*

_{PSD}, when combined with other parameters, such as the coefficient

*C*

_{0}from the PSD, carries more information about the surface structure and roughness than the value of roughness alone.

The obtained data contribute to the current very limited database on real soil particles, with potential use in numerical modelling for creating numerical particles and simulating interfacial grain contacts. A variant of WM function was used successfully to reconstruct a fractal surface with the parameters identified experimentally as direct input.

## Data accessibility

Surface data measured by interferometry on 15 particles (five for each size) and which support this article can be accessed at http://discovery.ucl.ac.uk/1508516/.

## Authors' contributions

H.Y. collected some surface data, performed the analyses and drafted the manuscript. T.Y. collected a large amount of surface data and has contributed to some figures. Both were working under B.A.B.'s supervision. All authors have read and approved the final version of the paper.

## Competing interests

There are no competing interests in this research.

## Funding

The authors acknowledge the financial support provided by the Research Grants Council (RGC) of HKSAR (grant nos. GRF 17200114 and TR22-603-15N).

## Acknowledgements

Prof. B. N. J. Persson is gratefully acknowledged for helping with early calculations of the PSD and for his advice. We also would like to thank Prof. Matthew Coop and Prof. Jidong Zhao for fruitful discussions.

- Received June 28, 2016.
- Accepted September 5, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.