Pulsed-laser source characterization in laboratory seismic experiments

Abstract

The present study aimed to characterize the properties of a laser-generated seismic source for laboratory-scale geophysical experiments. This consisted of generating seismic waves in aluminum blocks and a carbonate core via pulsed-laser impacts and measuring the wave-field displacement via laser vibrometry. The experimental data were quantitatively compared to both theoretical predictions and 2D/3D numerical simulations using a finite element method. Two well-known and distinct physical mechanisms of seismic wave generation via pulsed-laser were identified and characterized accordingly: a thermoelastic regime for which the incident laser power was relatively weak, and an ablation regime at higher incident powers. The radiation patterns of the pulsed-laser seismic source in both regimes were experimentally measured and compared with that of a typical ultrasonic transducer. This study showed that this point-like, contact-free, reproducible, simple-to-use laser-generated seismic source was an attractive alternative to piezoelectric sources for laboratory seismic experiments, especially those concerning small scale, sub-meter measurements.

1 Introduction

Reproduction of seismic wave propagation at the laboratory scale (De Cacqueray et al. 2011; Barriere et al. 2012; Bordes et al. 2015; Valensi et al. 2015; Holzhauer et al. 2017; Pageot et al. 2017; Devi et al. 2018) is a promising approach that could lead to significant progress in imaging complex media and monitoring at the near-surface and crustal scales. The results obtained at the laboratory scale may be applied at the field scale using upscaling methods (Backus 1962; Capdeville et al. 2010; Dvorkin and Wollner 2017). In the present study, we investigated the applicability of a pulsed-laser seismic source (Martin et al. 1994; Rasolofosaon et al. 1994; Lebedev et al. 2011; Mikesell et al. 2012) in experiments as an alternative to the mechanical or piezoelectric sources frequently used at the laboratory scale.

A pulsed laser, or ablation laser, has been considered an interesting ultrasonic wave generator since the 1960s. The physical mechanisms involved in elastic wave generation via a pulsed-laser impact of high incident power density ($P_d$) were first described by White (1963). An initial thermoelastic regime exists at low $P_d$, followed by an ablation plasma regime when $P_d$ increases; however, in most cases, the two wave-generation regimes coexist when waves are created by a pulsed-laser impact (Hutchins et al. 1981; Dewhurst et al. 1982; Aussel et al. 1988).

The thermoelastic regime dominates when the effective $P_d$ is below a certain threshold; for example, Aussel et al. (1988) estimated that the theoretical threshold is approximately 15 MW/$\hbox {cm}^2$ in aluminum materials. In the low-power-density regime, transient surface heating under the laser impact leads to local thermoelastic deformations (Scruby et al. 1980), which give rise to elastic wave propagation inside the medium; in this case, a thermoelastic source can be modeled with a shear-stress dipole on the surface (Rose 1984; Dewhurst et al. 1982; Hurley and Spicer 1999; Bernstein and Spicer 2000).

The ablation regime occurs when the $P_d$ threshold is exceeded (Dewhurst et al. 1982; Aussel et al. 1988; Audoin et al. 1996); the laser-generated seismic source can then be approximated by a point source if the laser beam is well focused (Hutchins and Nadeau 1983). In this regime, the ejection of plasma matter from the material surface due to the laser impact generates a shock wave as well as an excitation force (Hosoya et al. 2014). The focused laser source is thus contact-free and impulsive in the ablation regime. If additional limiting conditions are created by surface modifications such as coatings (Davies et al. 1993; Dewhurst et al. 1982), the laser source may enter a constrained surface regime, which induces normal stress on the impacted surface. These additional boundary conditions enhance elastic compressional P-wave amplitudes for both the thermoelastic (Dewhurst et al. 1982) and ablation regimes.

Aussel et al. (1988) studied laser point sources thoroughly, both theoretically and experimentally, highlighting differences between the thermoelastic and ablation regimes. Audoin et al. (1996) applied the ablation regime to anisotropic materials to generate elastic waves and derived the stiffness coefficients from the first-arrival times; they highlighted the potential application of this contact-free method in characterizing composite materials under unusual conditions, such as high temperatures (Scruby and Moss 1993).

There have been a few attempts to use laser-generated sources in laboratory experiments following the pioneering study of Scales and Malcolm (2003), who used a laser source on natural rock samples. Campman et al. (2004, 2005) applied a linear thermoelastic laser source to undertake laboratory-scale seismic investigations in an ideal medium. Malcolm et al. (2004) showed the potential of using a ballistic laser pulse to retrieve Green’s functions from waves in the diffusion and equipartitioning regimes. This motivated Mikesell et al. (2012) to analyze the coda generated by a laser impact on an aluminum block, including tiny scattering points in the thermoelastic regime, with the purpose of retrieving scatter locations using correlation techniques. Blum et al. (2013) described a contact-free method to precisely estimate the anisotropy parameters of shale cores, while Mikesell et al. (2017) applied the same method to study Antarctic ice-core layers. Van Wijk et al. (2019) confirmed the potential of a thermoelastic laser source to probe anisotropy in pine timbers. However, these various experiments are not sufficient to validate the applicability of a pulsed laser as a seismic source for either solid-media imaging or geophysics where a large number of shots (stacking) is necessary (Chai et al. 2011; Matonti et al. 2017) to produce high-resolution, intensive data sets.

In order to gain some insight into the types of seismic waves expected under the thermoelastic and ablation regimes, we present in Fig. 1 the analytical solutions of wave propagation following a laser-generated seismic source in the specific case of aluminum. We suppose that the medium is spatially infinite, the P- and S-wave velocities are $V_p = 6350 \,\hbox {m/s}$ and $V_s = 3130 \,\hbox {m/s}$, respectively, and the receiver is located 50-mm away from the source. For both regimes displayed in Fig. 1, we observe an impulsive first arrival (P-wave) near 7.9 $\upmu$s followed by a low-frequency transient displacement (near-field), followed again by another impulsive arrival (S-wave) signature near 16.3 $\upmu$s which is dominant in the thermoelastic regime. Figure 1a) shows the seismic displacement $u_z$ (normal component) at the receiver position under the thermoelastic regime according to the formulation given by Scruby et al. (1980) and Aussel et al. (1988). This analytical formulation used a Heaviside step function to model the time dependency of the thermoelastic force dipole. Note that a negative $u_z$ after the theoretical S-wave arrival time results from the step model because the thermoelastic effect is considered endless over time. Figure 1b) shows the impulsive response yielded by a normal point source (Aki and Richards 2002), as the seismic source wavelet in the ablation regime can be modeled by a triangular function.

Fig. 1

Analytical seismic displacement $u_z$ as a function of time computed at 50 mm from the seismic source in an infinite aluminum medium. The component $u_z$ of the displacement is the normal component on the epicenter point. In case (a), the laser-generated seismic source is in the thermoelastic regime, and the displacement is computed by combining the wave equation and thermoelastic equations (Scruby et al. 1980; Aussel et al. 1988). In case (b), the laser-generated seismic source is in the ablation regime, and the source is modeled as a point source (Aki and Richards 2002)

We anticipate observing experimental measurements in aluminum similar to those computed in Fig. 1a and b. We propose in the present paper to transfer well-known results (Scruby and Drain 1990) from the field of laser ultrasonics to exploration geophysics at the laboratory scale. Specifically, we propose to experimentally and numerically characterize a pulsed-laser ultrasonic source from the viewpoint of a geophysical seismic source. The role of the input parameters governing the properties of the wave-field generated by a pulsed-laser impact is investigated throughout the paper through a series of experiments and simulations performed mainly on standard 2017A aluminum alloy, as well as one set of experiments on a carbonate core. A thorough understanding of the laser source wavelet properties under different regimes controlled by various experimental parameters is indeed essential for future quantitative studies on seismic data.

In the present study, we focused first on understanding the generation and propagation of seismic waves induced by a pulsed-laser impact in a simple case – a purely elastic, homogeneous, and isotropic medium. For this purpose, we chose to perform experiments on aluminum samples while ensuring purely elastic propagation, which can be modeled through elastodynamic approaches. Then, we demonstrated through a set of measurements on the carbonate core that we could also use the pulsed-laser seismic source in natural rocks (Lebedev et al. 2011).

We will provide a detailed description of the experimental setup in Sect. 2 of this work, alongside theoretical and numerical interpretations of some typical seismic records following a laser-generated impact on the aluminum blocks and carbonate core. A thorough characterization of the pulsed-laser source is provided in Sect. 3 and was achieved by varying the experimental input parameters. Section 4 is devoted to the measurement and modeling of the directivity pattern of the seismic displacement following a laser-generated impact, which will be compared to that of a classical piezoelectric transducer. Finally, conclusions and perspectives are given in Sect. 5.

