ABSTRACT

Over the past few decades, the potential of collocated measurements of 6C data (3C of translational and 3C of rotational motion) has been demonstrated in global seismology using high-sensitivity, observatory-based ring laser technology. Proposed applications of 6C seismology range from tomographic reconstruction of near-receiver structure to the reduction of nonuniqueness in seismic source inverse problems. Applications to exploration problems have so far been hampered by the lack of appropriate sensors, but several applications have been proposed and demonstrated with array-derived rotational motion estimates. With the recent availability of, for example, fiber-optic-based high-sensitivity rotational motion sensors, widespread applications of 6C techniques to exploration problems are in sight. Potential applications are based on, for example, the fact that the extended set of combined translational and rotational motion observations enables carrying out array-type processing with single-station recordings such as wavefield separation and surface-wave suppression. Furthermore, measuring the rotational component (curl) of the seismic wavefield enables direct isolation of the S-wave constituents and could significantly improve S-wave exploration. Rotational measurements provide estimates of the spatial wavefield gradient at the free surface that allow carrying out analyses such as local slowness estimation and wavefield reconstruction. Furthermore, rotational motion measurements can help to resolve wavefield infidelity introduced by seismic instruments that are not well-coupled to the ground.

INTRODUCTION

Seismological observations of rotational motions induced by earthquakes and recorded with ring lasers (Figure 1; e.g., McLeod et al., 1998) as well as rotational motion measurements on exploration scale estimated from receiver arrays (e.g., Muyzert et al., 2012) are testimony that rotational motions induced by seismic waves are measurable and interpretable. An example of rotational motion recorded by a ring laser following the Mw 7.9 earthquake in Nepal on 25 April 2015 is shown in Figure 2; Figure 3 displays an exploration-scale data example of rotational-motion measurements computed using geophone arrays and finite-difference techniques. Such rotational-motion measurements extend the set of seismic observables beyond conventional translational measurements and open up possibilities to improve existing methods as well as develop novel techniques for subsurface exploration at all scales.

Conventionally, seismic exploration of earth’s internal structure is carried out with sensors that output up to 3C of translational motion due to a seismic wavefield (translational motion of the sensor casing relative to an inertial mass). The sensor motion is then interpreted as an estimate of the ground motion described by a vector field of up to 3C of displacement, particle velocity, or particle acceleration. However, the rotational components of the vector field and, hence, the fact that a ground-mounted sensor is also subject to seismically induced rotational motions have so far largely been ignored. In analogy to translational motion, measures of rotational motion are the angular position (angle), angular velocity (rotation rate), and angular acceleration. Rotational-motion measurements may be interpreted as estimates of the curl of the displacement field. Analogously, hydrophone (or pressure) measurements may be interpreted as the divergence of the displacement field. Seismic experiments measuring rotational motions are still rare because sufficiently reliable and sensitive enough rotational sensors useful for applications in global seismology and seismic exploration in the required frequency bands did not exist until recently.

Measuring in addition to 3C of translational motion also 3C of rotational motion extends the set of seismic observables to 6C. Such 6C measurements potentially allow (1) an improved understanding of the local ground motion, for example, due to the inherent link between the curl of the displacement and S-waves. (2) Rotational-motion point measurements at the free surface provide local estimates of the spatial wavefield gradients and, hence, provide spatial information on the seismic wavefield that would otherwise only be measurable with large and dense receiver arrays. Spatial gradients of the seismic wavefield appear, for example, in spatial wavefield filters and wavefield reconstruction algorithms. (3) Additional seismic observables enhance statistical wavefield analyses and filtering techniques. (4) Rotational motion measurements can help to resolve vector-measurement infidelity introduced by seismic instruments that are poorly coupled to the ground. Such infidelity can be caused by poorly coupled sensors that tilt or rotate upon the arrival of a seismic wave or due to external forces such as ocean currents.

As indicated in several recent reviews (e.g., Igel et al., 2015; Schmelzbach et al., 2016, 2017; Li and van der Baan, 2017a), there is a significant potential for collocated observations of translational and rotational motion in improving seismic inverse problems for the earth’s structure and sources at all scales. In regional and global seismology, rotational-motion data have been used to, for example, improve the characterization of earthquake sources (Bernauer et al., 2014; Donner et al., 2016; Reinwald et al., 2016) and inverting for subsurface structure (Bernauer et al., 2009; Fichtner and Igel, 2009), observing and investigating ambient noise sources (Hadziioannou et al., 2012; Tanimoto et al., 2015), and strong-motion earthquakes (e.g., Stupazzini et al., 2009), as well as in long-period seismology (Igel et al., 2011; Nader et al., 2014). Furthermore, understanding rotational motion induced by earthquakes is of great importance to seismic hazard assessment (e.g., Abdel-Ghaffar and Rubin, 1984; Castellani and Boffi, 1989; Basu et al., 2015).

The emerging field of rotational seismology has greatly benefited from high-resolution measurements of the vertical component of rotational motion using ring-laser technology (e.g., Igel et al., 2005, 2007). Even though these instruments are primarily dedicated to geodesy measuring the earth’s rotation rate (Figure 1), they are also sensitive to rotational motions induced by earthquakes and ocean waves (e.g., Schreiber et al., 2006, 2009). Ring-laser observations have led to several studies investigating the potential of rotational-motion observations for solving seismic inverse problems (see the recent review by Igel et al., 2015) with a view to earthquake seismology (for a historical review of rotational seismology, see Lee et al., 2009). A ring-laser data example of rotational motion observations combined with conventional translational motion recordings following the Mw 7.9 earthquake in Nepal on 25 April 2015 is shown in Figure 2, together with the rotation-based determination of the phase velocity and back azimuth (see also the analysis presented in Figure 1). Note that the phase velocity and back azimuth were determined from the displayed single-station recordings only.

Recent progress in rotational seismology related to measuring and interpreting rotational motions was reported at the meetings of the International Working Group on Rotational Seismology (IWGoRS, 2017). A seismic event database has recently been set up with access to Python/ObsPy scripts for standard processing of translational and rotational motion recordings (Krischer et al., 2015; Salvermoser et al., 2017).

