Revisiting Piezoelectric Sensor Calibration Methods Using Elastodynamic Body Waves

The application of absolutely calibrated piezoelectric (PZT) sensors is increasingly used to help interpret the information carried by radiated elastic waves of laboratory/in situs acoustic emissions (AEs) in nondestructive evaluation. In this paper, we present the methodology based on the finite element method (FEM) to characterize PZT sensors. The FEM-based modelling tool is used to numerically compute the true Green’s function between a ball impact source and an array of PZT sensors to map active source to theoretical ground motion. Physical-based boundary conditions are adopted to better constrain the problem of body wave propagation, reflection and transmission in/on the elastic medium. The modelling methodology is first validated against the reference approach (generalized ray theory) and is then extended down to 1 kHz where body wave reflection and transmission along different types of boundaries are explored. We find the Green’s functions calculated using physical-based boundaries have distinct differences between commonly employed idealized boundary conditions, especially around the anti-resonant and resonant frequencies. Unlike traditional methods that use singular ball drops, we find that each ball drop is only partially reliable over specific frequency bands. We demonstrate, by adding spectral constraints, that the individual instrumental responses are accurately cropped and linked together over 1 kHz to 1 MHz after which they overlap with little amplitude shift. This study finds that ball impacts with a broad range of diameters as well as the corresponding valid frequency bandwidth, are necessary to characterize broadband PZT sensors from 1 kHz to 1 MHz.