2 Experimental setup and interpretation of seismic displacement measurements

2.1 Experimental setup

We developed an experimental setup composed of a multiple-incident-power pulsed laser, a convergent lens, a laser Doppler vibrometer (LDV) (Audoin et al. 1996; Nishizawa et al. 1998; Bodet et al. 2005; Bretaudeau et al. 2011; De Cacqueray et al. 2011; Lebedev et al. 2011), and parallelepipedal aluminum blocks of various thicknesses, as seen in Fig. 2. In the experiments, multiple and controllable-incident-power laser beams were emitted from the pulsed-laser apparatus in the normal direction with respect to the surface of the aluminum samples; the impact of the pulsed-laser beam on the aluminum acted as a seismic source and therefore emitted elastic waves throughout the volume of the block. The convergent lens was used to modify the laser impact diameter and consequently the incident power density $P_d$ (Dewhurst et al. 1982) on the impacted surface. The normal (with respect to the surface) components of the seismic displacement field $u_z$ were then precisely measured by the LDV on the opposite surfaces of the aluminum blocks along the direction of the unit vector $\varvec{e_z}$, as shown in Fig. 2.

Fig. 2

Experimental setup dedicated to wave propagation investigations of parallelepipedal aluminum blocks. $\textcircled {1}$: pulsed Q-switched laser; $\textcircled {2}$: convergent lens; $\textcircled {3}$: target, a 50 mm thick parallelepipedal aluminum block; $\textcircled {4}$: laser Doppler vibrometer (LDV). The pulsed-laser beam emitted by the Q-switched laser is shown in yellow; the laser beam of the LDV measuring the seismic displacements is shown in red. The red spot right in the middle of the target is the epicentral point where most of the measurements were performed, directly behind the pulsed-laser impact on the aluminum target. The unit vector $\varvec{e_z}$ is normal to the aluminum surface

The laser was a Q-switched pulsed Nd:YAG (neodymium-doped yttrium aluminum garnet) laser radiating at energies reaching 850 mJ (Q-smar.t 850 from Quantel, $\textcircled {1}$ in Fig. 2). The laser emitted pulses with a diameter of 9 mm over a few ns (approximately 7-20 ns depending on the power requested) at a wavelength of 1064 nm and maximum pulse repetition frequency of 10 Hz. The total incident laser power could be modulated by the internal $\hbox {FLQS}_{{delay}}$ (time delay in nanoseconds between the flash lamp and Q-switch pulses) parameter of the instrument. In order to quantify the incident power density $P_d$ of the Q-smart 850 as a function of $\hbox {FLQS}_{{delay}}$, we used experimental estimates derived from the following expression:

$$\begin{aligned} P_d = \frac{P \times T}{FWHM \times S}\, \, , \end{aligned}$$

(1)

where P is the power of the laser beam in W, T is the shooting period between successive laser pulses during measurement in s, FWHM is the Full Width at Half Maximum of the laser impulse in s and S is the surface of the laser beam in $\hbox {m}^2$. We used a thermal sensor (S370C from Thorlabs) specifically designed to measure the emitted power P of the short and energetic laser pulses from the Nd:YAG laser to calibrate the incident laser power as a function of $\hbox {FLQS}_{{delay}}$. We also used a silicon photodetector (EOT-2030 from Electro-Optics Technology) to measure the FWHM of the laser pulses as a function of $\hbox {FLQS}_{{delay}}$. We concluded from the measurements that the power density of the laser typically varied within $P_d = 0.3$-$470 \hbox {MW}/\hbox {cm}^2$ as a result of adjusting the $\hbox {FLQS}_{{delay}}$ internal parameter of the Q-smart 850.

The original laser beam emitted by the Q-switched laser in the experiments could then be partially or totally focused onto the target surface by a convergent lens with a focal length of 14 cm (see $\textcircled {2}$ in Fig. 2). Generally, the Q-switched laser emitted at low repetition rates (between 1 and 2 Hz), and a few tens of consecutive laser-generated seismic raw traces were stacked to improve the signal-to-noise ratio. However, a higher shooting frequency of 6 Hz (Expression (1)) was applied during the laser-power calibration in order to guarantee a continuous heat flow. The measurement frequency of the thermal sensor was also set to 6 Hz to ensure good synchronization between laser shooting and laser power assessment.

The pulsed-laser beam was directed towards the “center”  of one of the two largest faces of the parallelepipedal aluminum blocks (see $\textcircled {3}$ in Fig. 2). We investigated three aluminum parallelepipeds with dimensions of $280 \times 280 \,\hbox {mm}^2$ and different thicknesses, namely 10, 50, and 100 mm, which corresponded to approximately 3, 15, and 30 times the P-wavelength under a dominant frequency of 2 MHz in aluminum, respectively. The various sizes of the aluminum blocks were chosen such that the direct seismic waves measured on the opposite side of the samples were not disturbed by lateral reflections.

The seismic displacements generated by the pulsed-laser were then measured by the LDV (see $\textcircled {4}$ in Fig. 2) on the opposite side of the aluminum blocks. More precisely, the vibrometer measured normal displacements $u_z$ on the surfaces of the aluminum blocks. The LDV was composed of a single-point laser head (OFV-505 from Polytec), a vibrometer controller (OFV-5000), and an auxiliary decoder (DD-300) that exported the seismic displacements. Stepper motors guided the the spatial movements of the LDV laser head, allowing millimeter-scale displacement measurements along three orthogonal axes with a spatial precision of less than a tenth of a mm. The LDV signals were recorded by a digital oscilloscope (DSO-S 054A from Keysight) at a sampling frequency of 100 MHz.

2.2 Seismic displacement measurements at the epicenter of the aluminum blocks

In order to characterize of the waveform of the pulsed-laser seismic source, we started with a strong laser beam ($P_d \simeq 147\, \hbox {MW}/\hbox {cm}^2$) and hit the aluminum blocks directly in the center of the plate. We measured $u_z$ via the LDV at the epicentral distance (see the middle red spot in Fig. 2). Figure 3 shows three epicentral seismic records measured for the aluminum samples of different thicknesses after stacking 20 successive laser-generated signals.

Fig. 3

Measurements of the seismic displacements $u_z$ vs. time, obtained by LDV at epicentral positions. Experiments were performed with $P_d \simeq 147 \,\hbox {MW}/\hbox {cm}^2$ on three different parallelepipedal aluminum blocks of a 10 mm, b 50 mm, and c 100 mm thickness. L1 stands for the first longitudinal arrival, L3 for the first two-way reflection, L5 for the second reflection, ans so on. Theoretically, the first arrival in (a) is expected to be at $1.57 \,\upmu \hbox {s}$, while the measured first break arrived at approximately $1.70 \,\upmu \hbox {s}$. This slight delay of $0.14 \,\upmu \hbox {s}$ was caused by slight desynchronization due to triggers of the electronic devices and was corrected by a simple time shift. The vertical dashed lines in (a), (b) and (c) are the theoretical arrival times for the L3 and L5 waves

Fig. 4

a Laser-generated seismic source $S_0(t)$ extracted from the recorded signal in the 50-mm-thick aluminum block shown in Fig. 3b before and after the filtering. The first arrival in Figure 3b has been first isolated in time and then deprived of its low-frequency components after Butterworth filtering. b) Frequency spectra $\widehat{S_0}(f)$ of the seismic source shown in (a). c 1$^{st}$-order time derivatives ${S_0'}(t)$ of the signals in (a); the values are thus the velocities of the surface vibration. d Frequency spectra $\widehat{S_0'}(f)$ of the seismic source time derivatives shown in (c). The spectra in b and d indicate a similar frequency peak at approximately 2 MHz for the source expressed in displacement and velocities, whereas the lower frequency band fades away in the spectrum of the source time derivative. Specially, (b) highlights the effect of the Butterworth filter on $S_0(t)$

We clearly identified the first longitudinal arrivals, denoted as L1 in Fig. 3, whose amplitude decreased as the thickness of the aluminum block increased. The multiples L3, L5, and L7 in Fig. 3 corresponded to waves bouncing back and forth inside the horizontal section of the aluminum block. The thinner the block, the more reflections occurred during the 60-$\upmu \hbox {s}$ acquisition time of the three experiments. These results are consistent with the works of Dewhurst et al. (1982) and Audoin et al. (1996). The signatures of the near-field waves (Aki and Richards 2002), as predicted in Fig. 1, are also clearly visible in Fig. 3b and 3c; these two cases were particularly remarkable owning to the slope breaks occurring at $\hbox {t} \simeq 17 \,\upmu \hbox {s}$ in Fig. 3b and $\hbox {t} \simeq 33 \,\upmu \hbox {s}$ in Fig. 3c. These events will be discussed later, in Sect. 2.3.