In exploration seismology, rotational observations have also recently gained increasing attention. This interest is driven by the expectation to overcome long-standing problems such as coherent-noise suppression and sampling requirements. Rotational-motion data have been used to, for example, facilitate S-wave interpretation (e.g., Sollberger et al., 2016a) and to develop single-station noise suppression schemes (e.g., Barak et al., 2014; Sollberger et al., 2018). Furthermore, the inherent link of rotational-motion measurements taken at the free surface to the spatial gradient of the seismic wavefield has been exploited to implement wavefield reconstruction schemes for spatially aliased data (Muyzert et al., 2012). In addition, the direct measurement of uncontaminated rotational motions can remove the well-known ambiguity between horizontal component recordings and tilts (rotational motions around horizontal axes), which deteriorates data quality of, for example, ocean-bottom recordings (e.g., Lindner et al., 2016). One particularly interesting aspect of collocated translational and rotational motion measurements is that the 6C observations with a single station enable the extraction of wavefield characteristics such as propagation direction, local slowness, and wave type that would otherwise require the analysis of data from spatially extensive sensor arrays.

An exploration-scale field-data example recorded on the ETH Zurich campus (Switzerland) is shown in Figure 3. The 6C data were acquired in a high-resolution seismic experiment using a sledgehammer as a vertically directed source and cross-shaped arrays of five closely spaced 3C receivers to record the particle velocity. A finite-difference approximation (e.g., Liang and Langston, 2009) was used to estimate rotation rates from the geophone array data. In the case of a horizontally layered subsurface, energy would be expected in the vx, vz, and ω˙y components only (x corresponding to the inline direction; see Figure 1c for naming conventions). The arrivals observed on vy and ω˙x are likely due to the subsurface topography and heterogeneity. Note that incident P-waves convert into downgoing S-waves at the free surface, which gives rise to rotational motion observations at P-wave arrival times.

Even though global and exploration seismology investigate the earth at different scales, there is a mutual interest in improving the acquisition and analysis of rotational-motion data. The aim of this paper is to demonstrate the potential of combined measurements of translational and rotational motion for seismic exploration as well as seismology in general. We want to foster the exchange of ideas between global and exploration-scale seismology. The paper is structured as follows: We first present some of the most pertinent theoretical foundations such as the link between strain and rotation. This is followed by a brief review of instrumental concepts to measure rotational motions. The core of the paper is the discussion of potential applications of 6C data. Finally, we conclude with an outlook on the currently most promising application fields and open issues.

THEORETICAL FOUNDATIONS

Displacement, strain, and rotation

Exploration and global seismology depend on recordings of the spatial components (x=[x,y,z]T) of displacement u(x,t), particle velocity u˙(x,t), or particle acceleration u¨(x,t), with dots denoting the time derivatives. Note that vectors appear in boldface (e.g., u), components of vectors (and tensors) are expressed with subscripts such as ui, with i=x,y,z, and Einstein’s summation convention applies on repeated indices. In the following discussion, we consider a particle that is located at x at time t0 (Figure 4). At time t, the particle has moved due to deformation to the coordinates ξ=x+u with respect to the same inertial reference frame. To describe the deformation, we will consider the relative motion of the particle at x to a nearby point at x+δx with the vector line element δx defining a line segment between the two points. Keeping the first two terms of a Taylor series expansion of u around x to express displacements as varying linearly over distance |δx| yields (e.g., Aki and Richards, 2002)  
ui(x+δx)=ui(x)+juiδxj+O(|δx|2),
(1)
where i denotes the spatial derivatives in the x-, y-, or z-direction. The first term ui(x) represents a rigid-body translation because all points in the neighborhood of x share the same displacement. The second term corresponds to the displacement gradient tensor, which is a second-rank tensor that can be further separated into symmetric and antisymmetric parts yielding when ignoring higher order terms:  
ui(x+δx)=ui(x)+12(jui+iuj)δxj+12(juiiuj)δxj.
(2)
Assuming that the displacement gradients are small (|jui|1), the symmetric and antisymmetric parts of the displacement gradient tensor can be associated with the infinitesimal strain tensor,  
εij=12(jui+iuj),
(3)
and the infinitesimal rigid-body rotation tensor,  
θij=12(juiiuj).
(4)
The infinitesimal rigid-body rotation has only three independent components that can be represented by a vector Θ with elements Θk=(1/2)ekijθij, where eijk is the Levi-Civita symbol (permutation tensor).
The strain tensor ε can be decomposed into volumetric (principal) strain and shear (deviatoric) strain:  
ε=(xux000yuy000zuz)+12(0yux+xuyzux+xuzxuy+yux0zuy+yuzxuz+zuxyuz+zuy0).
(5)
The trace of ε corresponds to a change in volume relative to the initial volume (dilatation), and it is equivalent to the divergence of the displacement field: Δ=εii=·u. Hydrophones and pressure sensors measure Δ. The off-diagonal elements of εij, ij, correspond to shear strains and represent a change in shape without a change in volume (the second term of equation 5). For infinitesimal deformations, shear strain is equal to the change in an initially right angle, which is associated with a rotation.
To determine the total rotation associated with deformation that will be measured with a rotational sensor, we consider the vector line element δx (see Figure 4). Deformation will cause a change in the length as well as the orientation of δx resulting in a deformed state δξ. Defining a unit vector δx^=δx/|δx|, the total rotation ω of δx after deformation is the cross-product of δx with δξ assuming infinitesimal strain (see Appendix A for the derivation; e.g., Segall, 2010):  
ω=δx×δξ|δx|2=δx^×(ε·δx^)+Θ(Θ·δx^)δx^,ωi=eijkεkmδx^jδx^m+ΘiΘjδx^jδx^i.
(6)
The total rotation ω is the sum of a strain component due to shear strain (the first term) and the rigid-body rotation (the second and third terms). The strain component is the cross-product of the unit vector δx^ and the relative displacement dui=εijδx^j. The strain component describes the relative rotation of the line element around its center of mass with the rotation sense following the right-hand rule (see Figure 1c), whereas rigid-body rotation is relative to some other point. For seismological studies that are concerned with seismic-wave propagation and, hence, strain in a continuum, the strain component of ω is of primary interest, whereas rigid-body rotation (Θ; e.g., earth’s rotation) can be ignored.

Elastic-wave equation and curl of the wavefield