Introduction
As brittle materials are subjected to external stress in a laboratory setting, localized and rapid inelastic deformation events occur that are associated with the growth or appearance of small defects at the grain-scale (from micrometers to millimeters), which can generate acoustic emissions (AEs) [1][2][3]. These emissions can cause high-frequency vibrations, at frequencies ranging from tens of kHz to several MHz, and are recorded by piezoelectric (PZT) sensors at known locations. In the last decade, great effort has been made to improve the understanding of laboratory generated AEs in a quantitative manner [4][5][6][7][8][9].
Fracture induced stress waves in rock formations at scale are also referred to as AEs that carry frequencies ranging from hundreds of Hz to tens of kHz. AE monitoring is utilized in industrial activities to manage the induced seismic risk, such as hydrocarbon storage and extraction, shale gas exploitation, mining operations, and impoundment of water behind dams [10]. To quantitatively study induced seismicity, i.e. AEs caused by anthropogenic forcing, local down borehole networks of PZT sensors were developed at scale from decameter to hectometer [11][12][13][14][15]. Kwiatek et al. [11] reported on picoseismic and nanoseismic events (M w > −4.1) caused by aftershock sequence following the postblasting activity inside a volume of 300 × 300 × 300 m at Mponeng Deep Gold Mine, South Africa. They installed PZT sensors that were relatively calibrated according to 3component accelerometers at frequency from 400 Hz to 17 kHz in boreholes along the tunnel and monitored the seis-mic events in a magnitude range from M w = −0.8 down to M w = −4.1 with corner frequencies from 800 Hz to 13.6 kHz (source length scales ranging from 8 cm to 1.3 m).
These studies differ from previous AE studies, in that they characterized the absolute mechanical energy released by fracturing processes due to the radiated waves in the stressed solids instead of the traditional parametric analysis of the voltage measured by PZT sensors [16,17].
Recent advances in both laboratory and in situ AE monitorings during fracturing experiments has greatly improved our understanding of microcrack mechanisms over broadband ranges of source dimension (micrometers to meters) and frequency (hundreds of Hz to several MHz). Prior to exploring such seismic characteristics, it is essential to absolutely characterize or calibrate the PZT sensors utilized in both in situ and laboratory applications so that the information conveyed via the ground motions can be interpreted from the measured voltages [2,[18][19][20].
Green's functions, used to map active source to theoretical ground motion, are vitally important for PZT sensor characterization. Researcher [21] investigated elastic stress wave propagation within a semi-infinite homogeneous and isotropic elastic plate, which was first solved [22] and known as "Lamb's problem". More specifically, Lamb's problem focuses on calculating the elastic disturbance caused by stress waves due to a point force in/on a half space. To find the solution of Lamb's problem (or "Green's function"), researchers numerically solved Lamb's problem starting from generalized ray theory [21,[23][24][25].
Two main concerns limit the application of the generalized ray theory to calculate Green's functions. First, the corner frequency of amplitude spectra of in situ AE events could be as low as hundreds of Hz to several kHz [11,15]. To model this, we require the spectra of Green's functions at least down to 1 kHz to calibrate PZT sensors of the same frequency band. This requires large computational loads to obtain the huge number of possible ray paths of Green's functions. Second, sample finiteness makes the semi-infinite conditions associated with Lamb's problem unrealistic for laboratory investigations and, therefore, the ray paths of side reflections from a finite elastic plate are non-negligible.
The finite element method (FEM) is an alternative approach to obtain the numerical Green's functions with regard to elastodynamic wave propagation [26][27][28] where more realistic boundary conditions, similar to that of an experimental configuration, can be modelled. Few researchers [29][30][31] performed FEM analysis with idealized or simplified boundary conditions to study the instrumental response of PZT sensors utilized in laboratory experiments. They obtained the ground motion down to 1 kHz at the location of receivers created by ball impacts of various sizes on the top surface of an elastic plate. Their approach by using conventional FEM codes and fine grids, however, requires huge computational costs to obtain relatively flat amplitude spectral of numerical Green's functions from 1 kHz to 1 MHz where the specified spectral resolution needs to be satisfied.
In this study, we used a state-of-art FEM-based solver to numerically calculate true Green's functions from 1 kHz to 1 MHz excited by a Heaviside step force-time function. To improve the computational efficiency, we perform FEM modelling using a fine grid to compute Green's functions of high frequency, and using a coarse grid to obtain Green's functions of low frequency. In course of modelling the low-frequency ground motions from 1 kHz to 100 kHz, physical-based boundary conditions are utilized in the FEM modelling to mimic the realistic experiments: elastic body wave propagation, reflection and transmission in/on an elastic medium. The numerical Green's functions of a group of source-receiver pairs from 1 kHz to 1 MHz are obtained and then the corresponding displacement at the location of the receivers is derived.
We performed ball impact tests over a range of diameters and used PZT sensors to measure the vibration in terms of voltage signal. Our analysis characterizes the broadband instrumental responses through accurately cropping (primary, secondary and tertiary) and linking (primary and secondary) a group of segmented instrument responses over valid frequency bands in line with the accuracy of FEM solutions and quality of experimental data. Moreover, we extend the instrumental response analysis from a single sensor to an array of PZT sensors at take-off angles taking account into their unique source-receiver characteristics. Figure 1 shows six concepts from a source to the output voltage that are used in our analysis: (1) active source, (2) Green's function , (3) theoretical displacement, (4) instrumental response, (5) amplification and bandpass, and (6) voltage signal. An active source is used to produce a rapid transient force f j (ξ , τ ), in the j direction at point ξ and delayed time τ on the top surface of medium (Ω). Elastic waves propagate through the medium and the theoretical displacement is represented as u k (x, t) in the k direction throughout the medium (Ω) at any location x and time t.

Instrumental Response of PZT Sensors
Since elastic wave propagation has linear time-invariant characteristics and due to the spatial reciprocity of the representation theorem [32], u k (x, t) can be expressed as wave propagation A D contact converter sensor cable Fig. 1 Schematic drawing of concepts to illustrate the operation principle of PZT sensor that link the active source to generated voltage signals measured using a data acquisition system where * denotes the convolution operator, g k j (x, t; ξ , τ ) denotes the Green's function between the location of the source and receivers, subscripts k,j reflect Einstein summation convention. Note that the assumption of point representation at the contact region of both the active source and PZT sensor is used, which is not exactly true due to finite area of the active source and the aperture area of the PZT sensor. This effect can be minimized by using a ball impact source and a conical-frustum PZT crystal with minimal contact area; this is discussed further in Sect. 2.2.
The time-varying displacement u k (x, t) measured by the PZT sensor is then distorted into a voltage signal ψ(x, t), which is recorded by a connected data acquisition (DAQ) system. The instrumental response, i k (t), maps the true mechanical input u k in the k direction to the measured voltage ψ(x, t). The mapping is assumed to satisfy a linear timeinvariant system such that where a(t) is the documented response of used amplifier including signal amplification, bandpass filter, etc. We aim to quantify the instrumental response of the PZT sensor, which integrates the effects of contact manner, intrinsic sensor characteristics, cables and analog/digital (A/D) converter. Note that if a(t) is not available through experiments, the effect of a(t) can be superposed into the final instrumental response. We perform the deconvolution operation to Eq. (2) in the frequency domain to compute I k (ω) as where ω denotes the ordinary frequency. I k (ω), Ψ (x, ω), U k (ω) and A(ω) are variables in the frequency domain corresponding to i k (t), ψ(x, t), u k (x, t) and a(t), respectively.
To obtain an accurate instrumental response I k (ω), we present the detailed analysis of acquiring voltage Ψ (x, ω) in laboratory experiments (see Sect. 2.2) and numerically calculating Green's function g k j (x, t; ξ , τ ) (see Sect. 2.3).
A(ω) used in this study is provided from [33]. The force-time function is determined from Hertzian impact theory [34] (see Appendix A).

Experimental Setup: PZT Sensor Calibration Station
In this section, we introduce the laboratory experiments performed on a sensor calibration station (see Fig. 2a). We characterize the performance of PZT sensors to measured kinematic motion excited by body wave propagation through an elastic isotropic, homogeneous steel transfer plate due to the active source produced by a steel ball impact. A detailed comparison of characteristics, (dis)advantages and trade offs of laboratory active sources was explained in [35]. An extruded aluminum framework is used to support an overlaid square steel transfer plate (35 cm × 35 cm × 5 cm). Pads (4 yellow patches, 3 cm × 3 cm × 0.5 cm) made of nitrile butadiene rubber (NBR), with low mechanical impedance (∼ 1/11 of steel), are used to effectively block elastic wave refraction into the aluminum framework. An electromagnetic holder is built in the upper crossbeam of the aluminum framework. Once the power is turned off, the steel ball (green) is released and drops down freely until impact at the central location of the top surface of the steel transfer plate.
In Fig. 2b, we show 13 PZT sensors with the same KRN Services (model KRNBB-PC) mounted at the bottom surface of the steel transfer plate. The elastic waves generated by ball impact propagated throughout the steel transfer plate and were recorded at the opposite side of the steel plate at various incident angles. Due to the relatively large thickness of plate and the locations of source-receiver pairs, we highlight that the proposed calibration scheme as a "body wave calibration" in which body waves are primarily investigated. The array of PZT sensors has unique source-receiver characteristics that need to be accounted for [36,37]. The sensor mounting plate has 7 × 7 locations with 23 mm spacing in both the 1 and 2 directions, respectively. Since body waves caused by ball impact in the 3-direction are symmetric about the 1-3 and 2-3 planes, we have focused most of our PZT sensors converge to one quadrant of the sensors mounting plate. In Fig. 2b, similar colors for the PZT sensor locations (A1 to D1) represent the redundant epicentral locations with respect to the source. The darker to lighter colors represent the increasing epicentral locations. We use the take-off angle θ to characterize 7 possible seismic ray paths between the source and PZT sensors, θ = 0 • (A1); 24.7 • (A2, B1); 33.0 • (B4); 42.6 • (A3, B2, C3, C4); 52.4 • (C1, D1); 54.0 • (A3, B4); 62.9 • (C2).
KRNBB sensor used a conical-frustum PZT crystal with an aperture diameter of 1.5 mm [38]. We assume that this aperture roughly approximates the tip contact area between the circular tip of the KRNBB-PC sensors and the lower  surface of the steel transfer plate, which is small relative to the dimension of the steel transfer plate. Also, the contact area between the ball and the upper surface of the steel transfer plate at the time of impact is negligible. Thus we assume that the point characteristics of the source and receiver in Sect. 2.1 is satisfied. No couplant (e.g. hot glue, Vaseline) is used through this study. These sensors are sensitive to ground motion in the 3-direction over a wide frequency range (1 kHz to 1 MHz) and their spectral characteristics have been well documented [9,25,38,39]. A small account of historical development of the KRNBB sensor is given in Appendix B.
To acquire voltage ψ(x, t), a DAQ system (Elsys Instruments TraNET, 32 Channel TPCE-2016-4/8) is connected with the PZT sensors with a sampling frequency F s of 20 MHz per channel. The Nyquist frequency is 10 MHz so that the F s is adequate enough to perform sampling. Elsys AE amplifiers (AE-Amp) were used to provide the internal amplification [JFET, see 38] to sensors with 25 mA excitation and the gain was 0 dB.
To obtain instrumental responses, we perform Fast Fourier Transform (FFT) of the measured signals [40] over the frequency band from 1/T w to F s /2, where T w is the time window and F s is the sampling frequency. In this case, when the frequency approaches the lower bound 1/T w , low resolution of amplitude spectral can occur if not enough lowfrequency cycles are analyzed. A proper time window T w is critical to accommodate trade offs between spectral resolution and computational costs in modelling the excited elastic waves. We suggest that there should be at least n = 8 bins from the lowest frequency limit ( f min = 1 kHz) to its adjacent frequency 2 kHz. Since the linearly spaced increment of frequency bins d f equals 1/T w , inequations can be given as where a window length of T w = 8000µs is used to crop caused voltage ψ(x, t) and displacements u 3 (x, t) with an identical length. To avoid spectral leakage, a symmetric window function with the same length as ψ(x, t) and u 3 (x, t) is essential. In Fig. 3, we show the raw data (gray line) of ψ(x, t) for typical acoustic events due to ball impacts with diameters of 0.3 and 3 mm. A 8000-µs Blackman-Harris window function win(t) (centered about the first P-wave arrival) is used throughout this study to obtain windowed signals (red and blue lines). The windowed signals have the maximum in the middle, and taper away from the middle. By performing FFT to the windowed ψ(x, t) and u 3 (x, t), and neglecting the phase information, we obtain where |U 3 (x, ω)| and |Ψ (x, ω)| are the amplitude spectra of windowed ψ(x, t) and u 3 (x, t), respectively. | · · · | represents the absolute value operator. g 33   Raw data (gray) and windowed data (colored) centered at the first P-wave arrival from typical acoustic events due to ball impact, with diameters of 0.3 (red) and 3 mm (blue), measured by the PZT sensor A1 located directly beneath the impact point. Window functions with two lengths 80 µs (see inset) and 8000 µs are used for high-and low-frequency analysis in Appendix C, respectively. where · denotes dot product, and F{} represents the FFT operation. Material parameters (e.g. elastic properties, density, geometry) of the steel ball, steel transfer plate and pad used in the calibration station are summarized in Table 1.

Green's Function
We now aim to obtain Green's function, g 33 (x, t; ξ , τ ), which reflects the displacement component in the 3-direction due to time-delayed Dirac delta (true impulse), δ(t − τ ) in the 3-direction. However, it is not easy to implement δ(t − τ ) numerically. Instead, we use the Heaviside step function, H (t − τ ), which is the integral of δ(t − τ ) with respect to time. The corresponding f 3 (ξ , τ ) can be expressed as where δ(x − ξ ) denotes the spatial source distribution of the Dirac delta function. In this study, δ(x − ξ ) is described as a limit representation of the Dirac delta where S is associated with the mesh size in the vicinity of the ideal loading point; its exact value is discussed in Sect. 3.2.
The u 3 (x, t) due to f 3 (ξ , τ ) is given by Recalling the linear time-invariant characteristics of convolution in Sect. 2.1 and the properties of convolution differentiation [40], we can conduct time differentiation oper- Therefore, Green's function g 33 (x, t) is derived from the particle motion velocity in the 3-direction, v 3 (x, t), caused by the force-time function H (t − τ ) with a spatial distribution described by Eq. (10).

Numerical Modelling of Green's Function
In this section, we present the FEM-based methodology to calculate numerical Green's functions. We validate the methodology against the reference approach over a highfrequency band (from 100 kHz to 1 MHz) and then extend it to low-frequency analysis down to 1 kHz.

Governing Equations
We simulate the elastic wave propagation using a stateof-art discontinuous Galerkin (DG) FEM by COMSOL Multiphysics [41]. The steel transfer plate and pads are modelled explicitly and their material properties are provided in Table 1. The particle motion velocity, v, and strain, E, with the implementation for steel (k = 1) and rubber (k = 2) domains, obey the first-order elastodynamic equations where ρ k stands for the material density, S denotes the Cauchy stress tensor, and f is the applied loading described by Eq. (9). C k represents the isotropic stiffness tensor characterized by Young's modulus, E, and Poisson's ratio, ν. To guarantee the uniqueness of v and E, the zero-traction condition is enforced on the free surface Σ F of the steel transfer plate, that is where n is the outward unit normal vector. By incorporating numerical flux, DG FEM weakly imposes mechanical continuity of v and E across the interior boundary, which is especially computationally efficient for studying 3Dtransient wave propagation problems [27].
To model elastic stress wave transmission and reflection along the interface, Σ RA (see Fig. 2a) between the pads and aluminum framework, we impose a velocity-dependent traction on Σ RA . This traction is caused by the non-zero mechanical impedance of the aluminum framework (k = 3) and can be written as a combined effect of P and S waves [42] where ρ k , c k p , c k s are the density, P-and S-wave velocity of aluminum, respectively, and t is the unit tangent vector. By adding this traction along Σ RA , aluminum framework is not explicitly modelled.

Modelling Parameters
In this section, we describe how the numerical calculations are performed to obtain numerical Green's functions, including the meshing schemes, time step and simulation procedures.
For meshing schemes, one of most critical issues is to choose an optimized mesh size; a relatively coarse grid works in a similar manner to a high-pass filter, leading to inaccurate solutions whereas a very fine grid results in huge computational consumption (in terms of memory and CPU time).
To produce a satisfactory solution, the maximum size h of the mesh elements should be lower than the wavelength of interest λ. Since higher-order Lagrange interpolation functions, up to the forth order, are utilized, this ratio could be set as a larger value [41]. In this study, the ratio is assumed to be 1, that is h λ = 1, where λ = c p ω 0 , λ, c p denote the wavelength and velocity of P wave, and ω 0 represents the upper limit of the bandwidth of interest. Even when using the efficient DG FEM, it is still time-consuming to calculate high-frequency numerical Green's functions. To reduce the computational cost, the adaptive grid scheme in FEM, can be adopted; however, this, is a non-trivial task. To remedy this, we introduce two separate simulations, which include a lowfrequency model with a relatively coarse grid (c-FEM), and a high-frequency model with a relatively fine grid (f-FEM).
Due to the geometric symmetry of the calibration station, only 1 8 of the steel transfer plate and pad was modelled. We set ω c 0 = 300 kHz (>100 kHz) and ω f 0 = 1.2 MHz (>1 MHz) to control h, where the superscripts c and f stand for c-FEM and f-FEM, respectively. The mesh around the center of loading area (red arrow) was locally refined to ensure a correctly applied f 3 (ξ , τ ). We determine the S value from Eq. (10) as 1e-6 (units: 1/m 2 ) by comparing the numerical integration of normal stress around the loading against f 3 (ξ , τ ). Grid discretization scheme of models simulated in c-FEM and f-FEM are shown in Fig. 4a and b, which have 762 and 104624 unstructured tetrahedral quartic elements, respectively. Details of material properties and the geometry information are given in Table 1.
The Courant-Friedrichs-Lewy condition is adopted to determine the stable time step, Δt, which is automatically optimized by COMSOL Multiphysics. During our simulations, the 4 th explicit Runge-Kutta method is adopted. Note that we use a "pseudo-Heaviside" loading scheme in the modelling: when we enforce a Heaviside step force-time function to the loading point; however, there is a very short rise time, τ r > 0, to reach the peak loading. To prescribe the rise time exactly, we specified Heaviside step functions at the same loading point (see red arrow in Fig. 4a) to both c-FEM and   f-FEM and probed the applied force history, respectively. We found that τ r is around 5 and 10 ns for c-FEM and f-FEM, respectively. Considering the fact that the force-time function of glass capillaries fracturing has a τ r of approximately 200 nanoseconds while most of its spectral energy does not fall below 2 MHz [9], this pseudo-Heaviside loading scheme has little effect on the spectrum over the bandwidth of interest. For modelling procedures, we simulate c-FEM for 4000 µs and f-FEM for 40 µs (half length of suggested time window in Eq. (5)) such that the spectral resolution is satisfied down to 1 kHz and 100 kHz, respectively. The computational efficiencies of one FEM simulation are summarized. To perform the simulation of c-FEM for 4000 µs, 6 cores (41 hours for each core) and average memory 4.6 GB are required. Meanwhile, regarding simulating f-FEM for 40 µs, 16 cores (13.1 hours for each core) and average memory 13.3 GB are needed.
We obtain time-varying v 3 (x, t) probed at the location of the virtual seismometer (red square in Fig. 4a) and derive the numerical Green's functions from both models by means of the transformation described by Eq. (12).

Model Validation: High-Frequency Analysis of Short-Duration Excitation
In above section, we propose using ω c 0 and ω f 0 to control the mesh size such that the frequency component of numerical Green's functions can be expected to be kept below 100 kHz and 1 MHz. Our aim is to validate the capabilities of these models and find the corresponding valid frequency band beyond which the numerical Green's functions deviate from reference results. Note that the reference ("true") results are computed by an approach based on the generalized ray theory (GRT ) [21,24] that governs transient wave propagation in a semi-infinite purely elastic plate. Figure 4c and d show the magnitude of the particle motion velocity (units: m/s) to qualitatively illustrate the performance of the f-FEM model in simulating elastic wave propagation. We see that the elastic wave initiates from the loading point and geometrically spreads. Figure 4c shows the first P-wave ray around the virtual seismometer and the first S-wave ray that follows at time 8.6 µs. The spatial distance between the wavefront of the first P and S waves is controlled by the wave propagation times and the speed difference between them.
In Fig. 4d, we observe multiple reflections of different rays of the propagated elastic waves between the upper and lower surface of the steel plate at time 25.8 µs: first P wave, first S wave, PPP wave, etc. The PPP-wave ray is shown around the virtual seismometer. These rays are important since they offer important information and frequency content in the velocity seismogram. We note that here we are only showing validation efforts for scenarios where the virtual seismometer is located directly below the source; however, this methodology can be extended to more source-receiver pairs due to their unique source-receiver characteristics [36] (see Sect. 4.3).
We now validate the numerical Green's functions at the virtual seismometer over a short duration. The seismometer starts calculating v 3 (x, t) from the instantaneous loading until the initial side reflection back. In our case, this duration is not allowed to exceed 60 µs so that numerical Green's function should, theoretically, be the same as that of the semiinfinite elastic plate since the wavefront has not reached the edges of the transfer plate. We simulated both c-FEM and f-FEM and the GRT for 40 µs under the loading described by Eq. (9). Figure 5a shows the amplitude spectra of Green's functions in the frequency domain. The amplitude spectrum calculated using the reference approach GRT (black line) is almost flat below 2 MHz since the GRT solution is semi-analytical. For numerical Green's functions from FEM modelling, the amplitude spectra match well with these of the GRT when the frequency ranges from 100 to 324 kHz (for c-FEM, blue line) and from 100 kHz to 1. Note that to increase the accuracy of FEM solutions, we consider that the Green's functions from the c-FEM model are valid up to 174.6 kHz instead of 324 kHz. After validating the capabilities of our built model with a simple geometry boundary, now we have a tool for studying elastic wave propagation that can accommodate realistic geometries with a validated level of accuracy. To have a better visualization in the accuracy of FEM solutions, Fig. 5b presents the theoretical displacement, u 3 (x, t), in the 3-direction from GRT and f-FEM due to 0.5 mm ball impact (other parameters refer to Table 1) in the time domain. The u 3 (x, t) from f-FEM is essentially equivalent to that of GRT. First P, S and reflected PPP phase of the elastic waves were captured and windowed inside the 3 pink bars.

Model Extension: Low-Frequency Analysis of Long-Duration Excitation Down to 1 kHz
We extend the capability of c-FEM to perform low-frequency analysis by elongating the simulation duration. Due to the finite dimension of the calibration station, the simple Lamb problem, where only the elastic wave reflection between the bottom and top surface of the steel transfer plate is modelled, is no longer valid. To study the effects from boundaries on the elastic wave propagation problem, we conducted the following three modelling scenarios: 1. float-NF: modelling the elastic wave propagation through an unsupported or floating steel transfer plate. We adopt similar boundary conditions (free of stress over all surfaces of the tested specimen) to those used in [31] and construct a float model. 2. float-HP: in the float model, rigid body motion occurs because there are no constraints from the supported material (i.e., the pad beneath the steel plate) are modelled.
To correct the numerical Green's functions imposed by improper boundary conditions, we process the raw v 33 (x, t) using a minimum-order high-pass filter. 3. true BC: to better model the real elastic wave propagation, more physical boundary conditions, described by Eqs. (14) and (15) were already introduced in the c-FEM model to simulate the elastic wave reflection and refraction occurring at the interfaces of the steel transfer plate, pads and aluminum framework.
We extract v 33 (x, t) with a duration of 8000 µs (centered at the first P-wave arrival) in the above scenarios and use the transformation (see Eq. (12)) to obtain the numerical Green's functions for the same source-receiver pair. By performing FFT and neglecting the phase information, we obtain the amplitude spectra of low-frequency Green's functions, which are termed as float-NF (blue line in Fig. 6a), float-HP (black line in Fig. 6a) and true BC (red line in Fig. 6a).
In Fig. 6a, float-NF tends to decrease 'linearly' from ∼1e-12 to ∼1e-14 m/N from 1 kHz and 100 kHz. The large amplitude of low-frequency components is caused by the rigid body motion of the steel transfer plate. By performing a high-pass filter operation on float-NF, the rigid body motion is significantly removed. We see that the float-HP has a relatively flat spectrum with few distinct spectral peaks representing the low-frequency resonances and anti-resonances of mechanical vibration. The resonant frequencies from float-HP and true BC match well over 1 to 100 kHz, demonstrating that a consideration of the true physical boundaries has little effect on the shift of the resonant frequency.
For a better comparison, we plot the ratio of true BC to float-HP from 1 to 100 kHz in Fig. 6b. From 1 to 10 kHz, the ratio at the locations of resonant and anti-resonant fre- quencies deviates significantly from 1. These locations are labelled AR (anti-resonance) and R (resonance). We find that around the resonant frequency (i.e., R1 and R2), the true BC is nearly half the float-HP. Conversely, around the antiresonant frequency, the true BC can be 7 to 13 times than of the float-HP. This ratio tends to be flat and close to 1 for frequencies between 10 and 100 kHz. This suggests that both scenarios have similar Green's functions over 10 to 100 kHz. We modelled how the physical constraints resulting from the experimental configuration affected the elastic wave propagation problem. We interpret the results from the physical-based boundary conditions. Kinematic energy of body waves initially flows from the loading point, trav-els through the steel transfer plate and pads, and finally leaves the c-FEM model naturally where the spatiotemporal evolution of the energy field is dominated by properties of both media (steel plate, pads) and their interfaces. No extra high-pass filter operation is needed to remove the low-frequency component of the kinematic energy, which is fully trapped inside the steel transfer plate in the float model. Therefore, the constructed c-FEM model can be utilized to obtain numerical Green's functions in a more physical way, which becomes important when more complex geometries and boundary conditions are applied. To solve the multiphysics problem of elastic wave propagation, our constructed models can potentially be integrated with multiphysics couplings (i.e. temperature and fluids) from the COMSOL Multiphysics software whose abilities have been validated against theoretical solutions and/or laboratory experiments in the geomechanics communities [43][44][45][46].

Broadband Characterization of PZT Sensors
In the above sections, we described the FEM-based methodology used to obtain high-frequency numerical Green's functions from 100 kHz to 1 MHz as well as low-frequency counterparts from 1 to 100 kHz. By performing ball impact experiments, we presented the displacement |U 3 | s,l χ and voltage spectra |ψ| s,l χ in short-(denoted by s, 100 kHz to 1 MHz) and long-duration (denoted by l, 1 kHz to 100 kHz) excitation in Appendix C. Label χ denotes an unique ball diameter, χ ∈ [1, N ], where N is the number of ball sizes. In this section, we obtain a group of segmented instrumental responses I s,l 3,χ with different forces and duration of ball impact in accordance with Eq. (8). We develop an algorithm of spectrum cropping and linking to integrate a group of I s,l 3,χ into a truly broadband understanding of I 3 from 1 kHz to 1 MHz for single sensor as well as an array of PZT sensors with unique source-receiver characteristics.

Algorithm to Integrate Segmented Instrumental Responses
We are reminded of the fact that a specific member of the I s,l 3,χ group is not valid over the whole frequency bandwidth. Instead, the accuracy of I s,l 3,χ is regulated by the quality of the experimental data (i.e., corner frequency ω c of voltage spectra |Ψ | s,l , signal-to-noise ratio or SNR, variations among repeated tests) and the solution accuracy of the FEM modelling. If we perform a simple union operation of the I s,l 3,χ group, the broadband I 3 will be distorted, with various level of uncertainty, by introducing all the components of I s,l 3,χ . To obtain a more accurate I 3 , we need to perform several 'cropping' operations to constrain the I s,l 3,χ group. The accuracy of FEM solutions is closely correlated with the mesh discretization. In high-frequency analysis, I s 3,χ is picked from 100 kHz (the lower bound of the valid frequency band for the f-FEM model) to the corner frequency ω c of voltage spectra |Ψ | s . In low-frequency analysis, I l 3,χ is cropped from 1 kHz to the minimum between 174.6 kHz (the upper bound of the valid frequency band of the c-FEM model) and ω c of |Ψ | l . This primary cropping operation can be denoted as / · · · / such that / I s,l 3,χ / is obtained. The quality of the experimental data is mainly controlled by the SNR (or ball size: large ball impacts generate higher SNR) and the repeatability of the test data. To alleviate the effect of background noise, a secondary cropping operation is proposed to better select the frequency band of I s,l 3,χ based on the criterion of SNR > 1. We denote this as · · · and obtain / I s,l 3,χ / . In this case, the small ball impact has a limited contribution to the broadband I 3 due to its low SNR especially at the low-frequency range. We have illustrated the repeatably of ball impact tests in Appendix C, but there exists an amplitude offset of I s,l 3,χ during repeated tests. We therefore remove the frequency band where the relative standard deviation of |Ψ | s,l is larger than specified threshold (2 % used in this study). We then average I s,l 3,χ over the rest of the frequency band to get / I 3,χ / ; this tertiary cropping operation as · · · is denoted.
Once the above 'cropping' operations are implemented, the I s 3,χ group is well constrained. We now need to link all the cropped I s 3,χ together from the low-and high-frequency analyses over a broad range of ball diameters. The final broadband I 3 is written as where ∪ is the primary link operation between the low-and high-frequency cases, denotes the secondary link operation over different ball diameters. To help understanding the used algorithm, we present the general principle and flowchart of cropping and linking operations in Appendix D.
In Fig. 7, we present the caused seismic magnitude M w versus valid frequency band of the grouped ball impacts. Note that M w due to external ball impact is estimated from [47]. The region between ω min χ (black dashed line) and ω max χ (black solid line) suggests that these overlapped frequency bands of ball drops could fully cover the frequency scope of interest (from 1 kHz to 1 MHz). Here ω min χ and ω max χ represent the minimum and maximum valid frequencies of ball drops, respectively. We see that the results from ball impact tests with diameters of 0.4, 1 and 3 mm almost cover broadband frequency range. This provides researchers with straightforward instructions about how to choose the right parameters for ball impact tests in accordance with their frequency range of interest when performing PZT sensor characterization.

Single PZT Sensor: Broadband I 3
In this section, we illustrate the differences in the calculated broadband instrumental response I 3 of a single PZT sensor with and without cropping algorithm (see Eq. (16)). In the lower part of Fig. 8a, we show the original I d 3 in a marker symbol series from 1 kHz to 1 MHz (units: dB, the superscript d denotes displacement). We use 1 V/nm as the reference sensitivity of the PZT sensors to measured displacement. We see that, without cropping some components from the final I 3 , there are remarkable amplitude variations of I d 3 in the frequency bandwidth of 1 to 7 kHz and 14 to 170 kHz, respectively.
Following the proposed algorithm described in Eq. (16), we obtain the cropped broadband I d 3 . In the lower part of Fig. 8b, we see that the 12 data series (for 0.3 to 3 mm ball diameters) overlap with negligible vertical shift. This suggests a promising stability of the developed algorithm, as well as the characterization methodology. The linked I d 3 (gray thick line) is almost flat from 10 kHz to 1 MHz with a slope of nearly 0 dB/decade and shows a strong dependence on the measured displacement. We suggest that the PZT sensors used in this study can be used as a displacement-sensitive transducer from 10 kHz to 1 MHz, an important observation to unravel the source physics associated with the AEs [9,48].
In the left lower part of Fig. 8b, from 1 to 10 kHz, the linked I d 3 has a slope of 40 dB/decade. By performing the transformation of displacement into acceleration in the frequency domain, we rewrite the displacementbased I d 3 into the acceleration-based I a 3 as where the superscript a denotes acceleration. We add the original and cropped I a 3 into Fig. 8a and b, respectively. In Fig. 8a, we find remarkable amplitude variations for the data series of I a 3 from 1 to 7 kHz. In Fig. 8b, we find that the I a 3 segments overlap well with little vertical shift and the linked I a 3 (blue thick line in the upper part of Fig. 8b) is relatively flat from 1.2 to 6 kHz. The PZT sensor is sensitive to timevarying acceleration over this frequency range, where the sensor shows similar responses as an accelerometer with a sensitivity of around 65 mV/g (V: voltage, g: gravitational acceleration).
Villiger et al. [49] installed accelerometers (Wilcoxon 736T) with a sensitivity of 100 mV/g into boreholes along the tunnel wall for in situ AE monitoring below 25 kHz. To also record higher-frequencies AEs at the same observation point, additional PZT sensors were installed adjacent to accelerometers. In a similar manner to Kwiatek et al. [11], accelerometers were used to provide a relative correction for the measurements of PZT sensors when larger events were recorded on both sensors. Unfortunately, the amount of AEs recorded by both sensor groups were relatively low and this may affect the magnitude correction factor for smaller events. With the methods presented here, having a wider and continuous understanding of the instrument response, would allow for absolute understanding of the AEs and help check the consistency in the relative calibration methods used therein. A longer overlap in frequencies between calibrated PZT sensors and accelerometers will improve the interpretation of source parameters of in situ AEs through the cross-validation of these two measurements.
In summary, we suggest that our proposed broadband calibration methodology can help to classify PZT sensors with respect to their displacement, velocity or acceleration response over a broad range of frequencies. We note, however, that this motion-dependence might not be clearly observable in PZT resonant sensors which have specific breaks in their spectral response and higher sensitivity on a series of resonance lobes.

PZT Sensor Array: Broadband I 3 Group
We now extend the analysis of the instrumental response of a single PZT sensor to that of an array of PZT sensors of the same KRN Services model. The detailed arrangement of the PZT sensor array is given in Sect. 2.2. The single sensor performance over different take-off angles is discussed in Appendix E. In Fig. 9a, we show the voltage signals near the first Pwave arrival (blue square) for all sensors, with the take-off angle ranging from 0 • to 62.9 • caused by a 0.5 mm diameter ball impact. At time t = 0, the first P-wave arrives in the ray path θ = 0 • . We see that the time of the first P-wave arrival has a positive dependence (gray line) on the take-off angle while the first peak amplitude of the P wave decreases significantly as the take-off angle increases. Also see the similar positive dependence in the first S-wave arrival. Since there are several sensors at the same take-off angle (i.e., θ = 42.6 • -A3, B2, C3, C4), we find that the corresponding shape of the voltages matches well with each other. But, due to the I 3 vari- Fig. 9 a Raw voltage signals (colored) measured by a PZT sensor array at a group of take-off angles caused by a 0.5 mm diameter ball impact. b The maximum valid frequency versus take-off angle with ball diameters ranging from 0.3 to 3 mm Through these observations, we suggest that the proposed experimental configuration realistically captures the elastic wave propagation phenomena using an array of PZT sensors. We obtain a group of numerical Green's functions from all source-receiver pairs and derive the corresponding displacements. The same procedures used for characterization of a single PZT sensor are followed. By using the algorithms described in Eq. (16), we obtain the valid frequency bandwidth over a broad range of ball diameters for all PZT sensors. In Fig. 9b, we show the ω max χ of the valid frequency bandwidth versus the take-off angle θ under ball impacts with different diameters. We find that when θ ranges from 0 • to 54.0 • , the ω max χ of all PZT sensors from the same ball impact has minor variations. Moreover, ω max χ has a positive dependence on the ball size. However, when θ increases to 62.9 • , there is weak resultant displacement in the 3-direction measured by PZT sensor C2. The derived ω max χ group corresponding to small ball diameters (< 2 mm) are out of order, suggesting that in the characterization experiments, we should keep the epicenter of PZT sensors close to the location of the ball impact. Note that we used conical shaped PZTs throughout these experiments; the corresponding analysis could be quite different for disk shaped PZTs with a large diameter [39,50].
All the PZT sensors used in this study have been used in laboratory fracturing experiments for some time (more than 3 years). Therefore they are assumed to have a similar response to mechanical vibration but are not exactly identical due to damage of the crystal material as well as uncertain variations in the manufacturing. In Fig. 10, we show wellstacked I d 3 of all the PZT sensors from 1 kHz to 1 MHz . We see that all these PZT sensors have similar displacement-and acceleration-dependence from 10 kHz to 1 MHz and from 1 to 6 kHz, respectively. Remind that for individual calibration of single sensor, the calculated I 3 integrates the effect of the interaction (e.g. static pressure, elastic wave reflection and transmission) between this PZT sensor and measured  Fig. 10 Summary of displacement-based instrumental responses (colored) over 13 PZT sensors. A1 to D1 represent the label given to the PZT sensors. We have shown displacement-and acceleration-dependence (gray dash dot line) of I d 3 over 10 kHz to 1 MHz and 1 to 6 kHz, respectively object. Regarding the collective calibration of a sensor array, I 3 of single sensor at low frequencies might be affected by interactions between other sensors and measured object. This is not evaluated in the current work but will be considered in future.

Discussion
In this study, we presented a comprehensive FEM analysis of laboratory body wave modelling to obtain accurate numerical Green's functions between an active source and an array of PZT sensors where the modelling parameters are similar to the experimental configuration. To avoid expensive computational costs on separate simulations to calculate the EGFs of a group of ball impacts, we used a Heaviside step function to describe the loading applied at the same location for the ball impacts. The resulting theoretical displacement was readily calculated by performing the convolution between the numerical Green's functions and the force-time functions of the ball impacts.
To further improve computational efficiency, we performed two separate simulations; low-frequency modelling based on a relatively coarse grid (c-FEM), and highfrequency modelling based on a relatively fine grid (f-FEM). Both models were validated against the reference approach over the high-frequency band such that high-precision FEM solutions are obtained. We performed the low-frequency analysis of wave propagation phenomena which integrated physical-based boundary conditions among the utilized experimental components. We suggest that the results from the c-FEM model have better physical foundations. Both the high-frequency validation and low-frequency extension work successfully demonstrate the capabilities of the constructed model to solve the laboratory elastic wave propagation problem.
We performed impacts tests using balls with diameters changing from 0.3 to 3 mm and obtained a group of segmented instrumental responses I s,l 3,χ with various levels of bandwidth overlap. To obtain the broadband I 3 , we developed algorithms to accurately pick the bandwidth of I s,l 3,χ by considering the accuracy of FEM solutions and quality of the experimental data. We showed that, to rigorously understand the valid frequency bandwidth of I 3 , a broad range of ball diameters are needed. Finally we obtained an accurate I 3 for an array of PZT sensors at different take-off angles with their unique source-receiver characteristics.
We realized that the accuracy and frequency range of instrumental responses greatly rely on the manner in which the theoretical displacement and the measured sourcedependent voltage signal are obtained. In the standard calibration of PZT sensor [51], theoretical displacement is derived from the generalized ray theory of which deficiency has been described in the introduction. ASTM International [51] suggests that the accuracy of the theoretical displacement applies only before any wave reflected from the boundary of the calibration block returns to PZT sensors (around 100 µs if a big steel block of 0.9 m diameter and 0.4 m thickness is adopted). This limits the frequency bandwidth significantly and in accordance to Eq. (5), the instrumental response is only valid with high resolution until around 40 kHz. However, our study suggests following our modelling methodology, PZT sensors can be calibrated down to 1 kHz while using a much smaller steel transfer plate (0.35 m × 0.35 m × 0.05 m in our study).
ASTM International [51] also adopted the fracture of the glass capillary as the source to generate elastic waves until several MHz, which is suitable to calibrate PZT sensors at high frequency range. However, due to a very short rise time (between 200 and 300 nanoseconds), this source cannot generate signals with high SNR at low frequency range, e.g. below tens of kHz, such that broadband calibration presented here cannot be proceeded with high-fidelity. In our study, ball impact with a relatively large diameter (e.g. 3 mm) can pro-duce signals with high SNR around 1 kHz (See Fig. 13b). Moreover, capillary fracture also suffers from repeatability of the test data while ball impacts using the electromagnetic delivery device described here, can produce a very repeatable source as shown in Figs. 12 and 13.
When we place conical PZT sensors from the bottom surface to the top surface of the steel transfer plate where ball impact occurs, current methodology can be extended to calibrate these sensors but using "surface (or Rayleigh) wave calibration" instead of "body wave calibration". Although the Rayleigh wave calibration of PZT sensors is outside the scope of this study [39], we suppose the calculated instrumental response keeps almost the same as which we showed in Fig. 10 since same KRNBB sensors are utilized and the point contact assumption still stands.
AE monitoring at in situ scale using PZT sensors with a relatively large aperture diameter (e.g. few centimetres) instead of the conical design are utilized of which calibration also attracts the interests of community [15,49]. If the wavelength is not significantly larger than sensor aperture, the amplitude of elastic waves will be averaged at the contact area while the true ground motion at the small tip is masked. In this case, aperture effect decreases the calculated instrumental response where the amount can be estimated by a Bessel function of wave speeds, frequency and sensor aperture [39,52]. This effect is weakest at the takeoff angle of 0 • (PZT sensor directly below the impact) since all waves arrive the sensor aperture almost in the same phases; and should be strongest at the large takeoff angle (PZT sensor on the bottom surface but has an incremental epicentral distance) until 90 • (strongest, PZT sensor on the top surface) where waves arrive out of phases. The calibrated results without properly correcting aperture effects could induce the interpretation of source parameters deviating from the truth of AEs: radiated energy, source size, corner frequency, etc. To eliminate the aperture effect, we need to advance our current methodology in both aspects of numerical modelling and deployment of larger scale laboratory experiments in the future work.

Conclusion
This study focused on developing methods that can bridge the gap between qualitative analysis and quantitative characterizations of laboratory and in situ AEs. Following the proposed methodology, we can absolutely calibrate PZT sensors and thus properly interpret the messages of ground motion from AE monitoring. In future, to study the physics of dynamic failures in complicated conditions, PZT sensors need to be characterized in prior (i.e. sensors sit inside a pressurised, fluid filled triaxial cell at high-temperature). Our constructed models can potentially be integrated with state-of-the-art multiphysics couplings from COMSOL Mul-tiphysics such that a well-validated FEM model capable at solving (multiphysics) problems of elastic wave propagation can be developed. We can then provide accurate insights into the source properties of microcrack behavior and quantitatively study the damage evolution and fracture propagation in brittle materials over a broadband frequency range and source dimensions.
Funding Open Access funding provided by ETH Zurich.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Hertzian Impact Source
We used the single impact of spherical ball as the mechanical active source to excite body waves in our transfer plate while the multibounce scenario of ball drop is outside the scope of this study [53]. In Fig. 11a, when the steel ball instantaneously impacts the steel transfer plate, the forcetime function f 3 (t) derived from the Hertzian impact theory [34] can be written as where f max = 1.917ρ 0 denotes the maximum force of the ball impact, and t c = 4.53(4ρ 1 π(δ 1 + δ 2 )/3) 2/5 R 1 v −1/5 0 is the contact time between the ball and test specimen. R 1 and ρ 1 are the ball radius and density, respectively. δ i is a material constant depending on Young's modulus, E i , and Poisson's ratio, μ i , of the ball and the plate, 2gh d represents the impact velocity due to the free drop motion, where h d and g are the dropping height and gravitational acceleration, respectively.
The amplitude spectrum of f 3 (t), denoted by F 3 (ω), is expressed as where ΔP ≈ 0.5564(t c f max ) is the change in momentum that the ball imparts to the steel transfer plate. Γ is the Gamma function [54]. Due to the properties of Γ (z) at non-positive integers, there exists a group of local minima and maxima of F 3 (ω), where Here i, j are non-negative integers. ω i 0 and ω j p are the frequency group corresponding to the local minima F 3 (ω i 0 ) of the i th and local maxima F 3 (ω j p ) of the j th lobe of F 3 (ω). Figure 11b illustrates this. When F 3 (ω) deviates from a flat plateau, spectral energy attenuates rapidly with changes of 2 to 3 orders of magnitude, from 100 kHz to 1 MHz for the force-time function of the 1 mm diameter ball impact (grey lines). The marker symbols positioned by ω i 0 and ω j p show the separation of a series of lobes in the spectral energy falloff phase. These characteristics of a series of lobes are useful to validate the theoretical f 3 (t) against the calibration at experimental data discussed in Appendix C.1.

B Short Historical Note on KRNBB Sensors
We used KRNBB sensors as the receivers to illustrate the calibration methodology throughout this study. The sensors employed here have been uniquely designed to minimize the resonant features of the PZT sensor. By adding this literature review, we hope to illustrate the conception, construction and developments of this unique variant of PZT sensors. The idea of using a conical PZT geometry to reduce resonant features was first proposed by Proctor [55]. In this design, Proctor chosed PZT-5A (a lead zirconate titanate composition) which was confirmed to be an ideal choice in comparison to other variants of piezoelectric materials (e.g. lead metaniobate, X-cut Quartz, 36 • Y-cut LiNbO, and PVDF) using experimentally and analytically evaluation [56]. In an overview of PZT sensor calibration methods at the time [57], it was noted that standard commercial (resonant)style sensors were sensitive to a combination of mechanical inputs (displacement, velocity, or acceleration). In contrast, the "optimal" transducer they highlight (i.e. that described by Proctor [55]) was, by design, only sensitive to displacements. However, Proctor's design was bulky and not easily deployed. Innovation was continued by Glaser et al. [38], which detailed all the components used to miniaturize Proctor's sensor. The design was eventually commercialized [38] and used in this study.

C Amplitude Spectral Analysis in Ball Impact Experiment
To obtain a nearly flat spectrum of the force-time function f 3 (t) from 1 kHz to 1 MHz, we dropped a series of steel balls with different diameters from the same height. The impact from balls with a small diameter (i.e., 0.3 mm) has a rise time (∼ 0.5 µs) and imposes a high-frequency mechanical pulse resulting in spectral energy concentrated well below its corner frequency (e.g. ≥ 1 MHz). Conversely, the ball impact with relatively large diameter balls (e.g. 3 mm) has a much longer rise time (∼ 4.7 µs) where most of the spectrum energy is concentrated under the low-frequency bandwidth (from 1 to 100 kHz). Larger diameter (e.g. 10 mm) ball impacts were performed and we found that the response of the PZT sensor used was saturated. For this reason, smaller diameter balls were used. Researchers [30] dropped balls of varying diameter (1.5 to 10 mm) onto a disk plate with a diameter of 103.8 mm. They found that when the contact time is extended by performing larger ball impacts, the elastic waves resulting from the impact and subsequent boundary reflections interacted with each other. As a result, the convolution between f 3 (t) and g 33 (x, t) does not hold true and u 3 (x, t) cannot be obtained. Therefore, we needed to evaluate the contact time t c and ensure that the travel time of the elastic waves was (at minimum) twice the plate thickness. Sensor saturation and convolution ineffectiveness are two factors that are under reported in the literature and, considered together, constrained the largest ball diameter to 3 mm in this study.
The same drop height h d = 138 mm was used for all impact tests; this was high enough to generate mechanical vibrations that could be measured by PZT sensors over all take-off angles. We repeated dropping the steel ball 5 times with 12 diameters (0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 1, 1.5, 2, 2.5, 3 mm). We introduced data processing techniques described in Sect. 2.2 to determine u 3 (x, t) and ψ(x, t), which were used to obtain the segmented amplitude spectra of displacement |U 3 | and voltage |Ψ |, from Eqs. 6 and 7 , respectively. Note that all the voltage spectra |Ψ | used in this section comes from PZT sensor A1 which located directly beneath the impact point.

C.1 Short-Duration Excitation (100 kHz to 1 MHz)
In this section, we windowed the raw voltage caused by shortduration (80 µs, denoted by superscript s) excitation of 0.3 and 3 mm ball impacts as shown in the inset region of Fig. 3. The time t = 0 is when first P wave arrives.
The amplitude spectra of |Ψ | s and |U 3 | s due to the 0.3 mm ball impact from 100 kHz to 2.5 MHz are shown in Fig. 12a. Since we are interested in the plateau part of the amplitude spectrum, we need to determine the corresponding corner frequency ω c . To obtain ω c , we use the Omega-n model [58] to perform fitting of voltage spectra |Ψ | s : In this study, Ω refers to voltage spectra |Ψ | s . The ω c is determined as 1.33 MHz (green line) in Fig. 12a and the fitted result is plotted as the thick gray line. The standard deviation (black error bar) of |Ψ | s among repeated tests is shown. Minor variations suggest that the ball impact is a repeatable mechanical source and works well at high-frequencies below ω c . We see the SNR (signal-to-noise ratio) is continuously larger than 1 from almost 250 kHz to 1.62 MHz and this acts as an indicator to crop a valid frequency band of |Ψ | s in Sect. 4.1. Figure 12b presents the results from a 3 mm ball impact. We note that both |U 3 | s (blue line) and |Ψ | s (red line) fall off rapidly with an almost constant slope (gray line) fitted from a series of "lobes" while the fitted ω c is fixed at 100 kHz. This means that we could not crop valid |U 3 | s and |Ψ | s segmentations from 100 kHz to 1 MHz. We suggest that these lobes are caused by a group of local minima and maxima impact that forces itself in the phase of spectrum energy fall off in accordance with Eq. (20). By labelling the local minima (i.e., 1, 2, 3, 4, 5, 6, ...) of |U 3 | s (pink dash line) and |Ψ | s (green solid line), we find there exists a horizontal offset of frequency. Reed's elastic impact theory does not work accurately at high-frequencies for relatively large diameter ball drops.

C.2 Long-Duration Excitation (1 kHz to 100 kHz)
In this section, we extend the short-duration analysis to longduration case (8000 µs, denoted by superscript l). The raw and windowed ψ(x, t) l caused by 0.3 mm and 3 mm ball impacts are shown in Fig. 3.
Following the same data processing technique given in Appendix C.1, we obtain the amplitude spectra of voltage . We see that the amplitude spectrum of the noise floor almost overlaps with or even higher than |Ψ | l below 110 kHz or beyond 1.23 MHz. Therefore |Ψ | l valid from 110 to 174.6 kHz can be used with |U 3 | l for further analysis of I 3 , see Eq. (8). By incrementally increasing the ball diameter until 3 mm, we obtain the corresponding |Ψ | l and |U 3 | l as shown in Fig. 13b. We find ω min as the minimum between ω c = 65.7 kHz of |Ψ | l and 174.6 kHz and determine the valid frequency band from 1 to 65.7 kHz, where the the mean value of SNR is around 20. Both |U 3 | l and |Ψ | l are relatively flat below ω c ; they have a distinct and similar spectrum shape regarding the local minima and maxima, which indicates the resonance and anti-resonance of the steel transfer plate due to the ball impact. By comparing Fig. 13a and b, we find that, with the longer contact times (i.e., from 0.5 to 4.7 µs) associated with larger (i.e., from 0.3 to 3 mm) ball drops, the interaction between elastic waves and the Hertz impact force history could enhance the resonance and anti-resonance observed in the experimental data. The similarity between the experimental data (voltage) and theoretical estimation (displacement) strengthened the reliability of proposed methodology for deriving low-frequency numerical Green's functions.

D Spectral Cropping and Linking
In Fig. 14, we present the general principle (left) and flowchart (right) of cropping and linking operations used in Eq. (16). Regarding a ball impact test with certain diameter, top five lines show the transformation from raw |I 3,χ | into well-constrained /| I l 3,χ | / ∪ /| I s 3,χ | / . By taking account of a broad range of diameters, we can get the instrumental response from 1 kHz to 1 MHz.

E Effect of Take-Off Angle on Broadband I 3 of Single PZT Sensor
In this section, we aim to evaluate the effect of take-off angle on broadband instrumental response I 3 while maintaining the other factors (e.g. contact, sensor, cable, A/D converter). We perform ball impact tests where the caused ground motion is recorded by the same sensor at 3 take-off angles, θ = 0 • , 42.6 • , 52.4 • , respectively. Ball diameter of 0.4, 1 and 3 mm are selected to cover broadband frequency range (from 1 kHz to 1 MHz) in accordance with the conclusion in Sect. 4.1. Same cables and A/D converter are utilized through each test. Other experimental parameters refer to Table 1. Figure 15 shows the contact manner between the measured object (upper) and PZT sensor. By measuring the length of linear springs, we maintain the pressure applied on the PZT sensor tip through tests at different take-off angles. The application of spring-loaded PZT sensors also refers to [59].
We follow the same data processing technique given in Sect. C.1 to analyze voltage spectra |Ψ | and displacement spectra |U 3 | and the proposed algorithm in Sect. 4.1 to calculate displacement-based instrumental response I d 3 . Figure 16 shows the I d 3 and associated standard deviations over three take-off angles from 1 kHz to 1 MHz. We notice these I d 3 match well at most frequencies and the standard deviation of I d 3 is mostly below 4 dB. Slight differences in I d 3 could be caused by: -additional physics not modelled but existing in the realistic elastic wave propagation due to ball impacts. For example, seismic attenuation (e.g. anelastic attenuation, scattering) causing dissipation of energy as elastic waves propagate through the steel plate. Advanced modelling methodology could be expected in future work. -experimental repeatability. Although we attempt to maintain the potential experimental setting, unexpected variations still occurred through tests.
We take an example to illustrate the effect of slight differences in I d 3 on the source characteristics of impact events. Assume this sensor detects a seismic magnitude of -5.849 (reference) due to 3 mm steel ball impact using the I d and -5.833, respectively. Comparing these seismic magnitudes, we conclude that there exist negligible differences in the analyzed source characteristics resulted from different takeoff angles while the same sensor is used through tests.