We concentrated on the case of the 50-mm-thick block in Fig. 3b in order to analyze the waveform generated by the pulsed-laser impact. We isolated the first arrival L1 and applied a $1^{st}$-order Butterworth filter cutting at 0.2 MHz. The residual low-frequency component was then subtracted from the whole signal to retain only the high-frequency component, yielding the pulsed-laser source wavelet; the time signal $S_0(t)$ and frequency spectrum $\widehat{S_0}(f)$ of that source (before and after the filtering), as well as the corresponding time derivative $S_0'(t)$ and its frequency spectrum $\widehat{S_0'}(f)$, are shown in Fig. 4. We concluded that the laser-generated seismic source in aluminum had a dominant frequency bandwidth between 0.8 and 2.5 MHz, as shown in Fig. 4b and d. Given that the longitudinal P-wave velocity was 6350 m/s in the aluminum alloy, the wavelengths of the P-waves were approximately 3 mm in our experiments at 2 MHz; this explained why it was more difficult to distinguish the different bouncing waves within the 10-mm-thick block in case (a) of Fig. 3 compared to cases (b) and (c). Strictly speaking, since this “wavelet”, $S_0(t)$ in Fig. 4a, was not zero-mean, it was not truly a mathematical wavelet (Mallat 1999); nevertheless, we chose to use this experimental pseudo-wavelet in the following section, where numerical simulations will be compared to experimental measurements. Figure 4c displays the derivative $S_0'(t)$ of the wavelet, which is, in contrast, approximately zero-mean. The main difference between $\widehat{S_0}(f)$ in Fig. 4b and $\widehat{S_0'}(f)$ in Fig. 4d is found in the lower-frequency band. Before Butterworth filtering, $\widehat{S_0}(f)$ has a dominant lobe centered at 0 Hz, while $\widehat{S_0'}(f)$ features one main peak at approximately 2 MHz; after the filtering, the lower-frequency components in the initial $\widehat{S_0}(f)$ are considerably attenuated, preserving the similar main peak at approximately 2 MHz.

2.3 Theoretical interpre.tation of the seismic epicentral record

In order to interpret the experimental signals shown in Fig. 3, we applied an analytical approach of wave propagation in an elastic, homogeneous, isotropic, and unbounded medium. A point dislocation source was used to represent the pulsed-laser source. The analytical formulation of the i-component of the displacement ${\varvec{u}}$ at any point inside the medium is given by Aki and Richards (2002):

$$\begin{aligned} \begin{aligned} u_i({\varvec{x}},t)&= S_0 (t) *G_{ij}\\&=\frac{1}{4\pi \rho }(3\gamma _i\gamma _j-\delta _{ij})\frac{1}{r^3}\int _{r/V_p}^{r/V_s}\tau S_0(t-\tau )d\tau \\&\quad +\frac{1}{4\pi \rho V_p^2}\gamma _i\gamma _j\frac{1}{r}S_0\left( t-\frac{r}{V_p}\right) \\&\quad -\frac{1}{4\pi \rho V_s^2}(\gamma _i\gamma _j-\delta _{ij})\frac{1}{r}S_0\left( t-\frac{r}{V_s}\right) \end{aligned} \end{aligned}$$

(2)

where ${\varvec{x}}$ is the vector position, t is time, $S_0(t)$ is the point source time function in direction j at the origin, $G_{ij}$ are the Green’s functions, r is the distance between the source and the receiver, $\gamma _{i,j}$ denotes the direction cosines of the receiver positions, $\delta _{ij}$ is the Kronecker delta function, $\rho$ is the volume mass density of the medium, and $V_p$ and $V_s$ are the velocities of the longitudinal and transverse seismic waves, respectively. We applied Formulation (2) according to our experimental setup, meaning that the source displacement $S_0(t)$ was applied along the same direction i as the measurement $u_i$. Note that the boundary effects occurring in the parallelepipedal samples were not taken into account in this formulation.

The integral term in Formulation (2) describes the near field, which can be considered a convolution between the source function and a straight-line function defined between $t=r/V_p$ and $t=r/V_s$. The second and third terms in Formulation (2) give the P-wave and S-wave displacements, respectively, whose waveforms are identical to that of the unidirectional point source.

The blue curves in Fig. 5 are the computed epicentral traces from Formulation (2), and the waveform of the point source $S_0(t)$ is shown in Fig. 4. The medium was the standard aluminum alloy used previously in the computation shown in Fig. 1, and the amplitudes of these signals were rescaled such that the amplitudes of the simulated first arrivals in Fig. 5 match the experimental ones.

Fig. 5

Comparison between the experimental epicentral $u_z$ (black curves) and an analytical solution of the wave generated by a point dislocation source in an elastic, homogeneous, isotropic, and unbounded medium. The comparison is presented for the three distinct experiments: a 10-mm-thick, b 50-mm-thick, and c 100-mm-thick blocks. The seismic source $S_0 (t)$ in Formulation (2) is the retrieved laser-generated seismic source $S_0(t)$ shown in Fig. 4a. The simulated amplitudes are rescaled by a common coefficient to approximately match the experimental first-arrival amplitudes in (a), (b), and (c)

We observed good correlations between the modeled and experimentally derived first arrivals in terms of frequencies and waveforms, as seen in Fig. 5. The less convincing modeling in Fig. 5 for case (a), that of the 10-mm-thick aluminum, was the result of comparing a theoretically unbounded model with experimental displacements that were strongly affected by multiple reflections. However, the recorded near field was very well reproduced by the analytical solution of the three aluminum plates of different thicknesses. It was also remarkable that the so-called near field is still visible after propagating 100 mm in Fig. 5c, which corresponds to approximately 30 times the P-wavelength (or approximately 60 times the S-wavelength). Finally, we note that the slope breaks in the experimental records shown in Fig. 5b and c and mentioned in Sect. 2.2 were qualitatively reproduced through the modeling and thus clearly originated from the near-field wave interactions.

2.4 Numerical interpretation of the seismic epicentral record

In order to further develop the quantitative interpretation of the experimental data, we conducted numerical simulations of wave propagation in the aluminum blocks. We used the Hou10ni numerical code developed by the Inria team project Makutu (ancient Magique3D), which can simulate wave propagation in elastodynamic media using the interior penalty discontinuous Galerkin method (IPDGM). We refer to Barucq et al. (2014) for a description of the features of Hou10ni for elastodynamic media in the frequency domain and to Riviere (2008) and Grote et al. (2006) for more details on the IPDGM. This approach introduced discontinuous basis functions instead of the continuous ones used in classical finite element methods (Ciarlet 2002; Raviart et al. 1998). One of the advantages of this method is the block-diagonal mass matrices yielded by the discontinuous functions, which can be more easily inverted. Another advantage, hp-adaptivity (Houston and Süli 2001; Houston et al. 2002; Baldassari et al. 2011), allows for adding additional finite elements or locally increasing the order of the discretization to improve accuracy in some critical areas if necessary, while coarser meshes can be implemented in simpler zones. This method is particularly efficient for highly heterogeneous natural rocks such as carbonates.

The equations (Baldassari 2009) that we solved numerically with Hou10ni to simulate elastodynamic wave propagation in the aluminum blocks are shown below:

$$\begin{aligned} .&\displaystyle {\rho \frac{\partial ^2 {\varvec{u}}}{\partial t^2}} - \nabla \cdot (\sigma ) = {\varvec{0}} \quad \text {in} \quad \Omega \times I \end{aligned}$$

(3)

$$\begin{aligned}&\sigma \cdot {\varvec{n}} = \varvec{S_0(t)}\delta ({\varvec{x_0}}) \quad \text {on} \quad \partial \Omega \end{aligned}$$

(4)

$$\begin{aligned}&{\varvec{u}}({\varvec{x}},0) = {\varvec{0}} \quad \text {in}\quad \Omega \end{aligned}$$

(5)

$$\begin{aligned}&\frac{\partial {\varvec{U}}}{\partial t}({\varvec{x}},0) = {\varvec{0}} \quad \text {in} \quad \Omega \end{aligned}$$

(6)

$$\begin{aligned}&\text {with}\quad \sigma = C\varepsilon ({\varvec{u}}) = \frac{1}{2}C(\nabla {\varvec{u}}+\nabla ^T{\varvec{u}}) = C\nabla {\varvec{U}} \end{aligned}$$

(7)

where $\rho$ is the volumetric mass density, $\varvec{{u}}$ is the displacement vector, ${\sigma }$ is the stress tensor, ${\Omega }$ is the space domain of the medium, I is the time domain, ${\varvec{n}}$ is the vector normal to the boundary surface $\partial \Omega$, $\varvec{S_0(t)}$ is the experimentally measured source vector (see Fig. 4) perpendicular to the boundary surface and applied at position $\varvec{x_0}$, $\delta ()$ is the Dirac function, C is the high-order stiffness tensor (which links the stress tensor to the strain tensor), and $\varepsilon ({\varvec{u}})$ is the strain tensor. Equation (3) is actually the microscopic form of Newton’s second law and Equation (4) is a Neumann boundary condition, where the source term $S_0(t)$ simulating the impact of the pulsed laser is null at the boundary except at $\varvec{x_0}$.

Fig. 6

Schematic view of $3^{\text {rd}}$-order tetrahedron meshes used by Hou10ni for 3D simulations of the aluminum blocks. There were approximately 2.4 million tetrahedrons in total

2D and 3D numerical simulations of wave propagation were performed with triangle (for 2D) and tetrahedron (for 3D) elements and 3$^{\text {rd}}$-degree polynomials. Figure 6 is an example of a bulk portion of the aluminum blocks meshed with 3$^{\text {rd}}$-order tetrahedrons.

We carefully considered the properties of the mesh used for the numerical calculations. Due to the practical issues of parallel computing center memory limits, the 3D computations were not performed within the real volumes of the aluminum used in the experiments but with cubes with side lengths equal to the thickness of the blocks. The two largest dimensions of the aluminum blocks were then reduced to the scale of the thickness in the calculations in order to obtain reasonable total number of finite elements. Consequently, absorbing boundary conditions (ABCs) were applied to the four lateral sides of the cubes during computations in order to avoid sidewall reflections, which were of little interest in this study.

Thus, we calculated the wave propagation in the $280 \times 280 \times 50 \,\hbox {mm}^3$ experimental aluminum block using the reduced dimensions $50 \times 50 \times 50 \,\hbox {mm}^3$ to numerically simulate the elastodynamic wave propagation in the cube with ABCs at four of its six surfaces. The average side length of the tetrahedrons was of the order of 1 mm in those simulations. We specifically tested the influence of the size of the tetrahedrons in the 3D calculation and concluded that a tetrahedron side length of approximately 0.5 mm was necessary for accurate simulation results; this is similar to 1/6$^{th}$ of the P-wave wavelength. This was in accordance with the results obtained by Baldassari et al. (2012), where it was reported that 10 points per wavelength allowed for a relative accuracy of 1% in 2D acoustic simulations. The admitted mesh featured 2.4 million tetrahedrons and met the side-length criterion to ensure accuracy.

Fig. 7

Comparison between the experimental records and the 2D/3D simulations for the 50-mm parallelepiped. L stands for longitudinal wave and T for transverse. a The epicentral displacements $u_z$ in different configurations. b The derivatives $u_z'$ of the records in (a). All of the data (experiment and simulations) were recorded at the epicentral distance

Figure 7a displays the experimentally measured displacements at the epicentral point for the 50-mm-thick aluminum block in black, as well as the corresponding simulated traces obtained from Hou10ni for the 2D and 3D contests in red and blue, respectively. The computations were performed with the experimentally measured seismic source displayed in Fig. 4a, written as $\varvec{S_{0} (t)}$ in Equation (4). The amplitude of the raw 2D simulation data followed a typical 2D geometrical expansion, which were not directly comparable with other data in a 3D context. Operto et al. (2006) applied a time-dependent correction to 3D data in order to perform quantitative comparisons with 2D data; hence, we rescaled the amplitudes of the 2D simulation data by a factor of $\frac{1}{\sqrt{t}}$ so that they became comparable with the 3D data. The raw 3D simulation simulation signal contained significant numerical noise due to the finite element size, error accumulation upon iteration, and imperfect ABCs. The frequency content of the main part of these numerical noises intersected the signal frequency band, making it difficult to filter them via low-pass filters. A noise gate technique based on the energy profile of the signal (Inouye et al. 2014) was thus applied in addition to the low-pass filters. This technique works well when the time signal is approximately zero-mean within an arbitrarily chosen signal window. Figure 4 shows that the velocity (derivative) profile should meet this condition. Therefore, Fig. 7b shows the 1$^{st}$-order time derivatives of Figure 7a, namely the velocities of the epicentral vibration for the aluminum block. It is obvious that the traces in Fig. 7b are zero-mean within relevant windows. Before calculating the derivatives, all of the traces expressed in displacement (Fig. 7a) were low-pass filtered with an infinite-duration impulse response (IIR) filter with attenuation starting at 1.5 MHz and stopping at 3.5 MHz in order to eliminate high-frequency noise above 3.5 MHz. Meanwhile, the amplitudes of the 2D and 3D numerical simulation data were rescaled such that the maximum vibration displacement of the first P-arrival matched the experimental records. In particular, the time-derivative signal of the 3D numerical simulation in Fig. 7b underwent a noise-gating process. The final displacement data of the 3D simulation displayed in Fig. 7a were ultimately derived from the integral of the processed velocity signal.

Both the 2D and 3D numerical simulations of the displacement gave satisfactory results compared with the experimental epicentral data in Fig. 7a, particularly in terms of primary longitudinal arrivals, followed by the near field, which sprawls in time. The main experimental features were clearly obtained in the 2D and 3D simulations. However, we note that, although the 2D simulations were relatively noise-free numerically, they contained a continuous positive component of amplitude; the various numerical tests we conducted to understand this anomaly indicated that this trend originated numerically from the seismic source, which did not fulfill the zero-mean condition of a mathematical wavelet. Despite these amplitude anomalies in the 2D calculations, the simulations led to remarkably clear waveforms in accordance with the experiment signals recorded by the LDV.

Furthermore, the time derivatives of the simulated traces in Fig. 7b fit the experimental data even better, owing to removal of the pervasive near-field signatures and continuous positive component of the 2D simulation. Meanwhile, the relative amplitudes of the principal arrivals (L1, L3, L5, and L7 in Fig. 7b) in both the 2D and 3D simulations were in good agreement with the experimental data. Obviously, the realistic feature of the 3D simulation in terms of wave amplitude was one of its major advantages over the 2D simulation, the latter requiring time-dependent correction in order to be comparable with the experimental data.

At first glance, the 3D simulation data in Figure 7a and b seem noisier than the 2D computations, as the waveforms corresponding to experimental data are difficult to recognize. However, the issue of the continuous positive component of amplitude in 2D was negligible in the 3D simulations. Our systematic 3D tests revealed that the ABC method used in Hou10ni was not completely satisfactory, since the current numerical scheme still yielded a considerable number of side-reflected waves (L2 for example), particularly in the case of S-wave reflection, such as the strong peak (T2) in Fig. 7a and b. These imperfect ABCs, which were the origin of most of the undesired high-frequency arrivals (L2, T1, L3, and T2) in the 3D simulation, also explain why the simulation data (both displacement and velocity) became increasingly oscillatory as time progressed.

Fig. 8

3D representation of an elastodynamic calculation performed with Hou10ni for a cube with a side length of 50 mm. Wave-fronts on three orthogonal sections crossing the center of the cubic aluminum model are shown at time t = 6 $\upmu$s after the initial pulsed-laser impact. The color-coded amplitudes are absolute values of the displacement vectors ${\varvec{u}}$. Reflections on the four lateral faces (parallel to the z direction) are largely canceled due to the ABCs. The other two sides have Neumann conditions. This type of 3D representation inside the medium clearly offers a comprehensive view and better understanding of wave propagation compared to the single-point (epicenter) data shown in Fig. 7a for the same calculation. We can highlight for example in this calculation P- and S-bulk wave-fronts, head waves (Aussel et al. 1988; Ilan and Weight 1990; Gridin 1998), and surface waves, as well as the spatial distribution of the near field

As an illustration, Figure 8 shows a snapshot of wave propagation for three representative sections inside the 50-mm-side-length aluminum model at approximately $6 \,\upmu \hbox {s}$ after the initial seismic impulse. We can clearly see the spatial amplitude distribution of wave fronts at this moment and identify the different waves highlighted in Fig. 8. Comparisons between the experimental and numerical data helped us to validate our numerical methods and models, as well as thoroughly understand the origin and kinematic behavior of the different arrivals during the experiments.

2.5 Seismic displacement measurements on a carbonate core

In order to explore the applicability of the seismic pulsed-laser source to heterogeneous, natural rocks, we tested the source on a carbonate core extracted from the Urgonian geological formation (Shen 2020). The core had a diameter of 10.7 cm and a height of 30 cm. When obtaining the set of measurements presented in Fig. 9, the pulsed laser impacted the cylindrical core at mid-height, and the LDV measured the displacement of a horizontal section at the periphery of the core over a total angle of 300$^\circ$, which was symmetrical with respect to the impact of the pulsed laser. We protected the surface of the carbonate core with an aluminum coating, which also acted as the backing material of the pulsed-laser source.

Figure 9a shows a seismogram of one laser shot, with the laser source and 109 LDV receivers all surrounding a single section of the core, while Fig. 9b highlights 9 of the 109 seismic traces. Owning to the particularity of the cylindrical core geometry, the surface waves arrived earlier than in the parallelepipedal aluminum blocks and dominated the seismogram in terms of energy, which clearly exhibited a linear dependence on angle, as shown in Fig. 9a. The first P-wave arrivals were clear, despite having weaker amplitudes than the surface waves. Other arrivals, such as the first S-wave arrivals and the reflected body waves, were also identifiable but more difficult to distinguish.

The first P-wave arrival times are marked with green crosses in Figure 9a and b. Interestingly, the heterogeneity of the carbonate core was revealed, as the sample exhibited deviation from the regular hyperbolic shape of the first P-arrivals; moreover, the heterogeneity of the core was seen in the arrivals following the P-arrivals, particularly the global pattern of the seismogram in Fig. 9a which is quite asymmetric with respect to the central receiver. The spectral content of the pulsed-laser wavelet following the first break peaked at a maximum frequency of approximately 1 MHz, which was less than half of the maximum observed in the aluminum blocks. This can be explained by the fact that natural rocks attenuate higher frequencies more effectively than elastic media. Finally, we demonstrated that the first seismic arrival times, which carry essential information for seismic tomography, anisotropy estimation, and waveform inversion, were perfectly recorded using a seismic pulsed-laser source on a natural carbonate rock sample, as seen in Fig. 9.

Fig. 9

Application of the seismic pulsed-laser source on a natural, cylindrical carbonate sample: a the seismogram of one shot, i.e., the arrival times of the successive seismic waves vs. positions of receivers around the core for a horizontal circular section; b highlighting 9 of the 109 seismic traces. We clearly distinguish first P-wave arrivals as well as strong and delayed surface waves

3 Pulsed-laser seismic source characterization

3.1 Repeatability and stability of the pulsed-laser source

A key advantage of using a Q-switched laser seismic source in geophysical laboratory experiments, where the amplitudes of seismic waves are critical, is the stability of the incident laser power through the repetition of impulses over time. In order to study this stability, we first quantitatively calibrated the level of the power emitted by the pulsed laser. First, the thermal sensor (see Sect. 2) was directly hit with the pulsed-laser beam at a repetition rate of 2 Hz for a total of 1000 impacts for one $\hbox {FLQS}_{{delay}}$, which itself varied from $0\, \upmu \hbox {s}$ to $255 \,\upmu \hbox {s}$ by step sizes of $10 \,\upmu \hbox {s}$. The thermal sensor measured the integrated power on its 25-mm-diameter surface for each laser impulse; the integrated power was then divided by the laser spot size multiplied by the FWHM of the pulse following Expression (1) in order to convert it to a power density $P_d$ in units of $\hbox {MW}/\hbox {cm}^2$.

Fig. 10

Distribution of incident power densities $P_d$ deduced from measurements. Four different levels of energy were studied: $\hbox {FLQS}_{{delay}} = 200$, 150, 100, and $50 \,\upmu \hbox {s}$ in (a), (b), (c), and (d), respectively. Each distribution was obtained with a set of 1000 successive pulsed-laser impacts repeated at a frequency of 2 Hz, with the average value $\upmu$ of $P_d$ being highlighted by a red line and the standard deviation $\delta$ around the mean value by two dotted blue lines

Figure 10 presents the probability density distribution, arithmetic mean, and standard deviation of the measured $P_d$ for 4 of the 26 measurement series of the different investigated energy levels. The standard deviations of the four histograms in Fig. 10 were very small (less than 1% in cases (b), (c), and (d) compared to their arithmetic mean values. These four chosen series served as representative examples of all the studied energy levels of the pulsed laser; these measurements demonstrated that, for all pulsed-laser energies, the incident power density of the laser beam was generally stable during a series of successive laser pulses.

Fig. 11

Distribution of first-arrival seismic wave amplitudes $u_z$ deduced from LDV measurements performed at the epicentral distance for the 10-mm-thick aluminum block. Four different levels of incident power density were studied: $P_d$ $\simeq 0.3$, 28, 90, and $323 \,\hbox {MW}/\hbox {cm}^2$ in (a), (b), (c), and (d), respectively. Each distribution was obtained with a set of 300 successive pulsed-laser impacts. The average value $\upmu$ of the measured $u_z$ is indicated by a red line, while the standard deviation $\delta$ around the mean value is highlighted by two dotted blue lines

For a second time, we quantified the impact of 300 successive laser pulses on the 10-mm-thick aluminum block sample for four different incident power densities, $P_d \simeq 0.3$, 28, 90, and $323\, \hbox {MW}/\hbox {cm}^2$, with a focused 3-mm-diameter laser beam on the impacted surface of the aluminum block. The seismic displacements $u_z$ were recorded at the epicentral position (see Fig. 2), and the analyses of the amplitudes were performed for the first P-wave arrival of each shot. Figure 11 shows the distributions of $u_z$ measured by the LDV for 300 shots at each incident power density, as well as the average values and standard deviations.

For the three highest $P_d$ in Fig. 11, the measured amplitudes showed relatively small standard deviations with respect to the average values; however, the standard deviations were much larger (of the order of 5%) for cases (b), (c), and (d) compared to the distribution of $P_d$ shown in Fig 10. For these three cases, the relatively dispersed measurements clearly demonstrated that the impacted surface state evolved through the successive pulsed-laser impacts, and that this temporal evolution contributed to the amplitude variations. While each dataset in Fig. 11 was obtained by fixing the laser-shot position, the aluminum surface could be modified under the ablation regime during an acquisition requiring thousands of shots, which contributed to the volatility of the pulsed-laser output and can partially explain the relatively wide amplitude distribution in Fig. 11. The case of a relatively small $P_d$ shown in Fig. 11a demonstrated that the measured seismic displacements were obviously not well distributed. Furthermore, the sensitivity limit of our LDV measurements had been reached, compounded by being in a low energy state, which was unfavorable for the stability of the Q-switched laser.

Fig. 12

Evolution of the amplitudes $u_z$ of the first-arrival seismic wave measured via the LDV at the epicentral position for the 10-mm-thick aluminum block as a function of 800 successive laser shots. The incident power density $P_d \simeq 0.3 \,\hbox {MW}/\hbox {cm}^2$, and the diameter of the focused laser beam was 3 mm. The blue curve was obtained with a fixed position of the laser shot on the aluminum block, whereas the red curve was obtained by changing the impact position at shot numbers 199, 350, 450, 550, 650, and 750, as indicated by the vertical gray lines

We performed another set of experiments in order to test the hypothesis that surface deterioration due to the successive laser pulses was responsible for the variations in amplitudes measured by the LDV. The amplitudes in Fig. 12 were obtained for two different configurations at $P_d \simeq 0.3\, \hbox {MW}/\hbox {cm}^2$. The blue curves show the amplitude evolution measured at the epicentral position with continuous shots on the 10-mm-thick aluminum block, during which the position of the laser impact did not change with time. In the other case, namely the red curves, we changed the laser impact position periodically; thus, each vertical gray line in Fig. 12 marks a fresh source spot for the red curves. It is clear that, each time the impact position was changed, the original stable amplitude appeared again in the LDV measurements. The measurements shown in Fig. 12 indicate that it was possible to derive an empirical law predicting the decrease in $u_z$ as a function of shot numbers (fitting of the blue curve) when the targeted position was fixed. This empirical law could then be used to correct the values of $u_z$ (red curve) in order to recover the “real stable amplitude level” , as highlighted by the horizontal line in Fig. 12.

We concluded from this stability study that the Q-switched laser generated stable and reproducible pulsed-laser beams within an optimal range of emitted power. We showed that the heat of successive laser pulses could alter the surface state through the thermoelastic regime, and this alteration may induce amplitude variations of the associated seismic waves measured by the LDV setup. We also verified experimentally that a moderate ablation regime caused surface alterations mainly via material consumption, yet this induced relatively smaller variations in the measured displacements compared to the surface heating.

3.2 Induced seismic displacement vs. the pulsed-laser incident power density

In this section, we establish the quantitative relationships between the incident power of the pulsed-laser beam $P_d$ estimated experimentally using Expression (1), and the induced seismic displacements $u_z$ measured by the LDV. We used the same experimental procedure as that described in Sect. 3.1 with a focused laser beam with diameter of 3 mm, where the pulsed laser was successively shot at the three aluminum blocks of varying thickness 20 times while measuring $u_z$ at the epicentral distance; this operation was repeated for $P_d$ values ranging from approximately 0.3 to $470\, \hbox {MW}/\hbox {cm}^2$ with step sizes ranging from 20 to $50\, \hbox {MW}/\hbox {cm}^2$. Note that we changed the impact location on the aluminum blocks (after averaging every 20 shots) at the beginning of each new series in order to minimize the uncertainty in amplitude measurement mentioned in the previous Sect. 3.1.

Fig. 13

a Normalized amplitudes of the displacement $u_z$ measured at the epicentral distance vs. the incident power density $P_d$ for the three aluminum blocks of different thicknesses. The amplitudes have been normalized relative to that at 10 mm from the source in order to eliminate the geometrical expansion. b Fourier spectra performed on the first arrivals of the measured seismic displacements $u_z$ vs. $P_d$ for the 50-mm-thick aluminum block

We compiled all of the measurements in Fig. 13a. The amplitudes of the displacements measured by the LDV were normalized by the ray lengths (thickness of the aluminum block), thus eliminating the effect of geometry spreading on the amplitudes, which should vary at a rate of 1/r for the aluminum. The three curves in Fig. 13a exhibit approximately the same trend, thus validating the normalization and bulk-wave origin of the seismic signals measured by the LDV. These curves were of $1^{st}$-order accordance with the results of Dewhurst et al. (1982), who concluded that the longitudinal first-arrival amplitude was proportional to the laser output energy above a certain $P_d$ threshold in the ablation regime. However, it was rather difficult to observe a regime shift from these measurements above and below the ablation threshold near $P_d = 15 \,\hbox {MW}/\hbox {cm}^2$ (Aussel et al. 1988).

We also performed fast Fourier transforms (FFTs) of the first arrivals of the LDV measurements in order to characterize the frequency spectra of the pulsed-laser-generated seismic waves as a function of $P_d$. The spectra generated from the 50-mm-thick aluminum block are shown in Fig. 13b as an example. It is obvious that the frequency spectra remained similar for all $P_d$ within the ablation regime, in which the useful dominant frequency of the generated wavelet was approximately 1.7 MHz (1.2 MHz to 2.2 MHz), as already observed in Fig. 4b. We obtained similar results for the other two aluminum blocks.

Considering the results in both Figs. 11 and 13, we decided to adopt an $\hbox {FLQS}_{{delay}} = 150~\upmu \hbox {s}$ in the following experiments because it corresponded to an intermediate laser power level, which yielded relatively small standard deviations, as shown in Figs. 11b and 13a. Even in a well-focused case, as shown in Figure 11b, the ablation was not intense ($P_d = 28 \,\hbox {MW}/\hbox {cm}^2$).

3.3 Influence of the lens-target distance on the pulsed-laser regime

The distance d between the convergent lens and the target in our experimental setup, as shown in Fig. 2, controlled the laser spot size on the target surface and consequently the incident laser power density $P_d$ received by the aluminum blocks. Three distinct situations could occur, which depended on the value of d, as seen in Fig. 14:

• For d > f, where f is the focal length of the convergent lens, the laser beam focused before reaching the target (depicted in Fig. 14 on the left side of the central target). In this case, the laser beam collapsed at the focus point, releasing a large amount of energy and creating a plasma due to ionization (Couairon and Mysyrowicz 2007) of the ambient air. This was presumably a favorable situation for the thermoelastic regime since the incident power density $P_d$ reaching the target may have been below the ablation regime threshold due to a substantial loss of energy in the air during the collapse. Note that the sound wave generated by the laser collapse propagated in the air until reaching the target; that sound wave then converted at the air-target interface into a low-frequency P-wave propagating within the volume of the target. The low-frequency sound wave could then be fortuitously recorded by the LDV in our experiments. Although this wave trailed behind the direct laser-generated P-wave arrival, it could possibly obscure the identification of the direct wave.

• For $d < f$, the laser beam did not collapse since it encountered the target on its path before focusing (Fig. 14 on the right side of the central position). In this case, the laser beam directly hit the surface of the target after partial focusing, ejecting plasma out of the surface. The propulsive plasma acted like a mechanical source, leading to high-frequency elastic wave propagation within the block. When the plasma ejection was intense enough to generate a sound, such as in the d > f case, a strong, low-frequency P-wave originating from the sound accompanied the high-frequency elastic waves and arrived at the other side of the target at the same time as the first P arrivals. The sound wave, which was not included in the numerical simulations, contributed to the differences between the laboratory measurements and the 3D simulations (Fig. 7).

• For d $\simeq$ f, the laser beam focused near the surface of the target (central target in Fig. 14). Depending on whether focusing occurs right before or right after the target, the collapse may occur or not. In the no-collapse case, the sub-millimeter laser beam diameter reached the target with a maximum $P_d$. Similar to the $d > f$ and $d < f$ cases, a low frequency sound wave coexisted with the laser-generated P-wave, both of them propagating at the same velocity.

Fig. 14

Schematic of the laser beam trajectory originating from the Q-switched laser (on the right). The parallel beam encounters a convergent lens (in blue) with focal length $f = 150\, \hbox {mm}$. The targets, depicted as gray planes, are at a distance d from the lens; they can then be impacted before the focusing of the laser beam if $d < f$, right at the point of focusing with a possible collapse of the beam if $d \simeq f$ or after the focusing of the beam if $d > f$. A substantial part of the incident laser power is lost in the air in the $d >f$ case, where a collapse occurs before the beam reaches the target

Fig. 15

a Seismic displacement $u_z$ measured by the LDV at the epicentral distance for the 50-mm-thick aluminum block, with various lens-target distances d from 30 to 270 mm. b, c, and d are the seismic displacement records seen in (a) at distances of 255, 75, and 125 mm, respectively. The seismic records are depicted in blue at negative displacements and red at positive values. The background colors indicate different regimes and the sound wave arrivals. b, c, and d highlight a typical thermoelastic record, mixed record, and typical ablation record, respectively

Gradually changing the values of d in the experimental setup covered the three distinct situations presented in Fig. 15; we modified the lens-target distance d in our experiments step by step to explore a wide range of laser-generated source regimes (thermoelastic or/and ablation). We hit the 50-mm-thick aluminum block with the laser while varying d varying from 30 to 270 mm with step sizes of 5 mm and a lens focal distance f of 150 mm. For each d, $u_z$ was recorded by the LDV at the epicentral position after stacking 30 successive shots as shown in Fig. 15a.

As shown in Sect. 1, the first-break wavelets of the seismic data in the thermoelastic regime were negative displacements (Scruby et al. 1980; Dewhurst et al. 1982) forming a transposed right-trapezoid (Fig. 1a), while those of the ablation regime featured an impulsive, positive displacement (Fig. 1b). In order to identify and highlight these regimes, the negative and positive displacements are filled in with blue and red, respectively, in Fig. 15a. Thereafter, we identified three different regimes:

• when the lens-target distance was either sufficiently small ($d<45 \,\hbox {mm}$) or large ($d>235\,\hbox {mm}$) compared to the focal length of the lens f, we mainly observed the thermoelastic regime, which is indicated by light blue in Fig. 15a. Figure 15b highlights the seismic record at d = 255 mm, which in this case appeared to be very close to the analytical example for the thermoelastic regime, as shown in Fig. 1a.

• At intermediate values of d, namely 50 mm $<d<95\,\hbox {mm}$ and $195 \,\hbox {mm}<d<230\,\hbox {mm}$, the thermoelastic and ablation regimes coexisted in Figure 15a, as indicated by gray. A typical example of this regime is shown in Fig. 15c for $d = 75 \,\hbox {mm}$. In this dual regime, the first positive arrival clearly originated from the ablation regime, while the thermoelastic regime yielded the subsequent signal.

• When d was close to f, as indicated by the pink areas of Fig. 15a, the laser beam was well focused, and the incident power density of the laser was well above the ablation threshold. As shown in Fig. 15d for $d = 125 \,\hbox {mm}$, the seismic signal was dominated by a clear, impulsive, positive first break, which was characteristic of a pulsed-laser source under the ablation regime. Furthermore, we identified a group of sound waves characterized by a relatively low frequency at approximately 40 kHz in Fig. 15a indicated in green. As mentioned earlier, the origin of these sound waves was either the ablation sound on the surface of the target when $d < f$ or the laser beam collapsing in the air due to focusing when $d > f$; in the latter case, when $d > f$, the sound wave arrived at the aluminum target at $t = S_{air} \times (d-f)$ after the laser beam collapsed ($S_{air} = \frac{1}{340}$ s/m being the slowness of sound in air), which agrees with the slope of the arrivals in the green area of Fig. 15a.

Fig. 16

a First arrival amplitude vs. incident power density for different regimes. Points marked by circles were obtained with a laser beam focusing in the air, thus incurring energy loss. Triangle points were achieved under shorter lens-target distances without energy loss in air. Colors indicate different regimes, as identified in Fig. 15. b Approximate regime mapping following the “impact radius vs power density” plan. The regime limits are supposed to be linear and pass through the origin of the plan. The two lines (curves under the logarithmic x-axis scale) separate the three regimes while approximately adhering to the regime color codes of the two scatter series in (a)

We also studied the different regimes by considering the amplitudes of the first arrivals measured by the LDV as a function of the incident power density of the laser beam, as seen in Fig. 15a. When $d < f$, as we knew the power emitted from the Q-switched laser and the focal distance f of the lens, it was straightforward to compute the size of the laser impact from geometrical optics rules and the incident power density $P_d$ on the target as a function of d. Accordingly we show the amplitudes of the first arrivals and the size of the laser impact as a function of $P_d$ for the $d<f$ case in Fig. 16a and b, respectively, denoted by triangles. The curve composed of triangles in Fig. 16a is in good agreement with that of Dewhurst et al. (1982). as mentioned previously, in the $d>f$ case, part of the incident power was lost in the air during the collapse (Couairon and Mysyrowicz 2007); this part was estimated experimentally to be as much as 40% of the initial laser power emerging from the Q-switched laser in the present study. This rough estimation enabled completing Fig. 16a and b for the $d>f$ case with circles.

At first glance, the colored zones in Fig. 16 allows identification of the different regimes already characterized in Fig. 15 (according to the same color coding). Thus, thresholds can be inferred to separate the different regimes, expressed in terms of $P_d$: the thermoelastic effect dominated for $P_d < 5 \,\hbox {MW}/\hbox {cm}^2$, the mixed thermoelastic-ablation regime occurred when $5 \,\hbox {MW}/\hbox {cm}^2< P_d < 20 \,\hbox {MW}/\hbox {cm}^2$, and the significant ablation regime occurred when $P_d = 20 \,\hbox {MW}/\hbox {cm}^2$. Figure 16a and b also indicate distinct differences in the trends of the upper (triangles) and lower (circles) branches as a function of $P_d$: these results clearly demonstrate that the first P arrival amplitude of the seismic data depended not only on the incident $P_d$ or the absolute laser energy (Dewhurst et al. 1982), but also on the effective impact area. It is worth noting that, for both $d<f$ and $d>f$, the first arrival amplitudes $u_z$ increased until a maximum level and then decreased despite the continuously rising power density. In general, the amplitude was weaker for $d>f$ compared to $d<f$ under the same incident power density, as expected due to the laser-beam collapse in the $d>f$ case. We also note that the sound-wave disturbance portion marked by the green stars in Fig. 16 is characterized by very low amplitudes for the first seismic arrival due to a very high incident power density but an extremely small impact on the target. Finally, we present an approximate map in Fig. 16b that delimits the regimes considering two parameters: $P_d$ and the laser-impact radius. The thermoelastic regime dominated when a low incident power density was used with a relatively large impact radius, whereas ablation was predominant as soon as the incident power was above $20 \,\hbox {MW}/\hbox {cm}^2$, regardless of the impact size.

Fig. 17

Fourier spectra amplitude vs. lens-target distance. The spectra were obtained by computing the FFT of the first arrivals of the seismic traces shown in Fig. 15a. The spectra for the short and long lens-target distances are similar to cardinal sine function shapes, with the strongest lobe at approximately 30 kHz. Another remarkable lobe at higher frequencies ($\sim 2 \,\hbox {MHz}$) is visible between $\hbox {d} = 100\, \hbox {mm}$ and $\hbox {d} = 150 \,\hbox {mm}$, as well as a weaker one between $\hbox {d} = 150 \,\hbox {mm}$ and $\hbox {d} = 200 \,\hbox {mm}$. The cardinal sine-shaped lobes correspond to thermoelastic wavelets, while the two lobes centered at 2 MHz were produced by the ablation wavelets

We showed in Sect. 3.2 that the frequency content of the laser-generated waves in the ablation regime depended very little on the input power density. In the thermoelastic regime, this frequency depended on pulsed-laser regimes. We computed the absolute values of the FFT after extracting the dominant portion of the wavelets from Fig. 15 for each d, as seen in Fig. 17. Spectra of selected wavelets from the thermoelastic regimes were similar to cardinal sine function shapes with their main lobes located between 10 and 100 kHz, since the characteristic thermoelastic wavelets were close to gate functions (see Figs. 1a and  15b). The spectra of ablation-origin wavelets ($100 \,\hbox {mm}< d < 190 \,\hbox {mm}$) also had cardinal sine signatures of approximately 10–100 kHz due to the presence of the near field in these measurements performed for the 50-mm-thick aluminum block; note that the spectra of the sound waves were also merged with those of the near field, particularly for $100< d < 150 \,\hbox {mm}$, where the sound waves were generated simultaneously with the ablation-origin wavelet. Furthermore, we observed a higher frequency band between 1 and 2 MHz, originating from the ablation effect, which corresponds to the impulsive and broadband (Audoin et al. 1996) wavelet, as seen in Fig. 17. Note that the cardinal sine spectra are absent for certain lens-target distance ranges in Fig. 17, where the positive near-field contribution canceled out the negative thermoelastic contribution in the time domain.

In conclusion, the pulsed-laser source was broadband in our experiments performed on aluminum, ranging from a few kHz to a few MHz, as seen in Fig. 17. This property is particularly suitable for studying seismic dispersion, especially for poroelastic (single or double porosity) or anisotropic media, where their attenuation properties may be highly frequency-dependent (Dupuy et al. 2016). It may also be well adapted to the study of seismoelectric effects (Pride 1994), whose frequency dependence is not well known (Devi et al. 2018).

4 Radiation pattern of pulsed-laser seismic sources

A real seismic campaign deploys large amounts of receivers associated with different sources along the surface, producing datasets that can be viewed as seismograms. Thus, we extended our epicentral measurements of $u_z$ to an entire receiver line with a constant interval of 1 mm to study the radiation pattern of the different regimes obtained by the pulsed-laser seismic source and to compare them with the radiation pattern of a typical piezoelectric source, as shown in Fig. 2.

The preceding analysis, as described in Sect. 3.3, demonstrated that ablation was the main origin of the recorded seismic waves when the lens-target distance d was 140 mm; in contrast, the seismic waves originated predominantly from the thermoelastic regime when d was 260 mm. We chose these two distances to run the linear measurements depicted in Fig. 2 in order to study both regimes. An extra acquisition was done with a high-frequency piezoelectric transducer (C548-SM from Panametrics, $\hbox {diameter} = 10 \,\hbox {mm}$, $\hbox {nominal frequency} = 1 \,\hbox {MHz}$, Ricker waveform) as the seismic source located where the Q-switched laser beam hit the aluminum target.

Fig. 18

Seismogram ($u_z$) registered by the LDV along a straight line (see Fig. 2) on the surface of the 50-mm-thick aluminum block. a) $d = 260 \,\hbox {mm}$, the waves originating mainly from the thermoelastic regime. b) $d = 140 \,\hbox {mm}$, the waves originating mainly from the ablation regime. c) Using a P-wave piezoelectric transducer (1 MHz) as the seismic source