Whereas the foregoing discussion was concerned with kinematic aspects, in the following we will present some of the relevant dynamic aspects. For isotropic media, the components of the stress tensor σij in a source-free region are related to the strain tensor εkl (equation 3) by the isotropic elastic constitutive equation:  
σij=(λδijδkl+μ(δikδjl+δilδjk))εkl,
(7)
where λ and μ are the Lamé constants and δij is the Kronecker delta. The elastic-wave equation following from the equation of motion as well as equations 3 and 7 can be written as  
ρu¨=f+(λ+2μ)(·u)μ×(×u),
(8)
where ρ is the density and f denotes a body force. Equation 8 allows us to describe the propagation of a displacement perturbation through a homogeneous isotropic elastic medium.
Lamé’s theorem states that for a displacement field u, there exists potentials ϕ and Ψ such that (Aki and Richards, 2002)  
u=ϕ+×Ψ,
(9)
with ϕ and ×Ψ being called P- and S-wave components of u, respectively. Because ϕ is curl-free, we find that  
×u=××Ψ.
(10)
Inside a (locally) homogeneous isotropic elastic medium, rotational motions are thus exclusively exhibited by S-waves. Hence, by measuring ω with a rotational sensor, we can extract the curl of the wavefield and isolate the S-wave constituent of u, which reads within a medium as  
ω=12×u=12(yuzzuyzuxxuzxuyyux).
(11)
In anisotropic media, this separation no longer holds and even P-waves can have a rotational component (Pham et al., 2010).

Rotational motion at the free surface

At the free surface, the vertical components of the stress tensor vanish (σiz=0  for  i=x,y,z) and, therefore, also the strain tensor components εiz=0 for i=x,y (equation 3). Due to these constraints imposed by the free surface, equation 11 can be rewritten as (e.g., Robertsson and Curtis, 2002; Cochard et al., 2006)  
ωFS=(yuzxuz12(xuyyux)),
(12)
which relates rotation measurements at the free surface (denoted by the superscript FS) to the spatial derivatives of the displacement field. Note that no vertical derivatives appear in equation 12 because the free-surface condition allows us to convert vertical derivatives into horizontal derivatives (e.g., Robertsson and Curtis, 2002; Cochard et al., 2006). Equation 12 illustrates that measurements of the rotations around the two horizontal axes provide point estimates of the spatial wavefield gradients of the vertical displacement at the free surface, a quantity that is of relevance for wavefield reconstruction and imaging (e.g., Muyzert et al., 2012). Considering that combined displacement and displacement spatial gradient measurements allow reconstructing the wavefield away from the sensor location, rotational motions provide spatial information on a large spatial scale that would otherwise only be available through large receiver arrays.

Plane-wave analysis of rotational motion induced by body and surface waves

At the free surface, estimates of the rotational motion induced by body and surface waves can be deduced from acceleration measurements of monochromatic plane waves. Considering a wave traveling in the x-direction, the rotation rate around the y-axis and vertical acceleration for P-, SV-, and Rayleigh waves is (e.g., Li et al., 2004; Langston, 2006)  
ω˙y=pxau¨z,
(13)
where pxa is either the P-, S-, or Rayleigh-wave’s local apparent horizontal slowness in the x-direction. Note that when a P-wave impinges on the free surface under nonnormal incidence, a P-to-S-wave conversion occurs that gives rise to rotational motion (see Figure 3).
In the frequency domain, equation 13 reads as  
ω˙˜y=pxa(j2πf)u˙˜z,
(14)
where the tilde symbol marks variables in the frequency domain, j is the imaginary unit with property j2=1, and f is the frequency. From equation 14, we note that rotation amplitudes increase with increasing frequency and increasing slowness. These relations suggest that significant rotational energy can be expected in the typical seismic-exploration frequency range for similar translational-motion amplitudes compared with global seismology. Furthermore, slow-traveling surface waves are characterized by higher levels of rotational energy compared with body waves.
For SH and Love waves traveling along the x-direction, the induced rotational motion around the vertical axis at the free surface is given by (e.g., Igel et al., 2005; Lin et al., 2011)  
ω˙z=12pxbu¨y
(15)
with pxb being the local apparent horizontal S-wave slowness. Equations 1315 show that rotation rate and acceleration are related by the apparent horizontal slowness. Hence, amplitude ratios from collocated 6C measurements can be used to retrieve instantaneous slowness (velocity) information for isolated arrivals without using traveltimes, which has been exploited, for example, by Igel et al. (2005, 2007), Langston (2007a), Edme and Yuan (2016), and Wassermann et al. (2016). An example is shown in Figure 2. Bernauer et al. (2009) extend the analysis of translational and rotational ground motion amplitude ratios and derive finite-frequency kernels for a velocity tomography method that does not require traveltime measurements as an input (see also Fichtner and Igel, 2009).

In summary, a collocated point measurement of translational and rotational motions contains information that allows directly inferring on the subsurface seismic velocity structure and propagation direction. Thus, it is possible to reconstruct local, shallow velocity models in the vicinity of the receiver from amplitude information without using traveltimes.

Anisotropy

The curl operator separates divergence-free S-waves from curl-free P-waves inside a (locally) homogeneous elastic 3D isotropic medium (equation 10). This separation no longer holds in anisotropic media in which a quasi-P-wave also contributes to the curl and, thus, to the observation on a rotation sensor in an anisotropic medium. The analysis of rotational motions in anisotropic media remains unexplored except for the study by Pham et al. (2010) in which the effect of anisotropy on P-waves was investigated using the Christoffel equation. They conclude that, for local earthquakes and typical reservoir situations, quasi-P-rotation rates induced by anisotropy are significant, recordable, and can be used to constrain inverse problems on anisotropic parameters. Therefore, there is a potential to use rotational-motion observations to quantify anisotropic parameters, in particular, in downhole experiments (see the section on applications for vertical seismic profiling (VSP) examples and Horne, 2012).

INSTRUMENTATION: HOW TO MEASURE ROTATIONS

