We use surface wave measurements to reveal anisotropy as a function of depth within the Juan de Fuca and Gorda plate system. Using a two-plane wave method, we measure phase velocity and azimuthal anisotropy of fundamental mode Rayleigh waves, solving for anisotropic shear velocity. These surface wave measurements are jointly inverted with constraints from shear wave splitting studies using a Markov chain approach. 

We use Rayleigh wave data in the period range 20 to 125 s from Bell et al. (2016), augmented by data from the fourth year of the CI deployment, which greatly improves resolution in the southern half of the study area including Gorda, the Blanco fracture zone region, and the southern part of Juan de Fuca plate.  In brief, we select shallow earthquake sources in the epicentral distance range 30–120° with Ms > 5.5, correct the records for instrument response to displacement, and then remove tilt and compliance noise from the vertical component using the transfer function techniques of Crawford and Webb (2000) and Bell et al. (2015).  The records are then filtered with a set of narrow band-pass filters and windowed to limit contamination with higher modes and scattered waves (Figures S1 and S2).  If the filtered records from an event show evidence for excessive multi-pathing or signal-to-noise ratios < 10, those periods or stations are eliminated.  In total, we use data from 107 events, well distributed in azimuth, with  the number of seismograms ranging from 2232 at 20 s to 2910 at 50 s to 1750 at 125 s period.  

To obtain the phase velocities as a function of azimuth, period, and location, we employed the two-plane-wave method with two-dimensional finite frequency kernels (Forsyth and Li, 2005; Yang and Forsyth, 2006; Rau and Forsyth, 2011; Bell et al., 2016).  Beginning with the predictions of phase velocities for an isotropic starting model, we solved for corrections to the model linearly interpolated between grid points with a 0.5° spacing.  At each grid point, we assumed the phase velocity is a 2θ function of frequency and direction. We neglected 4θ terms because they are expected to be small for structure with orthorhombic symmetry (Smith and Dahlen, 1972), as has been demonstrated for young seafloor near the East Pacific Rise (Weeraratne et al., 2007).  

The starting isotropic phase velocity model was similar to that of Bell et al. (2016), but represented the result of iteratively inverting body and surface wave observations (VanderBeek et al., 2018), incorporating variations in water depth, sediment and crustal thickness, and age of the seafloor.  To obtain stable results that can be compared to the splitting results of Bodmer et al. (2015), we required that sine and cosine coefficients of phase velocity be uniform at each period for all grid points that fall into each of the four oceanic subregions illustrated in Figure 1b (i.e., describing structure that varies only in depth within each region). To avoid bias from neglecting structure outside the subregions, we introduced six other subregions surrounding the array in which anisotropy was also assumed to be uniform at each period. However, the anisotropic terms in those regions are poorly constrained because the limited azimuthal range possible for paths approaching the array from the outside leads to greater tradeoffs with the isotropic variations; we ignore those results in the rest of this paper.

We inverted the sine and cosine terms for perturbations to the shear velocity structure in many thin layers extending down to 410 km, with minimum length and minimum curvature regularization.   The results, in terms of one particular solution to each of these non-unique problems, are shown by the dotted lines in Figures 3c,d for the Juan de Fuca and Gorda subregions. Resolution analysis reveals that the resolution matrices are reasonably compact with rank ~ 3.5, indicating that there are approximately 3.5 independent pieces of information about the anisotropic vertical structure.  The average strength and direction of the fast axis can be resolved from Moho depth down to about 40 km, from 40 to 90 km, and from 90 to 210 km, with a little information about deeper structure.  We expect anisotropy beneath young seafloor to develop primarily from shear strain in the evolving low-velocity zone, so in what follows we assume that it is negligible deeper than 210 km (Figure 3) and adopt the averages and the covariances of the three resolved layers as the starting model for a joint inversion with the splitting data.

A joint inversion needs to account for the uncertainties in each method, and for the strong non-linearity of shear wave splitting in the case of multiple layers with different fast directions. To address this issue, we used a Monte Carlo Markov chain approach to simultaneously fit both data types through forward modelling.

Proposed models were created by normally perturbing sine and cosine anisotropic terms from the most recently accepted model (one pair of terms per layer). These models were used to forward-model SKS splitting by propagating synthetic pulses through each anisotropic model, as follows: For each layer, we calculated the full anisotropic tensor by assuming orthorhombic A-type fabric using a scaled version of the tensor computed for olivine by Anderson (1989), with horizontal a- and b-axes. Synthetic input waveforms comprised a 5 second Gaussian pulse on the radial component, and zero energy on the transverse component, with random white noise then added for a signal-to-noise ratio of 5 on the radial component. The synthetic pulses were propagated through layered model structure using source parameters and back-azimuths for the 16 most commonly observed earthquakes in the dataset employed by Bodmer et al. (2015), computing ray parameters using the 1D IASP91 model. We solved the Christoffel equations for fast and slow S-wave propagation velocities and polarizations in each layer in a ray-based coordinate system. Taking each layer in series, seismograms were rotated into fast/slow polarisation and delayed by the appropriate splitting time. Finally we output virtual north and east components, corresponding to an ‘observation’ at the top of the model. 

The resultant synthetic data were ‘un-split’ using the average splitting parameters observed in the Bodmer study for each of the study regions (Figure 1), rotating and shifting pulses according to the (single-layer) splitting time (δt) and fast azimuth (φ) measured in that study. In reference to the “Minimum Energy” method of Silver & Chan (1991), the transverse component energy, T, was used to define a misfit term that dictated the likelihood of this data having been observed given the proposed model (see paper). 

Proposed models were accepted (replacing the current model) or rejected using a Metropolis-Hastings acceptance criterion (Metropolis & Rosenbluth, 1953; Hastings, 1970). We ran 14 individual chains with 100,000 inversions per chain, saving the current model every 100 iterations; the concatenated ensemble of saved models yielded the empirical posterior model distribution. We conducted L-testing to determine weighting values that required the inversion to simultaneously fit both data types well.