Figure 18 shows the three seismograms measured for the 50-mm-thick aluminum block under different regimes or sources. At first sight, these seismograms were rather similar in terms of arrival times, and we easily identified the first longitudinal arrivals and their reverberations, which are labeled as L1, L3, and L5 in Fig. 18. Shear wave arrivals were also clearly identified in the seismograms based on their well-known travel times; they are labeled S1 and L2S1, the latter being the S-wave converted from the L2 longitudinal wave on the aluminum-air interface. Seismic reflections and conversions coming from the side of the aluminum blocks at offsets of $\pm 140 \,\hbox {mm}$ were also clearly observed as straight lines bouncing on the side boundaries of the three records in Fig. 18.

Looking now more closely at the three seismograms, where negative displacements are colored in blue and positive ones in red, we deciphered some differences among them. Thermoelastic depressions are visible, especially in the central part of Fig. 18a where wave propagation trajectories are the shortest, indicating that the negative displacements due to the thermal expansion attenuate quite rapidly in space with geometric spreading. Similarly, in Fig. 18b, we observe a central zone where the displacements are first positive, a signature of the positive near-field waves in the ablation regime, followed by the negative signatures of the sound wave mentioned in Sect. 3.3. Finally, as seen in Fig. 18c, the data recorded using the piezoelectric source had less continuous (low frequency) components than those of Fig. 18a or 18b. This can most likely be attributed to the quasi-zero-mean waveform of the piezoelectric source; in contrast, the waveform of the pulsed-laser source was not zero-mean, as mentioned previously in Sect. 2.4.

