Laboratory-scale characterization of saturated soil samples through ultrasonic techniques

The propagation of poroelastic waves in a soil specimen is dependent on the physical and mechanical properties of the soil. In the geotechnical practice, such properties are obtained through in-situ geotechnical testings or element soil testings in the laboratory. These methods require almost advanced equipment and both testing and sample preparation may be expensive and time-consuming. This paper aims to present an algorithm for a laboratory-scale ultrasonic non-destructive testing to determine the physical and mechanical properties of saturated soil samples based on the distribution of stress waves. The ultrasonic setup, in comparison to most conventional soil lab testing equipment, is low-cost and non-invasive such that it reduces the soil disturbance. For this purpose, a poro-elastodynamic forward solver and differential evolution global optimization algorithm were applied to characterize the porosity, density, and other mechanical properties for a soil column. The forward solver was developed based on a semi-analytical solution which does not require intensive computational efforts encountered in standard numerical techniques such as the finite element method. It was concluded that the proposed high-frequency ultrasonic technique characterizes effectively the saturated soil samples based on the output stress wave measured by the receiver. This development makes geotechnical investigations time-efficient and cost-effective, and as such more suited to applications in remote areas.

responses of geomaterials is neglected. In fact, the wave propagation in porous soil layers can be better represented by using dynamic poroelastic models instead of elastodynamic models, especially in fully saturated soils in which the pore water can significantly attenuate the stress waves, and in high frequency regimes. The dynamic poroelastic models consider the coupling effect between the pore water and solid skeleton, which induces three types of waves (fast P wave, slow P wave, and S wave in the solid skeleton). Under an impact load, those three waves travel at different speeds, which are captured by the receiver placed at the end of the soil specimen in an ultrasonic setup.
The problem of dynamic poroelasticity 10,11 has been solved using various analytical and numerical methods. A direct boundary element approach for solving three-dimensional problems of dynamic poroelasticity in the time domain was developed by 12 . Such a technique was based on an integral equation formulation in terms of solid displacements and fluid stress. The 2D and 3D fundamental solutions of dynamic poroelasticity was further developed by [13][14][15][16] . The solutions were obtained in both time and Laplace transform domain, and can be recovered to elastodynamics and steady-state poroelasticity. In layered saturated media, similar approaches have been reported by 17,18 . Other than the boundary element method, the finite element method has also been applied by 19 . The finite difference method is also used to simulate the wave propagation in heterogeneous poroelastic media by 20 .
In a conventional geotechnical apparatus used to determine the dynamic properties of a soil specimen, the focus is mainly on the estimation of shear wave velocity and the interpretation method is mostly based on the time interval difference between the input and output stress waves. To the best of our knowledge, there is currently no laboratory-scale ultrasonic setup which is able to determine a range of physical and mechanical properties of a soil sample. Furthermore, the development of cheaper, faster and portable means of soil characterization may significantly lower the cost of overall soil testing, and make better assessments possible in sensitive locations. This paper aims to present an ultrasonic-based poroelastodynamic algorithm, which can be used in an ultrasonic setup to determine a range of physical and mechanical properties of a soil sample such as shear wave velocity, compression wave velocity, density and porosity. Such a setup can also be used for in-situ geotechnical investigation on extracted soil samples. In this algorithm, the poro-elastodynamic forward solver for the characterization of soil samples in high frequency regimes is developed using the spectral element method. Such a meshless semi-analytical technique reduces significantly the computational efforts by avoiding unnecessary calculations for the entire domain. Instead, only the response at the receiver location is calculated, which will then be used during the optimization process. A robust global optimization algorithm is then applied to predict the soil properties given the stress signal measured by the receiver.

problem Statement
A general schematic of the problem is illustrated in Fig. 1. The domain is composed of a saturated porous medium. The transmitter located at one end of the sample generates the stress waves which travel through the specimen and is received by a receiver at the other end of the sample. The soil properties (Young's modulus, Poisson's ratio, density and porosity) will be captured by the proposed solver using the distribution of transmitted stress waves.