Recent approaches to develop portable instruments for the recording of rotational motions involve a variety of technical solutions. The most common techniques are liquid-based rotational motion sensors, magneto-hydrodynamic sensors, and fiber-optic gyroscopes (Figure 5). Liquid-based rotational motion sensors measure the differential movement in a liquid-filled torus by means of molecular electronic transfer transducers. These sensors were found to be suitable for strong motion measurements. Comprehensive overviews of liquid-based rotational motion sensors are given by Huang et al. (2013) and Egorov et al. (2015). Some of these sensors were successfully used in strong motion field campaigns (Lee et al., 2009; Nigbor et al., 2009; Wassermann et al., 2009), but they show high sensitivity to changes in ambient temperature (Bernauer et al., 2012).

Magneto-hydrodynamic sensors (e.g., Pierson et al., 2016) are another approach taking advantage of the inertia of a conductive fluid in a stable external magnetic field. The relative motion between the fluid and the magnetic field in the case of externally applied rotation generates a radial current that can be measured.

In a fiber-optic gyroscope, the difference between the optical path length of two counterpropagating laser beams in an optical fiber loop causes an interference pattern that is dependent on the applied rotation rate (the Sagnac effect, see Figure 1; Lefèvre, 2014). The most important advantage of this measurement principle is that it is based on massless photons. This makes it inherently insensitive to translational motion and guarantees the measurement of pure rotation. Several studies have already been conducted using fiber-optic gyroscopes for sensing rotational motion (Jaroszewicz et al., 2006, 2012; Velikoseltsev et al., 2012; Lindner et al., 2016). Recently, the first commercial fiber-optic gyroscope designed specifically for the needs of seismology has become available (BlueSeis3A by iXBlue, France; Bernauer et al., 2017).

An overview of sensor self-noise levels of different portable rotational ground motion sensors is shown in Figure 6. The solid and dashed lines represent half-octave integrated root power spectral densities. This time-domain representation of sensor self-noise is commonly used in so-called operating range diagrams (ORDs) (Evans et al., 2010). The self-noise spectra of BlueSeis3A (BS1), the fiber-optic gyroscope LCG-Demonstrator (LCG by Litef, Germany), and the liquid-based sensor R2 (Eentec, 2017) were recorded with the sensors being installed on an auxiliary monument of the seismic station at the Geophysical Observatory in Fürstenfeldbruck, Germany. For estimating the self-noise level of the seismic magneto-hydrodynamic prototype sensor ATA-smhd (ATA by Applied Technologies Associates), we took a 30 min lasting record with the sensor being installed in Oklahoma, United States (data available from the IRIS-DMC database, 2017; start time: 2017-05-07, 04:27:18). Note that different site conditions may bias the comparison of self-noise levels for different sensors. Our estimate for the self-noise level of the ATA-smhd is very close to the one reported by Pierson et al. (2016), which was recorded at the Albuquerque Seismological Laboratory underground vault. This circumstance supports the assumption that the contribution of site effects to our comparison of self-noise levels is negligible.

The dot- and star-shaped symbols in Figure 6 represent rotation rate amplitudes and frequency ranges recorded from sweep (10–120 Hz) and blast events during an active-source seismic exploration campaign in Glattfelden, Switzerland. In the ORD plot, the source examples are abstracted in the following way: First, the peak amplitude is extracted from the event. Second, the power-spectral density (PSD) is computed from a record that contains the event. The half-peak-amplitude width of the PSD is used as the “bandwidth” of the event. Finally, the peak amplitude scaled by 2 is used as the “root-mean-squared” amplitude of the event. For more details on ORD representation, the reader is referred to Evans et al. (2010).

Other techniques proposed to measure rotations are based on microelectromechanical system technology (e.g., D'Alessandro and D'Anna, 2014; Liu and Pike, 2016), on mounting conventional geophones on a fixed frame and use differencing of the geophone recordings to estimate rotation rates around all three coordinate axes (i.e., rotaphone; Brokešová and Málek, 2010), or on instruments including a rotational transducer to measure the instantaneous rate of particle rotation, so-called “rotational geophones” (Cowles, 1981).

Array-derived rotations

Because reliable rotational-motion instruments have not been available until recently, arrays of conventional single or multicomponent sensors have been used in several studies to estimate rotational motion by computing the horizontal spatial gradients of the recorded seismic wavefield (e.g., Spudich et al., 1995; Huang, 2003; Suryanto et al., 2006; Langston, 2007a, 2007b, 2007c; Liang and Langston, 2009; Spudich and Fletcher, 2009; Sollberger et al., 2016b). The rotational motion at the free surface around all three axes can be estimated using equation 12. Furthermore, arrays of 3C stations enable the estimation of the full displacement gradient tensor and, hence, strain and dilatation (see also the “Theoretical foundations” section) (e.g., Spudich and Fletcher, 2009; Donner et al., 2017a). Estimating spatial gradients and rotations using finite-difference techniques is based on the assumption that the subsurface is homogeneous and deformation is linear over the array area. In addition, uniform receiver-to-ground coupling and receiver characteristics (e.g., transfer functions, noise) are assumed to be uniform within the array.

Computing the spatial gradients by simple differencing of the recordings of neighboring stations may be particularly prone to noise. Therefore, inversion schemes are usually used to stabilize the spatial gradient estimation (Liang and Langston, 2009; Spudich and Fletcher, 2009). Even when using inversion schemes, the accuracy of the spatial gradient estimates depends on the relation between the apparent wavelength and the arrival azimuth of a seismic wave and the array geometry (Schmelzbach et al., 2016), introducing a pseudo-anisotropy in the gradient estimates (e.g., de Ridder and Curtis, 2017).

We carried out a field experiment on the ETH Zurich campus (Switzerland) with the goal of comparing array-derived rotation estimates with point measurements using rotational seismometers (METR-03). In the first experiment, spatial gradients were estimated using a regular cross-shaped layout of five conventional 3C geophones with four different rotational seismometers (Figure 7a). Because the instrument response of the geophones and rotational seismometers were not sufficiently well-known, a second calibration experiment was carried out in which two circles of 32 3C geophones each were laid out (Figure 7b). The geophone-array data from the calibration experiment allowed extraction of high-quality spatial gradient estimates due to the large number of receivers and the uniform angular coverage. These rotational-rate estimates were then used to design a filter to match the rotational-sensor data. The rotational sensor data from the first experiment filtered with the inverse transfer function from the calibration experiment were then compared with the array-derived rotations from the geophone-cross array showing an overall close match between the two measurements (Figure 7d). The somewhat increased mismatch at large offsets could be due to the fact that the match filter was designed for data recorded with offset <20  m. The quality of the array-derived gradient estimates depends on the apparent wavelength of an arrival, and it is based on the assumption of plane waves (Schmelzbach et al., 2016). Considering that the calibration data were recorded relatively close to the source in relation to the dominant wavelengths and that the dominant wave types change as a function of offset could explain the increasing mismatch with increasing offset observed in Figure 7d. In comparison with the vertical-component recording (Figure 7c), the rotation rate around the transverse axis (Figure 7d) shows enhanced surface and S-wave arrivals, whereas the P-wave first arrivals are significantly reduced in amplitude. This experiment confirms that it is indeed possible to estimate reliably rotational motion from arrays of conventional geophones for exploration-scale settings, which has also been shown by Muyzert et al. (2012) in a similar experiment.