Fig. 19

Radiation patterns of the normal seismic displacement $u_z$ extracted from the onset P-wave arrivals measured by the LDV with a a pulsed-laser seismic source in the ablation regime, b a piezoelectric transducer, and c a pulsed-laser seismic source in the thermoelastic regime. The radiation patterns are expressed in nm vs. angles in degrees. The experimental data are shown as solid black lines with rms values shown as error bars, accompanied by modeled radiation patterns in red lines

First-arrival amplitudes were then extracted from the three seismograms displayed in Fig. 18 and processed to obtain radiation patterns for the various seismic sources. We first filtered all the seismic traces in the seismograms by applying a low-pass IIR filter with a pass-band upper frequency of 4 MHz, which eliminated the high-frequency noise. Then, the amplitudes (maximum positive values) of the first P-wave arrivals were measured directly for each record. We then normalized the amplitudes by the geometric spreading factor (1/r). We show the experimental radiation patterns obtained from the LDV measurements of $u_z$ for the three different cases as black curves in Fig. 19. In terms of amplitudes, it is clear in Fig. 19 that the displacement $u_z$ generated by the piezoelectric transducer was at most three times larger than that generated by the pulsed laser in the ablation regime. The displacements in the ablation regime are one order of magnitude larger than those generated in the thermoelastic regime; this agrees with our experimental results in which the thermoelastic regime was obtained by the post-collapse laser beam.