Dynamic poroelastic forward Solver
By assuming the infinitesimal deformation of solid skeleton, the dynamic poroelastic governing equations are written as follows:¨ü where u is the displacement vector of the solid skeleton; w is the fluid displacement relative to the solid skeleton; λ and μ are the Lamé constants; α is the Biot coefficient; p is the pore-water pressure; M is 1/ in which K f is the bulk modulus of the fluid; K s is the bulk modulus of the solid skeleton and φ is the porosity. λ c = λ + α 2 M; m = ρ f β∕φ in which β is the tortuosity which is used to describe the diffusion properties in porous media, and ρ f is the density of pore-water, taken as 1000 kg∕m 3 . The drag-force damping coefficient b is calculated as 21 : where η is the fluid dynamic viscosity and κ is the permeability coefficient; F is the viscous correction factor 22 : s c c f in which M s is taken as 1; = − i 1 and ω is the angular frequency. The governing equations can be written in frequency domain through the Fourier transform by performing convolution with e −iωt in which = − i 1; ω is the frequency and t denotes time variable. The governing equations in Laplace domain can be obtained by replacing ω with − is where s is the Laplace variable.
To obtain the analytical solution, the Helmholtz decomposition is used to decouple the P and S waves. The displacement vector is usually expressed in terms of a scalar potential (φ) and a vector potential (ψ ψ ψ ψ → = θ [ , , ] r z ), as shown in Eqs. (4a and 4b). In axisymmetric conditions, only the components in r and z directions are considered. Since P wave exits in solid skeleton and fluid, two P wave potentials are used, φ s and φ f , respectively.
The governing equations in frequency domain in terms of potentials are finally obtained as shown in Eqs. (5a, 5b, 5c and 5d): where ρ m = m − ib∕ω;  represents the terms in frequency domain.
Solution of dilation wave (p waves) using eigen decomposition. The equations in terms of P wave potentials (Eqs. (5a) and (5b)) in a matrix form is shown as: It can be seen from Eq. (6) that φ  s and φ  f are coupled in the governing equations. The diagonalization of such a matrix is required to decouple the system. The Eq. (6) is then rearranged into: , ( 2 ) The K matrix can be rewritten using the Eigen decomposition method: www.nature.com/scientificreports www.nature.com/scientificreports/ 1 where P is the eigenvector matrix and D is the eigenvalue matrix of the K matrix: It should be noted that Eq. (8) is still valid after neglecting the term k 1 21 in the eigenvector matrix P due to the existence of the term P −1 . Introducing Eq. (8) into Eq. (7) and by multiplying P −1 and P in the left and right sides, respectively, we can obtain: p p 1 2 , the system is finally decoupled as: 2 Under axisymmetric conditions, Eq. (10) for in cylindrical coordinates is written as: are a function of r and z in the cylindrical coordinates, the separation can be used. By setting the both sides equal to − k 2 where k is the wavenumber in the radial direction, we can obtain the following equations: The solutions to Eqs. (12a and 12b) are: k D z 2 2 11 in which J 0 is the Bessel function of the first kind; C 1 and C 2 are the coefficients to be determined from the boundary conditions. Similarly, the solution for φ  p1 can be obtained. The solution for where A and B are the coefficients to be determined from the boundary conditions. For simplicity, the term + k D 2 11 and + k D 2 22 is denoted as k p1 and k p2 , respectively.
Py , the solution for φ  s and φ  f can be finally obtained as: www.nature.com/scientificreports www.nature.com/scientificreports/ Solution of rotational wave (S wave). The rotational wave is governed by Eqs. (5c) and (5d). By replacing Under axisymmetric conditions, the solution for Eq. (16) in the cylindrical coordinates is obtained as: where C is the coefficient to be determined from the boundary conditions and J 1 is the Bessel function of the first kind of order one. For simplicity, the term + is denoted as k s .
Displacement, stress and pore-water pressure in terms of potentials. In the cylindrical coordinates (r, θ, z), considering the axisymmetric conditions = θ ∂ ∂ ( ) 0 , the vector potential ψ →  has only the component in the θ direction that does not vanish. For simplicity, the vector potential ψ →  in the θ direction is denoted as φ  s and φ  f for solid skeleton and porewater, respectively. This property reduces the displacement to the following forms: The effective stress and pore-water pressure are written as: Spectral element formulation for dynamic poroelasticity. In u-w formulation (displacement of solid and relative displacement of porewater), the displacement components w r and w z are linearly dependent. In this paper, only w z is used in the stiffness matrix. For two-node elements where a layer has a finite thickness, the matrix for the displacement components are written as follows:  www.nature.com/scientificreports www.nature.com/scientificreports/ Similarly, the matrix for effective stress components and porewater pressure in frequency domain is shown in Eq. 21 in which the components for matrix ′ S 2 can be found in Appendix A.  1 will be denoted as the G i matrix, in which i denotes the layer number.
After obtaining the stiffness matrix for each element, the global stiffness matrix can be obtained by applying the continuity conditions between the layer interfaces. The stiffness assembling method is shown in Fig. 2. The global stiffness is denoted as H matrix for simplicity. An example of the global stiffness matrix for a two layer system is provided in Appendix B.