Limitations of array-derived rotations

Processing recordings from different receivers with some areal distribution can further be affected by near-receiver variations in the subsurface, different receiver-to-ground coupling, and variations in the receiver characteristics. To overcome these issues, Sollberger et al. (2015) propose a receiver-coupling correction scheme. An example of the receiver-coupling correction procedure is illustrated in Figure 8. Note how the receiver coupling corrections improve the continuity of the spatial gradient estimates, yielding more reliable array-derived local apparent phase velocity estimates (see equation 13). For example, the apparent velocity estimated for each array independently better matches the apparent velocity of approximately 270  m/s estimated by “classic” linear-moveout analysis for the Rayleigh wave (marked with R). Further issues with array-derived rotation estimates can be caused by strain-rotation coupling from near-receiver heterogeneity leading to variations in rotational motions experienced by the individual sensors of an array (van Driel et al., 2012).

POTENTIAL APPLICATIONS

Applications of rotational-motion measurements and 6C data in exploration seismology can roughly be grouped into three main categories with the rotation data being interpreted as (1) direct measurements of the S-wave component of the seismic wavefield, (2) estimates of the spatial wavefield gradient, and (3) additional seismic observable extending the set of observables for deterministic and statistical wavefield analyses. In the following section, we discuss a series of applications grouped into exploration and global seismology examples and the potential links between the two fields.

Applications in exploration seismology

Improved S-wave identification

The wave equation expressed in terms of potentials in isotropic (locally) homogeneous media consists of a curl-free P-wave potential and a divergence-free S-wave potential (see equation 9). Hence, applying the curl operator to this formulation of the wave equation will isolate the S-wave potential (see equation 10). Because measuring rotations can be understood as measuring the curl of the wavefield, rotational measurements provide a direct estimate of the S-wave component of the wavefield. Sollberger et al. (2016a), for example, show that array-derived rotation estimates from a dense high-resolution multicomponent profile can be used to isolate the direct S-wave, highlighting the potential of rotational data to improve S-wave exploration. Array-derived rotations may be estimated from vintage pure vertical-component translational data (e.g., geophone data proving particle velocity data), opening up new ways to reevaluate existing data.

The analysis of the Apollo 17 data by Sollberger et al. (2016b) is an impressive example of reevaluating vintage data with rotational seismology tools and resulted in the first near-surface S-wave profile of the moon. This study demonstrates that rotational data can provide critical information in cases with very limited site access and resultant very small data sets of exceptionally complex character. In 1972, seismic data generated with explosives at eight locations were recorded with four geophones in a star-like arrangement during the Apollo 17 mission (Figure 9). Although the original analysis of the data allowed only to establish a P-wave velocity profile at the landing site, the computation of array-derived rotation estimates enabled, for the first time, to identify S-wave arrivals and to construct an S-wave velocity profile of the shallow lunar crust from a total of only 32 seismic traces.

Spatial-gradient estimates and wavefield reconstruction

At the free surface, the rotational motion around the two horizontal axes corresponds to the horizontal gradient of the vertical particle velocity component (see equation 12). Point measurements of the spatial gradient open up new avenues to address some of the long-standing problems in seismic exploration such as the requirement of dense spatial sampling to avoid spatial aliasing of surface waves. Collocated measurements of the seismic wavefield and its spatial gradient allow relaxing the sampling requirements and reconstructing the wavefield from data that are spatially aliased according to Shannon’s sampling criterion. Such wavefield reconstruction techniques are successfully used with multicomponent marine data (i.e., measurements of the pressure and its spatial gradient in three directions; Robertsson et al., 2008; Vassallo et al., 2010). Muyzert et al. (2012) discuss an application to land seismic data using combined particle velocity and rotational measurements to reconstruct aliased land seismic data. A recent extension of the wavefield-reconstruction concept by inverting for the wave parameters and the near-surface properties was presented by de Ridder and Maddison (2017). Maeda et al. (2016) use seismic gradiometry to reconstruct the seismic wavefield using earthquake recordings of the Japanese Hi-net seismography network.

Local velocity estimations

Interesting applications arise from the so-called seismic gradiometry technique proposed by Langston (2007a, 2007b, 2007c) and Liang and Langston (2009) that allow estimating the local and instantaneous slowness of an isolated plane wave by combining measurements of the spatial and temporal gradient of a wavefield (see equations 13 and 15). Liu and Holt (2015) apply this technique to estimate phase velocities for USArray data.

An application of the instantaneous local slowness estimation is shown, for example, in Figure 9 for the Apollo 17 array data. Edme and Yuan (2016) present a workflow to estimate local surface-wave dispersion curves based on instantaneous local slowness estimates of ambient noise using array-derived rotation data. The study by Edme and Yuan (2016) using small arrays of three geophones in 1.5 crossline and 2 m inline separation highlights that spatial gradient and rotational data can potentially replace the multichannel analysis of surface waves technique, which requires long receiver arrays. Other spatial-wavefield gradient-based approaches for estimating local subsurface parameters have been proposed by de Ridder and Biondi (2015) using a wave-equation-based technique to estimate phase velocities from short time windows of seismic noise recordings; de Ridder and Maddison (2017) propose a full-wavefield inversion techniques that exploits constraints on the spatial wavefield gradients.