The experimental radiation patterns were then compared with different models, depicted as red curves in Fig. 19:

• The measured radiation pattern of the ablation regime in Fig. 19a was compared with that of a point source simulated using the numerical code Hou10ni, as presented in Sect. 2.4. The agreement between the point-source model and the measurements was perfect at small incident angles, while the two patterns tended to differ slightly above an aperture angle of 30$^{\circ }$. However, it should be noted that the reflections and conversions occurring at the aluminum-air boundary at large angles influenced the amplitudes of the experimental first arrivals considerably according to the Zoeppritz equations (Castagna and Backus 1993), whereas the numerical simulation allowed for measuring the inside of the medium, thus eliminating boundary influences. This difference may explain the mismatch between the experiment and the point-source model in Fig. 19a for angles beyond 30$^{\circ }$.

• The piezoelectric radiation pattern appeared to be rather anisotropic, with a strong radiation in the normal direction and weaker ones on the two sides, as seen in Fig. 19b. This pattern was expected for a collimated piezoelectric transducer. We used a polar, anisotropic, analytical expression of type $A=A_0\cdot \cos ^6(\theta )$ to fit well the experimental radiation pattern in Fig. 19b where $A_0$ is the amplitude at the epicenter position.

• The radiation pattern of the laser thermoelastic regime in Fig. 19c differs noticeably from the two previous cases, showing minimum values at small incident angles and two lobes with maximum values pointing at approximately $45^{\circ }$. A theoretical P-wave radiation pattern was derived by Bernstein and Spicer (2000), who represented the thermoelastic regime as a line source composed of a shear-stress dipole. Our data agree well with this model at small incident angles, as seen in Fig. 19c. For incident angles larger than 30$^{\circ }$, the agreement is still qualitatively satisfactory, even if the experimental data are rather noisy and the theoretical modeling did not take into account the free boundary condition of the experiment. The radiation pattern in Fig. 19c is rather similar to patterns observed from S-wave seismic sources (Fukushima et al. 2009), suggesting that the waves originated from a shear-stress motion that was directly converted into longitudinal displacement at the source-targeted interface under the thermoelastic regime.