Soil response under dynamic load (boundary conditions).
In the ultrasonic tests, a vertical impulse load f(t, r) is applied to one end of the soil specimen. The surface is assumed to be permeable, which implies the porewater pressure at the surface is zero. Under such conditions, the displacements in the frequency domain can be written as: The impulse load f is firstly defined in time domain and can decomposed into two independent functions in terms of time variable f n (t) and radial variable f r (r): n r The mathematical expression for the function f n (t) depends mainly on the type of impulse loads created by the signal generator. In this paper, a sinusoidal impulse function is used as the external load to simulate the applied load. The load with amplitude of one is mathematically described in Eq. (26).  where r 0 is the radius of the contact area; k m is the mode number; n is the total mode number; r ∞ is the diameter of the soil specimen. The displacement obtained in Eq. (24) is in the frequency domain. To obtain the soil response in time domain, the numerical Durbin inverse transform method is applied 23 :

Results and Discussion
The characterization of porosity has been a challenge because soil porosity can not be captured through traditional low-frequency tests. Such limitations can be explained by comparing the size of pore space and wavelength. A sensitivity analysis of the soil porosity is performed to verify such limitations. In this study, a soil column with a height and radius of 0.1 m is studied. The impulse load is applied to an area with a radius of 1cm at the center of the top end of the soil column. The displacement at the center (r = 0) in the other end is recorded and compared.
The typical values of Young's modulus, porosity, density, permeability and Poisson's ratio are well documented in the literature [24][25][26][27] . For example, high-plasticity clay (CH based on the Unified Soil Classification System (USCS)) has a Young's modulus ranging from 0.35 to 32 MPa and porosity from 0.39 to 0.59; Silts and clays of low plasticity (ML, CL) have a typical value of Young's modulus ranging from 1.5 to 60 MPa and porosity from 0.29 to 0.56; poorly graded sands (SP) normally have a Young's modulus from 10 to 80 MPa and porosity from 0.23 to 0.43; The Young's modulus of well-graded gravel (GW) is between 30-320 MPa and its porosity is from 0.21 to 0.32. The average dry density ranges from 1700 to 2300 kg∕m 3 . The average permeability varies from 5 × 10 −10 (clay of high plasticity) to 0.4 m/s (sand and gravel). The typical values of Poisson's ratio vary from 0.1 to 0.49 for clay and from 0.3 to 0.35 for silt.
In this paper, two groups of soils are studied: the first group includes clay, silt, sand and loose gravel which generally have a relatively low Young's modulus (lower than 100 MPa). The second group includes dense gravel which has a Young's modulus equal or greater than 200 MPa.
The effect of frequency and soil parameters on dynamic response. The effect of impulse load frequency and soil parameters on the dynamic soil response is studied in this section for the above-mentioned groups of soils. For the first group, the soil properties are taken as: Young's modulus is 20 MPa; Poisson's ratio is 0.35; dry density is 1800 kg∕m 3 . The wavelength can be calculated using the algorithm shown in Appendix C. Several sensitivity analyses under three impulse loads with various predominant frequencies are performed. The impulse load distributions in time and frequency domains are shown in Fig. 3. The loads 1, 2 and 3 have a predominant frequency of 0.05, 0.5 and 5 kHz, respectively. The amplitude of the input force is assumed to be 1 kN. The corresponding soil response at the receiver location is shown in Fig. 4.
As shown in Fig. 4, the different porosities (0.2, 0.4 and 0.6) give similar output displacement for load 1 and 2, which verifies that the size of pore space is not captured by the low-frequency impulse loads. In the inversion process, the porosity will be located at the shallow dimension, which makes the optimization algorithm difficult to be updated. Therefore, the characterization of saturated soil under low-frequency impulse load (below 5 kHz in this case) is nearly impossible. However, in the case of load 3 with a predominant frequency around 5 kHz, the effect of porosity is clearly triggered. The pore-scale of sand, for example, is around 760 μm as reported by 28 . Through the root search algorithm described in Appendix C, the wavelength under the load 3 is calculated around 1000-2000 μm, which is close to the poro-space scale of the studied soil. Therefore, the impulse load 3 is a good choice for the lab-scale characterization of soil specimens for group 1. www.nature.com/scientificreports www.nature.com/scientificreports/ Similarly, the sensitivity analyses are performed by considering different densities, Young's modulus and Poisson's ratios. The output displacement is shown in Fig. 5. The effects of Young's modulus, Poisson's ratio and density of soil are also shown in Fig. 5. A higher Young's modulus leads to a faster wave travelling speed and a smaller amplitude of the output wave. A higher density, on the contrary, leads to a lower travelling wave speed. Poisson's ratio that measures the tendency of material to expand in directions perpendicular to the direction of compression has an inverse relation with the wave speed. Therefore, it can be seen that the distribution of the output stress wave is a function of porosity, density, Young's modulus and Poisson's ratio.
In the case of soil group 2, dense gravel whose Young's modulus is up to 320 MPa, it is found that the load 3 (up to 5kHz) generates similar displacement outputs at different porosities (0.1, 0.3 and 0.5), as shown in Fig. 6. It means that load 3 can not trigger the effect of porosity. In order to characterize the porosity for very dense soils, one of the techniques is to further reduce the wavelength of the stress wave by increasing the frequency of the impulse load. It is found that an impulse load 4 with a higher predominant frequency (up to 0.5 MHz), as shown in Fig. 7, can effectively differentiate dense soils with various porosities. case study. In this section, a case study is presented to show the process of saturated soil characterization. For this purpose, a synthetic data is firstly generated to simulate real measurements. For simplicity, the results are only presented for soil group 1. The nature of this inversion problem and inversion algorithm selection are discussed in detail in the following sections. At the end, the inversion results (soil parameters) are given based on the synthetic data and selected inversion algorithm. www.nature.com/scientificreports www.nature.com/scientificreports/ Synthetic data. A synthetic data (the displacement measured by a piezoelectric receiver) is firstly obtained using the following settings: Young's modulus is 20 MPa; Poisson's ratio is 0.35; density of solid skeleton is 1800 kg/m 3 and porosity is taken as 0.3; The time interval is set to be 2 ms. Under the impulse load 3, as shown in Fig. 3, the snap shot of displacement contours (symmetric) at various time spans are shown in Fig. 8. The locations of impulse load and receiver are shown in Fig. 8. It is shown that the stress wave propagates through the sample and reaches the receiver at about 0.6 ms. The wave reflection at the bottom boundary is clearly visualized at time 0.8 ms and 0.9 ms.
The response measured at the receiver location is summarised in Fig. 9. In the laboratory ultrasonic test, the soil response is only recorded at the receiver location. Thus, in the following inversion process, only the results at the receiver location will be used as the input instead of the displacement at the entire domain.
Inversion algorithm. The inversion algorithm takes the measured displacement at the receiver location (shown in Fig. 9) as the input. The goal of the inversion process is to predict the soil properties including Young's modulus, Poisson's ratio, density and porosity based on the receiver signals. Given the initial guesses for the soil parameters, the inversion algorithm updates the prediction based on the difference between the displacement measured by the receiver and the predicted displacement response.
The update process can be achieved through the gradient-based and gradient-free optimization method. The gradient-based optimization is efficient in large convex problems such as linear least square problems and are commonly used in large optimization problems (e.g. deep learning and adjoint method). Therefore, the gradient based method is preferred in most cases, especially for convex optimization problems. However, such a method is highly likely to be affected by the local minimum since the gradient at any local minimum is zero. Thus, it is not favorable for non-convex problems.
An analysis was performed to show the nature of the soil characterization optimization problem. It is important to determine whether such application belongs to convex or non-convex problem. Then the corresponding optimization algorithm can be selected based on the nature of the problem. The aim (cost) function is defined as the Euclidean norm between the synthetic and predicted data. The optimization space can be visualized by performing parameter sweep. For example, the optimization space for the porosity and Poisson's ratio is shown in Fig. 10.
It is shown in Fig. 10 that a multiple local minimum exists in the optimization space. Therefore, the characterization of soil parameters is a non-convex optimization problem. If the gradient-based optimization algorithm is applied, the predictions will be highly dependent on the initial guess, which may leads to erroneous predictions in most cases. To make the estimation robust and accurate, a global optimization algorithm is favorable. In this work, the differential evolution algorithm that is designed for nonlinear and non-differential problems is used. Such an algorithm requires fewer control variables in comparison to other algorithms (e.g. genetic algorithm) and can be easily implemented in parallel computation 29 . www.nature.com/scientificreports www.nature.com/scientificreports/ A brief description of the differential evolution algorithm is given in Fig. 11. First, a population of candidate solutions are generated randomly; Then by moving around in the search space through a combination of the existing temporary solutions, a series of better solutions is expected to be obtained. In the differential evolution, the mutation constant is taken in the range of 0.5 to 1 and the recombination constant is recommended to be 0.9 30 .
Inversion results. Combining the synthetic data (as the input) shown in Fig. 9 and the differential evolution algorithm described above, the updates of the soil parameters and the corresponding values of the cost function are shown in Fig. 12. The iteration number shows the number of times that the forward problem is solved independently. After 200 iterations, the differential evolution algorithm stabilizes. The predicted soil parameters are as follows: Young's modulus is 20 MPa; Poisson's ratio is 0.35; density is 1800 kg∕m 3 ; porosity is 0.3 and loss function is 0. It can be seen that the prediction of soil parameters based on the transmitted wave measured by the receiver (as shown in Fig. 9) is exactly the same as the original input.  www.nature.com/scientificreports www.nature.com/scientificreports/ The differential evolution algorithm successfully finds the global minimum, despite of the existence of multiple local minimum. The spatial distribution of soil parameters updates are shown in Figs. 13 and 14. Through the projection of each parameter, it can be seen that Young's modulus is relatively easier to update. For the other three parameters (Poisson's ratio, density and porosity), there are multiple locations where cost function is close to zero. Thus, it took more number of iterations to update to the true values. However, it can be seen such a multidimensional optimization problem is well handled by the differential evolution algorithm.
Uncertainty analysis. The predicted soil properties (Young's modulus, Poisson's ratio, density and porosity) are likely to be affected by the noise level of the measurement data, which could be introduced by the sensor measurement errors and ambient noise. In this uncertainty analysis, random white noise is added to measured displacement data with targeted signal-to-noise (SRN) ratio. For example, the noisy data with 10 and 20 dB of SRN is shown in Fig. 15a. A normal distributed probability density function of SRN is used as the input to account the uncertainty introduced by noise, as shown in Fig. 15b. It is assumed that there is a 28% possibility to have a SRN of 20 dB in measured data.  www.nature.com/scientificreports www.nature.com/scientificreports/ In addition, the uncertainty can be introduced by the unknown coupling performance in the interface of piezoelectric sensors and soil specimens. The input electricity signal does not necessarily generate the desired input pressure. To account for such uncertainties, the magnitude of input load is assumed to be in normal distribution,   www.nature.com/scientificreports www.nature.com/scientificreports/ as shown in Fig. 16a. The uncertainty also comes from the inherent soil property assumptions made in soil specimen during the inversion analysis, such as hydraulic conductivity. Thus, a normal probability distribution is also applied to account such uncertainty, as shown in Fig. 16b.
The generalized Polynomial Chaos Expansions (PCE) method developed by 31 is used for the uncertainty analysis in this paper. The PCE technique, as a rigorous uncertainty quantification method, provides reliable numerical estimates of uncertain physical quantities. It was also reported that the PCE is much faster than Monte Carlo methods when the number of uncertainty parameters are lower than 20 32 . The 90% confident interval of the displacement at the receiver location is calculated through the PCE technique, shown in Fig. 17.
Then, based on the inversion analysis, the predicted soil properties in the 90% confidence interval are shown in Table 1. Then, the variation ratio is calculated by comparing the mean values (obtained through uncertainty analysis) with the original predictions. It is found the prediction of porosity could be affected by the uncertainty introduced by the white Gaussian noise, coupling effect between transmitter and soil specimen as well as other factors. However, various signal processing methods can be used to improve the noisy measurements.