A novel tomographic scheme was suggested by Fichtner and Igel (2009) and Bernauer et al. (2009, 2012) as well as a generalization of this concept to optimal observables by Bernauer et al. (2014). Inspired by the simple relation between transverse acceleration and rotation rate (equation 15), the concept of “apparent velocity” as a wavefield observable was defined as the ratio of amplitudes of acceleration (or velocity) and rotation rate (or rotation angle). Applying the adjoint method to such observables allows the calculation of the corresponding sensitivity kernels. It is intuitive that the sensitivity with respect to structure near the source vanishes and only sensitivity adjacent to the receiver position becomes important. This allows the development of an inversion scheme based on the (e.g., frequency-dependent) apparent velocity observations without ever making use of traveltimes (Bernauer et al., 2009). Such schemes might prove useful to derive local subsurface velocity models in the absence of a seismic network (e.g., small-scale site effect characterization or subsurface velocity estimation on rocky planetary objects).

Wavefield separation

Robertsson and Muyzert (1999) propose using 3D (volumetric) receiver layouts to estimate all spatial derivatives to compute curl and divergence of the measured seismic wavefield to separate the wavefield into its P- and S-wave constituent. Using surface-based multicomponent receivers only, Woelz et al. (2009) achieve a wavefield separation of field data with approximated curl and divergence operators. Robertsson and Curtis (2002) show that dense-array measurements taken at the free surface are sufficient to estimate the horizontal and the vertical derivatives due to the free-surface boundary conditions (see equation 12). They then derive a set of wavefield separation filters to estimate the pure upgoing (incident) P- and S-wave constituents at the free surface, which require estimates of the spatial wavefield derivatives. The derivation involves a Taylor expansion with the resultant first-order terms containing spatial gradients of the wavefield.

Van Renterghem et al. (2016) present variants of these gradient-based wavefield separation filters to estimate the pure upgoing (incident) and/or P/S separated particle velocity components of land seismic data. Van Renterghem et al. (2017, 2018) extend the approach to ocean-bottom seismic data to provide, for example, the upgoing P- and S-wave constituents with receiver-side water-layer multiples (i.e., downgoing wavefield) removed. For the land-seismic case, the gradient-based filters to isolate the upgoing S-wavefield read as (e.g., Van Renterghem et al., 2017, 2018) 
vxSU12(vx1jν2βxvz),
(16)
where j is the imaginary unit with property j2=1, ν=2πf is the angular frequency with f marking frequency, and β is the S-wave velocity at the receiver location. Considering that xvz corresponds to the rotation around the y-axis that could be measured with a rotational sensor, equation 16 illustrates the potential of collocated velocity and rotation measurements for wavefield separation to enable “true” amplitude S-wave imaging. Figure 10 shows a synthetic and field-data example of P/S separation of land-seismic data.

Single-station wavefield characterization and surface-wave suppression

Taking a different viewpoint, rotational measurements provide additional data to extend the set of seismic observables for statistical wavefield characterization and processing. For example, Barak et al. (2014) present a surface-wave suppression scheme that is based on manually extracting the 6C statistical characteristics of surface waves in a sample time window and then use this so-called “wave signature” to remove surface waves in other time windows. Sollberger et al. (2018) derive analytical expressions for these wave signatures by providing 6C polarization models for Rayleigh waves, which allows us to identify and suppress ground roll in a single-station approach.

Based on the observation that recordings of the rotational motion around the horizontal axes are slowness-scaled versions of the vertical-component wavefield with surface-wave amplitudes amplified relative to reflections, Edme et al. (2013) propose a technique to suppress side-scattered surface waves by adaptively subtracting time-integrated rotational data from vertical-component recordings. Following field tests by Edme and Muyzert (2013), they use arrays of two vertical geophones to determine the crossline rotational motion based on finite differencing.

P-S separation in VSP experiments

In isotropic media, the curl operators perfectly separate the divergence-free shear motion from the curl-free compressional motion (equation 10). This is illustrated in a synthetic seismogram section of a VSP geometry in Figure 11 modeled on the calculations by Horne (2003) but with 6C rather than 3C downhole recordings. The velocity model consists of four isotropic layers. The source distance to the borehole is 500 m, the receivers are located at a depth range of 405–1995 m (vertical receiver separation of 15 m), and the dominant frequency is 10 Hz. The specific details of the velocity model are not relevant here. What is important is the fact that in this case the rotation records allow a direct separation of the S-wavefield without processing. An explosion source only generates P- and SV-waves with the SV-waves appearing as the only signal on the rotation component perpendicular to the sagittal plane (i.e., the ωy-component). Even though it is illustrated here on a simple subsurface model, there are numerous advantages of accessing shear-velocity information directly. In addition, the amplitude ratios between translational and rotational motion measurements, as discussed above, can be used to access apparent phase-velocity information (see also Horne, 2012).

As indicated by Pham et al. (2010), the isolation of S-waves by computing or measuring the curl of the wavefield no longer holds in anisotropic media because P-waves are no longer curl-free. A rotation sensor would measure a signal for a quasi-P-wave recorded in an anisotropic structure. In return, this effect could potentially offer additional constraints on anisotropic parameters measured in downhole experiments.

Vector infidelity correction of ocean-bottom-seismometer data

Ocean-bottom seismometers (OBSs) are subject to translational and rotational motions caused by seafloor currents as well as induced by the arrival of seismic waves (i.e., rocking sensor motions). In addition, OBSs are almost always placed on unconsolidated sediments (Sutton et al., 1981) resulting in severe resonances in the OBS-seafloor system. As a result, the horizontal OBS components often appear noisy and OBS multicomponent recordings suffer from vector infidelity that needs to be corrected for (Garmany, 1984; Dellinger et al., 2001). Part of the horizontal component contamination stems from rotational (tilt) motions of the OBS system (Gaiser, 2007). Lindner et al. (2016) recently demonstrate that the inclusion of rotational seismometers in ocean-bottom nodes allows to remove tilt contaminations from the horizontal translational-motion components of OBSs leading to improved S-wave imaging.

Seismological applications

Rotational-motion measurements have so far mainly been used in seismological applications to improve the localization of earthquakes and to study their source parameters. These applications are of direct relevance to investigate induced seismicity in hydrocarbon or geothermal reservoirs. Furthermore, global-seismology techniques using rotational data to characterize surface waves can be of interest for near-surface geophysical investigations, where conventional surface-wave analyses are important tools.

Back-azimuth estimation

To obtain the transverse and radial components of the wavefield needed to determine the apparent slowness of an arrival, the back azimuth ϕ of the incoming signal is needed. To solve this problem, Igel et al. (2007), Hadziioannou et al. (2012), and Salvermoser et al. (2017) (with accessible Python codes) apply a systematic grid search for ϕ by stepwise rotation of the horizontal components of acceleration to find the direction of maximum correlation with the vertical rotation rate (see Figures 2 and 12). The corresponding phase velocity is then determined. Wassermann et al. (2016) improve this method by performing an orthogonal distance regression on the simultaneous estimate of ϕ and phase velocity c using the equation  
2cω.z=sin(ϕ)u¨NScos(ϕ)u¨EW,
(17)
where wz is the vertical component of rotation and u¨NS and u¨EW are the translational acceleration in the north–south and east–west directions, respectively. Equation 17 allows incorporating errors in either of the measurements.

Earthquake source location

In principle, it is possible to determine an earthquake epicenter (distance and back azimuth) and origin time from a measurement of the three translational components of a single station assuming that the velocity model is known. However, this approach comes with large uncertainties, and depth estimations are only possible in the case that rare depth phases (e.g., pP, sS phases) exist in the data (e.g., Khan et al., 2016; Böse et al., 2017).

Li and van der Baan (2017b) suggest an approach to automatically detect event hypocenter and origin time incorporating the additional information on the spatial gradient of the particle velocity wavefields. They time reversed the recordings, forward propagate them through the medium, and finally stack all back-projected images. By applying a focusing criterion to the wavefield at each time slice, they solve the localization problem. However, their approach is involved and assumes a known velocity model.

Donner et al. (2017c) propose a different approach to determine event hypocenter and origin time from a collocated 6C measurement in a more straightforward way and without any prior knowledge on the velocity model. This approach includes the following steps: The origin time is determined in a conventional way from the arrival times of P- and S-waves. Using the ROLODE method (equation 17; Wassermann et al., 2016), the back azimuth is determined. Using this back azimuth, the translational acceleration is correlated with the vertical rotation rate in different frequency bands to obtain a phase velocity at each frequency resulting in a dispersion curve. From the dispersion curve, a velocity model can be derived, which is then used to estimate the distance of the signal. In combination with a ray tracer, the determined velocity model also allows us to obtain depth and distance estimates. Figure 12 shows an example localization derived from combined recordings of three translational components and the vertical rotation component of the Mw 6.2 central Italy earthquake of 24 August 2016. Being able to locate events precisely from a point measurement is of high interest for planetary seismology applications as well (Donner et al., 2017b).

Inversion for the earthquake source parameters

Seismic point source moment tensors are a general mathematical description of an earthquake as a point source (Jost and Herrmann, 1989; Aki and Richards, 2002). They provide information on the orientation of faulting (i.e., tectonic part) of an earthquake as well as on nontectonic parts (e.g., volume changes, fracturing/cracking, volcanic/geothermal activity), centroid depth, and earthquake strength. The seismic point source moment tensor is a symmetric tensor of second order with six independent components. According to Aki and Richards (2002), the relation between ground motion and moment tensor is given by  
un(x,t)=Mkj·jGnk(x,tτ),
(18)
where un is the displacement motion observed at a point x to a time t relative to the source time τ. The terms Mkj are the components of the moment tensor and jGnk are the spatial derivatives of the displacement Green’s functions (the Einstein summation convention applies). Rupture length, width, rupture velocity, and direction are additional inversion parameters when determining the kinematics of earthquake sources (Bernauer et al., 2014).

Although the theory of waveform inversion for source solutions is well established, the successful source solution retrieval is often hampered by several difficulties (e.g., Donner et al., 2016), especially in the local and regional distance range. Theoretical studies comparing inversion results using only 3C of translational motion data with data from only half the number of stations but recording 6C (i.e., 3C of translational and 3C of rotational motion) have shown that the inclusion of rotations leads to improved point and kinematic source solutions (Bernauer et al., 2014; Donner et al., 2016; Reinwald et al., 2016). In particular, the resolution of the depth-dependent components of the moment tensor and the centroid depth can be improved drastically. Donner et al. (2017a) have shown that the full moment tensor could be obtained from a single station with collocated measurement of translation and rotation.

Focusing on the double-couple part of the source, Cochard et al. (2006) derive the rotational radiation pattern that corresponds to the curl of the S-wave far-field radiation pattern. The rotational radiation pattern looks different in orientation compared with the translational one. Thus, the amplitude ratios of the different combinations of ground motion components contribute strongly to the better resolved mechanism.

The improved moment tensor results due to the inclusion of rotational motion measurements may also be helpful in various aspects of induced seismicity, for example, in hazardous cases such as the Rudna (Poland) mine collapse on 19 March 2013 (Rudziński et al. 2016) as well as in unclear cases such as the 2012 September-October seismic sequence offshore Spain, which is suspected to be induced by human activities (Cesca et al., 2014). Dahm et al. (2015) have set up a quantitative probabilistic approach to discriminate between natural and human-induced/triggered earthquakes. They state that the uncertainties in location and mechanisms influence the discrimination and need to be reduced as much as possible to come to reliable conclusions. The additional constraints added by the integration of rotational motion into inversion for the seismic source help to improve the inversion results and thus to better understand underlying deformation mechanisms.

Next to discriminating induced from natural seismicity, the analysis of microseismicity (ML<3.0) — induced or natural — is still a challenge. Furthermore, small events require a broader frequency range for waveform inversion compared with large events. Consequently, the successful inversion for moment tensors strongly depends on the resolution and reliability of the underlying velocity model. Especially, the resolution of non-double-couple parts suffers from inaccuracies in the velocity model (Sileny, 2004). A topic of current research is whether including the 3C of rotation, a 3D velocity model, or the usage of both will improve the resolvability of the moment tensor from microseismicity (Donner et al., 2017b).

CONCLUSION

Preliminary studies using 6C data highlight the great potential of direct measurements of the three rotational components induced by seismic waves. It can be expected that the rapid advances in rotational sensor technology that have been observed over the past few years will further progress and enable accurate 6C measurements with affordable, portable, and robust devices in the near future. Potential applications of collocated translational and rotational (6C) data can be divided into (1) exploiting the direct link of rotational motion to S-waves and surface waves for, for example, enhanced S-wave processing and ground-roll suppression; (2) using rotational-motion data at the free surface as spatial wavefield gradient estimates for local-slowness estimation, wavefield reconstruction, and wavefield separation; and (3) extracting near-receiver material properties by the analysis of relative amplitudes of rotational and translational data. One particularly attractive aspect of 6C measurements is that 6C data enable analyses with single-station recordings that conventionally require large receiver arrays. Therefore, 6C measurements may be particularly interesting for applications in environments in which only few and sparsely distributed seismic sensors can be installed such as in boreholes, at the ocean bottom, or for planetary seismology missions.

From the research directions in rotational seismology that we identified in this paper, we expect significant advances in the coming years in the following fields of exploration seismology with relevance for global seismology as well:

  1. 1)

    Full elastic wavefield decomposition into all wave modes: Because the 6C motion is characteristic for each wave mode, arrivals can be automatically identified and separated from a direct measurement of rotations from a single station.

  2. 2)

    S-wave and converted wave imaging: Due to the inherent link between rotational measurements and S-waves, we expect S-wave imaging to benefit the most from 6C data. The development of divergence sensors could further enhance the separation of P- and S-wavefields.

  3. 3)

    Near-surface characterization and velocity model building by tomography without traveltimes: Rotational data and seismic gradiometry can lead to improved near-surface velocity models that can help to suppress distorting near-surface effects from seismic images by improved static corrections, especially for S-waves.

  4. 4)

    Wavefield reconstruction: The inherent link of rotational-motion measurements at the free surface to the spatial gradient of the seismic wavefield bears the potential for powerful wavefield reconstruction techniques that could allow rethinking current high-channel count acquisition strategies.

  5. 5)

    Anisotropic medium characterization: The three additional, independent observations that 6C data provide compared with standard 3C data can help to better constrain multiparameter inversion problems, such as the estimation of anisotropy parameters. For example, P-waves can exhibit a rotational component when propagating through an anisotropic medium.

  6. 6)

    Microseismic monitoring: Given the demonstrated enhancement of earthquake moment tensor inversions by including rotational data, we expect that 6C data will help to better constrain the location and focal mechanisms of microseismic events. This could prove valuable to monitor hydraulic fracturing experiments and reservoirs.

  7. 7)

    Rotational seismic sources: Arrays of closely spaced, active seismic sources theoretically allow to simulate rotational and dilatational sources, which could further improve wavefield separation techniques.

  8. 8)

    Vector infidelity corrections for OBS and ocean-bottom cable-data: Tilt contamination of horizontal component OBS data can be suppressed by additionally taking into account rotational data, which could enhance ocean-bottom multicomponent exploration.

The emerging availability of portable rotation sensors fit for broadband observations in seismology and seismic exploration could bring the many potential uses of 6C data from the synthetic realm and small-scale field tests to widespread applications.

ACKNOWLEDGMENTS

The authors are thankful for the support by the CARNEVAL industry consortium (Nagra, OMV, Schlumberger Gould Research). They are also very thankful to the company iXblue for the cooperation in developing a portable, broadband rotation sensor. The research of this study was also partly supported by the European Research Council (advanced grant to HI: ROMY, number: 339991) and the Swiss National Science Foundation (grant 200021_156996). We acknowledge support from discussions within TIDES COST Action ES1401. We appreciate the constructive comments on an earlier draft of this manuscript by the associate editor and the journal reviewers O. Barak, S. Horne, and M. van der Baan.

DERIVATION OF THE TOTAL ROTATION OF A LINE SEGMENT

The derivation of the total rotation angle ω of a line segment δx presented here (equation 6) follows Segall (2010). The total rotation angle between the undeformed and deformed line segment δx and δξ is given by the cross-product (Figure 4)  
ω=δx×δξ|δx|2.
(A-1)
Here, we assume infinitesimally small deformation |δx||δξ| and small angles. Then, the magnitude of ω is  
|ω|=|δx||δξ||δx|2sinωω.
(A-2)
Equation A-1 reads in indicial notion as  
ωi=eijkδxj×δξkδxnδxn,
(A-3)
where eijk is the Levi-Civita symbol (permutation tensor), and we used the relation (a×b)i=eijkajbk as well as the Einstein’s summation convention. We can express the deformed state δξi using the displacement gradient  
δξk=(muk+δkm)δxm,
(A-4)
where δkm denotes the Kronecker delta. Inserting equation A-4 into equation A-3 yields  
δxnδxnωi=eijkδxj(muk+δkm)δxm,=eijkmukδxjδxm+eijkδxjδxk,=eijk(εkm+θkm)δxjδxm.
(A-5)
For the last step, we used the expansion of the displacement gradient tensor into the sum of a strain (εkm) and rotation tensor (θkm) (see equation 2). Furthermore, note that eijkδxjδxk corresponds to the cross-product of δx with itself and is zero. Defining a unit vector δx^j (δx^=δx/|δx|), we can express the total rotation as  
ωi=eijk(εkm+θkm)δx^jδx^m.
(A-6)
We can replace the rotational term of equation A-6 using ekmnΘn=θkm:  
θkmδx^jδx^m=eijkekmnΘnδx^jδx^m,=(δimδjnδinδjm)Θnδx^jδx^m,=ΘiΘjδx^jδx^i.
(A-7)
For these steps, we used the relation between the permutation tensor and the Kronecker delta as well as the fact that δx^mδx^m=1. This leads to the general expression of the total rotation of a line segment  
ω=δx^×(ε·δx^)+Θ(Θ·δx^)δx^,ωi=eijkεkmδx^jδx^m+ΘiΘjδx^jδx^i,
(A-8)
which corresponds to equation 6. The total rotation ω has a strain component (first term) and two rigid-body rotation terms (second and third term). The strain component of ω corresponds to the cross-product of the unit vector δx^ and the relative displacement dui=εijδx^j. Depending on the angle between the line segment and the rigid-body rotation vector, Θ contributes between zero (i.e., the rotation vector is parallel to the line segment) and Θ (i.e., the rotation vector is normal to the line segment) to ω.
Freely available online through the SEG open-access option.