5 Conclusions and perspectives

We have explored the possibility of using a pulsed-laser beam emitted by a Q-switched laser as a seismic source for geophysical laboratory experiments. Our experimental setup was composed of a pulsed laser, a convergent lens focusing the original laser beam (9 mm in diameter), aluminum blocks, and a laser Doppler vibrometer (LDV) to measure the seismic displacement at the surface induced by the pulsed-laser source.

We quantitatively studied different regimes generated by the impact of the pulsed laser by varying the incident power density of the beam $P_d$ on the aluminum blocks. Three types of regimes (Aussel et al. 1988) were identified: 1) when $P_d \le 5 \,\hbox {MW}/\hbox {cm}^2$, a thermoelastic regime dominated (thermoelastic deformation); 2) an intermediate regime occurred when $5 \le P_d \le 20 \,\hbox {MW}/\hbox {cm}^2$, where both thermoelastic deformation and the slightly ablated surface contributed to the wave generation; and 3) when $P_d \ge 20 \,\hbox {MW}/\hbox {cm}^2$, the ablation regime dominated, creating an excitation force oriented normal to the impacted surface (Aussel et al. 1988; Dewhurst et al. 1982). The third regime led to a remarkably broadband seismic source via ablation in the aluminum, which ranged approximately from kHz to MHz; the wide frequency band of this pulsed-laser source differed noticeably from the rather narrow bands of the other mechanical or piezoelectric seismic sources used in geophysical laboratory experiments. The 2D and 3D numerical elastodynamic modelings of the wave propagation in the aluminum blocks following the pulse-laser impacts were performed via finite element approach featuring the interior penalty discontinuous Galerkin method (Baldassari 2009). We demonstrated that the 2D modeling was sufficient to obtain good agreement between the experiments and simulations in terms of arrival times of the different wave fields (P-wave, S-wave, and near field). However, the 3D simulations were required to obtain reliable amplitudes that were comparable with those of the experiments.

In a dedicated set of experiments, we extracted the amplitudes of the first arrivals in the seismograms measured by the LDV and produced radiation patterns according to the seismic displacements linked to three types of excitation: thermoelastic or ablation regimes via pulsed-laser excitation, and piezoelectric transducer excitation. The modeling of the experimental measurements demonstrated that the thermoelastic pulsed-laser source behaved as a shear-stress dipole (Bernstein and Spicer 2000; Scruby et al. 1980). The ablation pulsed-laser source turned out to be similar to a point-source force normal to the impacted surface, and most of the energy in the case of the piezoelectric transducer source was rather collimated compared to that of ablation, which was also in the normal direction with respect to the impacted surface.

We conducted our study mainly in an ideal, homogeneous, elastic, and isotropic aluminum alloy to avoid non linear effects induced by heterogeneous and/or anisotropic materials on the source. Multiple characteristics of the pulsed-laser source used throughout the study (contact-free seismic source, easy handling and orientation of the pulsed laser, reproducible incident power density $P_d$, varying seismic source regimes vs. input $P_d$, large seismic spectrum of the pulsed laser-source vs. input parameters, and simple and regular radiation patterns of the seismic sources) demonstrated that it was truly suitable for geophysical laboratory experiments involving seismic wave propagation. As part of an applicability study of the pulsed-laser seismic source on natural rocks, we also conducted experiments on carbonate rocks (Shen 2020). The experimental measurements of the carbonates confirmed that the pulsed laser could be used as a seismic source on natural heterogeneous/anisotropic rocks at the laboratory scale.

The pulsed-laser seismic source opens new perspectives for various applications, such as precise non-destructive test on metals and micro seismic exploration of small- and intermediate-scale rock samples of random shapes in the laboratory; we are especially interested in its eventual geophysical applications in rock mechanics (Dupuy et al. 2016; Shen 2020), rocks, or digital rock imaging for geological reservoir exploration.

Finally, another important perspective resulting from the present study is the quantitative characterization of the displacements measured on an interface where all the wave reflections and conversions are considered in the analysis. That would require beforehand a numerical “amplitude vs angle/offset”  (Castagna and Backus 1993) study, where the amplitude attributes of seismic data would be quantitatively modeled and explored before establishing a simulated reference for experimental data processing. This will be of particular interest when efficiently and precisely exploring the experimental attributes of amplitude.