conclusions
In this paper, an ultrasonic-based characterization of soil specimens is developed for the instant measurement of soil properties including Young's modulus and Poisson's ratio (compression/shear wave velocity), density and porosity. The developed meshless semi-analytical algorithm reduces the computational effort significantly in comparison to standard numerical techniques such as the finite element method. In fact, the advantage of such a solution is that the dynamic response is evaluated at the receiver location only rather than the entire domain. The soil response in other locations is not measured in the real application and does not factor in soil characterization.
It is concluded that high-frequency impulse loads (with predominant frequency of up to 5 kHz) is required to trigger the effect of porosity for soils with relatively low Young's modulus (e.g clay, silt and sand). For stiffer  www.nature.com/scientificreports www.nature.com/scientificreports/ materials, such as very dense gravels, an impulse load with predominant frequency of 0.5 MHz is required to characterize their porous nature. The characterization of soil properties has been proved as a highly non-convex optimization problem in this paper. The differential evolution algorithm, as a global optimization method, is found efficient and effective in finding the optimum soil properties, such that the difference between the predicted Figure 14. Updates of Young's modulus and density through a differential evolution algorithm.  www.nature.com/scientificreports www.nature.com/scientificreports/ and measured stress waves is minimized. In conclusion, the developed method in interpreting dynamic response of saturated soil can be used for the immediate characterization of Young's modulus, Poisson's ratio, density and porosity for a given soil specimen.

Appendix A: components of Matrix
The components of the matrix ′ S 2 for effective stress components and porewater pressure in frequency domain is shown as follows: