The seismic signature of lithospheric deformation beneath eastern North America due to Grenville and Appalachian orogenesis

The formation of mountain belts at convergent plate boundaries is one of the fundamental processes in plate tectonics. Despite the importance of orogenesis to the plate tectonic system, our understanding of how deformation is accommodated in the crust and mantle lithosphere during orogenesis remains poorly understood. This is true both in present-day orogenic systems and in regions that have been affected by past orogenic events. Lithospheric deformation during orogenesis, while poorly understood, is crucial for our understanding of the evolution of topography, the partitioning of strain in collisional settings, and the evolution and modification of continental margins such as eastern North America. Furthermore, lithospheric deformation reflects the rheology of the crust and mantle lithosphere, which remain poorly understood despite extensive study.

One type of observation that can inform our view of deformation in the lower crust and upper mantle is seismic anisotropy, or the directional dependence of seismic wave propagation. This is because there is a direct, albeit complicated, link between strain and the resulting seismic anisotropy via the crystallographic preferred orientation (CPO) of anisotropic minerals. In the upper mantle, seismic anisotropy primarily reflects the CPO of olivine (e.g., Karato et al., 2008), while in the lower crust there might be contributions from amphiboles, micas, and/or quartz (e.g., Ward et al., 2012; Liu et al., 2015; Ko and Jung, 2015). There are many published observations of seismic anisotropy within continental settings in general (e.g., Fouch and Rondenay, 2006; Tommasi and Vauchez, 2015), and beneath eastern North America in particular (e.g., Barruol et al., 1997; Levin et al., 1999; Wagner et al., 2012; Long et al., 2010, 2016; Yang et al., 2017). A challenge in their interpretation, however, is that it is often difficult to obtain good constraints on the depth extent of anisotropy. For example, the splitting or birefringence of SKS waves (that is, waves that travel as shear waves through the mantle but as compressional waves through the core), perhaps the most common method for studying continental anisotropy, is a path-integrated measurement involving nearly vertically propagating shear waves, so its depth resolution is poor (e.g., Long and Silver, 2009). Complementary constraints on anisotropy can be obtained from surface wave dispersion analysis (e.g., Deschamps et al., 2008; Yuan and Romanowicz, 2010) and by combining different types of seismic data (e.g., Yuan and Levin, 2014; Bodin et al., 2016).

One analysis technique that can yield good resolution of the depth distribution of anisotropy within the lithosphere is anisotropic receiver function analysis, which can reveal sharp contrasts in anisotropic structure with depth (e.g., Levin and Park, 1997, 1998; Bostock, 1998; Frederiksen and Bostock, 2000; Liu and Park, 2017). This technique has recently been applied to data from continental settings; in particular, there are several recent studies that have applied it to study both crustal anisotropy (e.g., Porter et al., 2011; Schulte-Pelkum and Mahan, 2014a; Liu et al., 2015) and anisotropy within the lithospheric mantle (e.g., Yuan and Levin, 2014; Wirth and Long, 2014; Ford et al., 2016). This analysis strategy, which looks for P to SH wave conversions at a dipping and/or anisotropic interface, can provide unambiguous evidence for anisotropy (usually under the assumption of hexagonal symmetry) and constrains the timing of converted waves (and thus the likely depth of the interfaces). A disadvantage of the technique, however, is that the full anisotropic geometry in each layer across an interface cannot be easily deduced from the observations. Detailed forward modeling can identify plausible anisotropic models (e.g., Wirth and Long, 2012, 2014; McCormack et al., 2013), but forward modeling is generally nonunique and computationally expensive, and the tradeoffs among different model parameters are typically strong (e.g., Porter et al., 2011; Wirth et al., 2017). Simplified modeling using harmonic decomposition or similar techniques (e.g., Shiomi and Park, 2008; Bianchi et al., 2010; Schulte-Pelkum and Mahan, 2014a; Liu et al., 2015; Ford et al., 2016; Olugboji and Park, 2016; Park and Levin, 2016) can identify the presence and primary features of dipping and/or anisotropic interfaces without extensive computation of synthetic seismograms.

In this study we examine data from six long-running broadband seismic stations in the eastern United States (Fig. 1) that overlie lithosphere that has been affected by past mountain-building events, i.e., the Grenville and Appalachian orogenies. The goal of this work is to explore the nature and depth extent of deformation of the crust and mantle lithosphere due to orogenesis using the technique of anisotropic receiver function analysis. We apply a harmonic decomposition modeling approach to our computed receiver functions to identify interfaces within the continental lithosphere with a dipping and/or anisotropic character beneath our selected stations. Three of our stations are just to the east of the Grenville Front, which marks the mapped westward extent of deformation during the Grenville orogeny, and three of our stations are situated within the Appalachians. This particular study is motivated by previous work (Wirth and Long, 2014) that identified evidence for a contrast in anisotropy within the mantle lithosphere at stations in the Granite-Rhyolite Province of eastern North America. The comparison of our new results with these previous results, obtained at stations located to the west of the Grenville Front, allows us to investigate whether and how Grenvillian and Appalachian orogenesis modified the anisotropic structure of the lithosphere.

The eastern margin of North America (Fig. 1), today a passive margin in the interior of the North American plate, has been shaped by two complete supercontinent cycles over the past ∼1.3 b.y. of Earth history (e.g., Thomas, 2004; Cawood and Buchan, 2007; Benoit et al., 2014). This tectonic history includes multiple and protracted episodes of orogenesis, including the Grenville orogenic cycle that culminated in the formation of the supercontinent Rodinia (e.g., Whitmeyer and Karlstrom, 2007; McLelland et al., 2010) and the Appalachian orogenic cycle that culminated in the formation of the supercontinent Pangea (e.g., Hatcher, 2010; Hibbard et al., 2010). Dispersal of each of these supercontinents was accomplished via continental rifting, including a major diachronous rifting episode between ca. 750 and 550 Ma (e.g., Li et al., 2008; Burton and Southworth, 2010) that broke apart Rodinia and shaped the pre-Appalachian eastern margin of Laurentia (e.g., Allen et al., 2010). The later breakup of Pangea was also accomplished via a complex set of rifting processes that was complete by ca. 190 Ma (e.g., Frizon de Lamotte et al., 2015) but began substantially earlier (ca. 250 Ma), and was accompanied by voluminous magmatism associated with the Central Atlantic Magmatic Province, a large igneous province (e.g., McHone, 1996, 2000; Schlische et al., 2003). These multiple tectonic episodes were likely accompanied by substantial deformation of the crust and mantle lithosphere as the edges of the (progressively growing) North American continent were modified.

The term Grenville orogeny is used somewhat inconsistently in the literature (e.g., McLelland et al., 2010), but here we use the term to encompass both the Elzeverian orogeny ca. 1.3–1.2 Ga, which sutured the Elzevir block to the Laurentian margin (e.g., Moore and Thompson, 1980; Whitmeyer and Karlstrom, 2007; Bartholomew and Hatcher, 2010), as well as later (ca. 1.09–0.98 Ga; Whitmeyer and Karlstrom, 2007) continent-continent collision that formed Rodinia (e.g., McLelland et al., 1996; Rivers, 1997) as the Grenville Province was joined to Laurentia. The Grenville deformation front, which represents the westward extent of deformation due to Grenvillian orogenesis (e.g., Culotta et al., 1990; Whitmeyer and Karlstrom, 2007), cuts through eastern Michigan and western Ohio and extends to the south and west through Alabama (Fig. 1). In the context of the work presented here, the Grenville Front separates the seismic stations within the Granite-Rhyolite Province examined in Wirth and Long (2014) and those in this study. Of the stations we examine in this study (Fig. 1), station ACSO is located within the Granite-Rhyolite Province but to the east of the Grenville Front, while station ERPA is located near the suture between the Elzevir block and the rest of the Granite-Rhyolite Province and station BINY is located within the Grenville Province, as defined by Whitmeyer and Karlstrom (2007).

The later Appalachian orogeny, a detailed overview of which was given by Hatcher (2010), involved three distinct phases of orogenesis over a period of several hundred million years, and formed the present-day Appalachian Mountains. The first, the Taconic orogeny, involved the accretion of arc terranes onto the margin of Laurentia between ca. 496 and 428 Ma (Karabinos et al., 1998; Hatcher, 2010), while the later phases (the Acadian-Neoacadian and Alleghanian orogenies) involved superterrane accretion (Carolina, Avalon, Gander, and Meguma superterranes, all of peri-Gondwanan affinity) and continental collision. The Acadian orogeny began ca. 410 Ma and primarily affected the northern Appalachians with transpressive north to south collision (e.g., Hatcher, 2010; Ver Straeten, 2010). The Alleghanian orogeny (e.g., Geiser and Engelder, 1983; Sacks and Secor, 1990; Hatcher, 2010; Bartholomew and Whitaker, 2010), which began ca. 300 Ma and was complete by ca. 250 Ma, culminated in the final assembly of the Pangea supercontinent. The protracted series of orogenic events that make up the Appalachian orogenic cycle produced widespread deformation, volcanism, and metamorphism, and produced the topographic signature of the Appalachian Mountains that is still evident today (Hatcher, 2010). Three of the seismic stations that we use in this study (SSPA, MCWV, and TZTN) are located in the Appalachian Mountains, along the edge of the Laurentian core of the North American continent, in regions that were likely affected by deformation associated with Appalachian orogenesis.

In this paper we use anisotropic receiver function (RF) analysis to identify and characterize contrasts in seismic anisotropy with depth within the continental crust and mantle lithosphere, following our previous work (Wirth and Long, 2014; Ford et al., 2016). This technique involves the computation of RFs, i.e., a time series that reflects converted seismic phases due to structure beneath a seismic station (e.g., Langston, 1979), for both the radial (oriented in a direction that points from receiver to source) and transverse (oriented orthogonal to the radial) horizontal components of a seismogram. For the case of horizontal interfaces and purely isotropic structure, all of the P to S converted energy would be expected to arrive on the radial component, as the P and SV wavefields are coupled and (in this case) propagate independently of the SH wavefield. When dipping interfaces or anisotropic structure are present, some of the P to S converted energy will arrive on the transverse component due to P to SH conversions (e.g., Levin and Park, 1997, 1998). For the case of a dipping isotropic interface, these converted waves will manifest as an arrival on the transverse component RF trace that changes polarity twice across the full backazimuthal range (i.e., a two-lobed polarity flip). For the case of a contrast in anisotropy (with a hexagonal symmetry and a horizontal symmetry axis), these conversions will manifest as an arrival on the transverse component RF trace that changes polarity four times across the full backazimuthal range (i.e., a four-lobed polarity flip). For the case of a contrast in anisotropy with a plunging axis of symmetry, a mix of two-lobe and four-lobe polarity changes on the transverse component RFs is predicted. One way to distinguish between a dipping (isotropic) interface and a contrast in anisotropy with a plunging symmetry axis, at least in theory, is to examine the transverse components in the time range associated with the direct P arrival. For a dipping interface, a change in polarity on the transverse components at a time associated with the P to S conversion at an interface will be accompanied by an arrival with the opposite polarity at time t = 0; for a contrast in anisotropy at depth, this zero-time arrival is not present. One caveat, however, is that this analysis is complicated when there are multiple interfaces at depth, because multiple zero-time arrivals may interfere. Predicted RF patterns for a suite of simple models that illustrate these ideas were shown and discussed in detail in Schulte-Pelkum and Mahan (2014b) and Ford et al. (2016).

We selected six long-running broadband seismic stations in the eastern United States for analysis in this study, all part of the permanent U.S. National Seismic Network or Global Seismographic Network (Fig. 1). In order to achieve good azimuthal coverage, we focus on a small number of high-quality stations with long run times. Of these, stations ACSO, ERPA, and BINY are located in regions that were likely affected by Grenville orogenesis, as described here; we collectively refer to these as the Grenville stations in this study. Stations MCWV, SSPA, and TZTN are located in regions of present-day mountainous topography that were likely affected by Appalachian orogenesis (Mountain stations). We examined at least 10 yr of data at each station, selecting events of magnitude 5.8 and greater at epicentral distances between 30° and 100° (Fig. 2) for analysis.

Our preprocessing methodology is identical to that of Ford et al. (2016); each seismogram trace was cut to identical length, rotated into vertical, radial, and transverse components, and bandpass filtered between 0.02 and 2 Hz for P wave picking. We visually inspected each trace for a clear P wave arrival, discarding those traces that did not display an unambiguous P wave, and manually picked the P arrival using the Seismic Analysis Code (SAC) software (http://ds.iris.edu/ds/nodes/dmc/software/downloads/sac/). Before computing RFs, we rotated the components into the LQT (theoretical directions of P, SV, and SH motions) reference frame (using a near-surface P wave velocity of 3.5 km/s) to account for the fact that the incoming P arrivals are not perfectly vertical, again following Ford et al. (2016). We note that while all RFs were computed in the LQT coordinate system, for simplicity we refer to radial and transverse RFs to follow common terminology conventions (e.g., Levin and Park, 1997, 1998). We emphasize, however, that the LQT and RTZ (radial, transverse, vertical) coordinate systems are not strictly equivalent, and that our rotation into the LQT coordinate system is imperfect, because we used a single representative near-surface velocity to calculate the rotations. Finally, before RF computation, we applied a final bandpass filter to each waveform, using a highpass cutoff of 0.02 Hz and a variable lowpass cutoff of 0.5 or 1 Hz.

To compute radial and transverse component RFs, we used the frequency domain multitaper correlation estimator of Park and Levin (2000), which is a commonly applied RF computation technique (recent examples can be found in Wirth and Long, 2012, 2014; McCormack et al., 2013; Yuan and Levin, 2014; Liu et al., 2015; Levin et al., 2016; Olugboji and Park, 2016; Ford et al., 2016). Following Ford et al. (2016), we computed individual RFs and then stacked after corrections for variations in slowness of the direct P wave arrival, which depends on epicentral distance. We computed single-station radial component stacks over all backazimuths, as well as stacks of both radial and transverse component RFs as a function of both epicentral distance (which allows us to evaluate the effects of multiply reflected waves on our results) and backazimuth (which allows us to identify dipping and/or anisotropic structure). Our stacking and plotting conventions are shown in Figure 3, which illustrates a set of RF stacks and plots at station SSPA as an example. In addition to the RF stacks shown in Figures 3A and 3B, we also visualize the backazimuthal variability in transverse component RF energy within specific time ranges (corresponding to specific depth ranges within the lithosphere) using rose diagrams, following the convention in Wirth and Long (2014) and as illustrated in Figure 3D. This plotting convention allows for visual inspection of the azimuthal variability in polarity (represented with blue and red colors) at specific time (depth) ranges. The character of any variability (two lobed, four lobed, or a mixture) sheds light on the type of interface (dipping, anisotropic, or both). Distinguishing between a dipping interface and a plunging axis of anisotropy requires an inspection of the transverse component RFs at the arrival time of the direct P wave, as noted here, although in the presence of multiple interfaces this can be difficult.

We implement a harmonic decomposition technique (Shiomi and Park, 2008; Bianchi et al., 2010; Park and Levin, 2016) that allows for the straightforward identification of time ranges within the RF traces that are likely associated with conversions at interfaces of different character. Specifically, harmonic decomposition identifies coherent energy on the RF traces that correspond to arrivals that are constant across backazimuth, that vary as a function of cos(θ) or sin(θ) (where θ is backazimuth), that vary as a function of cos(2θ) or sin(2θ), or a combination of these. Harmonic decomposition modeling thus serves as a useful complement to the simple visual inspection of backazimuthal RF gathers, and has been applied by several recent studies (e.g., Liu et al., 2015; Olugboji and Park, 2016; Ford et al., 2016). A detailed description of our implementation of harmonic decomposition can be found in Ford et al. (2016) and a recent technical overview was given by Park and Levin (2016). Briefly, for each time window the stacked RFs (both radial and transverse components) are modeled as a linear combination of sin(kθ) and cos(kθ) terms (after application of a phase shift that depends on k), where k = 0, 1, or 2. The k = 0 term, with no backazimuthal dependence, suggests an isotropic velocity change across a flat interface; the k = 1 terms imply a dipping interface and/or a dipping axis of anisotropic symmetry; the k = 2 terms correspond to contrast in anisotropy. The harmonic stacking technique is illustrated for example station SSPA in Figure 3, which shows the modeled harmonic decomposition as a function of time (along with an error estimate derived from bootstrap resampling; see Ford et al., 2016 for details). Figure 3 also shows the unmodeled portion of the signal, which cannot be represented as a linear sum of sin(kθ) and cos(kθ) terms. Here the unmodeled components for each value of k represent the portion of the signal that has the opposite phase shift as that of the harmonic expansion model (Shiomi and Park, 2008).

For each of the six stations examined in this study, we computed stacked radial RFs over all backazimuths for two different lowpass frequency cutoffs (Fig. 4) as well as backazimuthal and epicentral distance gathers (GSA Data Repository Figs. S1–S12¹). Here we focus our presentation of the results and discussion of our interpretation on the single-station stacked radial RFs (Fig. 4) as well as the harmonic decomposition results for each station, grouped by region (Grenville stations in Fig. 5 and Mountain stations in Fig. 6). We use the harmonic decomposition modeling to guide our interpretation of major interfaces and their anisotropic or isotropic character beneath each station, as shown in Figures 7 and 8 and discussed in the following. While we do not discuss the epicentral distance gathers in detail, for each of the interfaces that we interpret in this paper, we have checked the epicentral distance gathers (odd numbered Figs. S1–S12) to ensure that these arrivals do not exhibit moveout that suggests that they are contaminated by multiply reflected waves from the crust. We note that while our use of the LQT coordinate system should theoretically remove the direct P arrival from our RF stacks, some energy at or near t = 0 is often visible (Fig. 4; Figs. S1–S12). This reflects minor inaccuracies in the velocity model used to predict the geometry of the incident P wave; an offset from t = 0 in this arrival (e.g., stations BINY and MCWV in Fig. 4) may reflect the presence of thick sedimentary cover beneath the stations or (more likely) shallow intracrustal interfaces.

The stacked radial RFs shown in Figure 4, which reflect events arriving across the full backazimuthal range and thus obscure any dipping or anisotropic structure, yield evidence for complex structure with a number of interfaces beneath each of the stations examined in this study. At our lower frequency range (lowpass cutoff at 0.5 Hz; Fig. 4A), we observe at all stations a strong, clear positive pulse (light blue lines in Fig. 4) ∼5–5.8 s after the P arrival, which we interpret as corresponding to the positive (isotropic) velocity change across the Moho. This time interval corresponds to a crustal thickness of ∼38–48 km for reasonable velocity models for continental crust, and is generally consistent with previously published estimates of crustal thickness in our study region (e.g., Abt et al., 2010). In addition to the Moho, we also see evidence for (typically two or more) negative pulses (corresponding to a decrease in velocity with depth) at later time ranges (∼8 s and greater; pink lines in Fig. 4), which may correspond either to mid-lithospheric discontinuities (MLDs; Abt et al., 2010; Ford et al., 2016) or, perhaps, to the base of the lithosphere. When higher frequency energy (as high as 1 Hz) is included, the radial RF stacks (Fig. 4B) have similar character but more detail, often hinting at intracrustal layering (pulses arriving before the Moho conversion).

Figures 5 and 6 show the results of the harmonic decomposition analysis for all stations, computed for our low-frequency data (cutoff frequency of 0.5 Hz, the same as the RFs shown in Figs. S1–S12). Focusing first on the Grenville stations (ACSO, BINY, and ERPA; Fig. 5), we find that the harmonic stacking approach is generally able to model most of the features of the RFs (note relatively small amplitudes on the unmodeled stacks) and yields evidence for a number of interfaces with different character. Beneath station ACSO, the Moho is expressed as a large-amplitude positive pulse on the constant term (at t = 5.8 s relative to direct P); in the same time window, there is some energy (above the error estimate, shown on each trace) on the k = 1 and k = 2 terms, although it is weaker. Still, the presence of energy on the k = 1 and k = 2 terms arriving at the same time as the constant Moho arrival suggests a contrast in anisotropy across the Moho, perhaps with a dipping symmetry axis, and/or a dipping Moho interface. We also see evidence on the constant term stack for an isotropic velocity drop at a depth internal to the mantle lithosphere (t = 9.5 s, corresponding to ∼90 km based on the ak135 Earth model of Kennett et al., 1995); the k = 1 and k = 2 terms in this time range suggest a contrast in anisotropy and/or a dip to this interface as well. This inferred MLD phase does not show a clear moveout in the epicentral distance stacks (Fig. S2), which argues against significant contamination from crustal multiples; a similar argument can be made for each of the mantle interfaces we discuss herein, with one exception noted, as follows.

Returning to the harmonic decomposition analysis in Figure 5 and focusing purely on the k = 1 and k = 2 components, we can identify peaks corresponding to interfaces at mantle depths beneath ACSO that suggest a contrast in anisotropy (perhaps with a dipping symmetry axis), but do not necessarily correspond to high amplitudes on the constant term (meaning an isotropic velocity change is not required). In the time range associated with the direct P arrival, the transverse component RFs show little or no energy at most backazimuths (Fig. S1); we do not observe any convincing zero-time arrivals that might result from a dipping interface, so contrasts in anisotropy with a plunging symmetry axis may be more likely. We emphasize, however, that the presence of multiple interfaces with large k = 1 terms in the harmonic expansion implies that an examination of the zero-time transverse component arrivals to distinguish between dipping interfaces and plunging anisotropy is not straightforward.

Beneath station BINY (Fig. 5), the general character of the inferred interfaces is similar to what we infer beneath ACSO, but the details are different. In the time window corresponding to the Moho pulse (t = 5.2 s) on the constant term, there is clear energy on both the sin(2θ) and cos(2θ) terms, again suggesting a contrast in anisotropy across the Moho. At later times, corresponding to depths internal to the mantle lithosphere, we see a number of interfaces evident on both the constant and the higher order (k = 1, 2) terms, but the exact timing of these arrivals, and the anisotropic and/or dipping geometry suggested by relative strength of the cos(θ), sin(θ), cos(2θ), and sin(2θ) stacks, are different than for station ACSO. Similarly, beneath ERPA, we see the same general pattern, with both isotropic and anisotropic/dipping interfaces suggested by the harmonic decomposition, but with different specific characteristics. It is striking that ERPA exhibits an interface at a time of ∼7.8 s (relative to direct P; ∼75–80 km depth based on ak135) that shows both a strong (negative) constant component and a strong cos(θ) component. This strong k = 1 term may correspond to either a dipping MLD or a contrast in anisotropy across the MLD that involves a plunging axis of anisotropy; the very weak signal at zero delay time (relative to P) on the transverse component RFs for this station (Fig. S5) suggests that the latter possibility is more likely. There is another strong negative arrival on the constant stack (t = 12 s; ∼120 km depth) that also shows large amplitudes for cos(θ), sin(θ), and sin(2θ), again suggesting a contrast in anisotropy with a plunging symmetry axis.

Our harmonic stacking results for the Mountain stations (Fig. 6) reveal a series of complicated interfaces, many with dipping and/or anisotropic character, beneath each station. At each of the three stations, we observe a clear Moho arrival on the constant term, and at two of the three stations (SSPA and TZTN), the higher order stacks also exhibit clear peaks in the same time window, suggesting an anisotropic and/or dipping character to this interface. (At MCWV, there is an anisotropic intracrustal interface whose conversions arrive ∼1 s ahead of the Moho Ps phase, but it appears to be distinct from the Moho.) The constant terms at each of the three stations at time windows corresponding to mantle depths suggest one or more interfaces with a velocity drop, implying the likely presence of multiple MLDs; typically, these are also associated with an anisotropic or dipping character, as revealed by the higher order terms. For the most part, the epicentral distance gathers argue against contamination of crustal multiples in our interpretation of mantle interfaces beneath these stations, although we observe some moveout for the negative arrival at times near t = 9 s for station TZTN (Fig. S12); its interpretation should therefore be treated with some caution. An examination of the transverse RFs gathers at t = 0 for the Mountain stations (Figs. S7, S9, and S11) reveals modest amounts of energy, with complex backazimuthal patterns. For each of the three stations, there is nonzero energy at t = 0, suggesting the presence of dipping interfaces at depth, but the azimuthal variability is not clear enough to identify individual two-lobed polarity changes and associate them with specific interfaces.

As with the Grenville stations, we are also able to pick out features on the k = 1 and k = 2 stacks at individual Mountain stations that do not obviously correspond to features on the constant stacks, suggesting that there is some anisotropic layering with the lithosphere that is not accompanied by an isotropic change in velocity. As with the Grenville stations, the collection of harmonic stacks for the Mountain stations, taken as a whole, suggests pronounced differences in the details of both isotropic and anisotropic layering among different stations. Visual inspection of the three stations shown in Figure 6 reveals major differences in the timing and character (i.e., the geometry of polarity reversals as revealed by the harmonic stacking technique) of individual interfaces among different stations.

A key advantage of the harmonic stacking technique is that it allows us to identify conversions associated with interfaces that exhibit a dipping and/or anisotropic character, a task that is somewhat difficult, relying only on visual inspection. From the harmonic stacks shown in Figures 5 and 6, we have picked a series of time ranges (indicated by gray lines in Figs. 5 and 6, with dashed lines for interfaces within the crust and solid lines for interfaces within the mantle) that indicate dipping and/or anisotropic interfaces beneath each station. Our methodology for picking these interfaces follows Ford et al. (2016) and involves a summation of each of the four nonconstant terms at each delay time in order to determine where coherent peaks in energy occur. Another way of visualizing the character of these interfaces, which lends itself well to comparison among different stations, is to plot rose diagrams (as illustrated in Fig. 3) for each of these time ranges, showing the variability in the transverse component RF polarities and amplitudes as a function of backazimuth. In Figure 7, we show rose diagrams interfaces at crustal depths for each station (either within the crust or colocated with the isotropic Moho), grouped by tectonic setting (Grenville versus Mountain stations). In Figure 8, we show similar rose diagrams for interfaces within the mantle lithosphere; this figure also includes a comparison with data from our previous study using stations within the Granite-Rhyolite Province (Wirth and Long, 2014). We have also used the harmonic decomposition results to calculate azimuthal phase information for individual interfaces (for both the k = 1 and k = 2 terms) that corresponds to the azimuthal location of polarity changes in the rose plots in Figure 7 and 8 and that contains information about the possible dip or plunge directions for a dipping interface or plunging symmetry axis (for the k = 1 components), and information about contrasts in the orientation of horizontal azimuthal anisotropy (for the k = 2 components). These orientations are shown in map view in Figures 9–11, and discussed in detail in the following.

Taken together, the collection of rose diagrams shown in Figures 7 and 8 illustrates that our stations overlie highly complex lithosphere, with a number of interfaces both within the crust and within the mantle lithosphere that dip and/or feature a contrast in anisotropy, perhaps with dipping symmetry axes. This is true both at crustal depths (Fig. 7) and within the mantle lithosphere (Fig. 8), and is similar to our previous findings within Archean cratonic regions of North America (the Wyoming and Superior cratons; Ford et al., 2016). In Ford et al. (2016), we identified a surprising amount of both vertical and lateral heterogeneity in lithospheric structure, suggesting that the complex tectonic processes that produce very heterogeneous structure in the surface geology also produce similarly complex structure within the deep lithosphere. It is strikingly evident from Figures 7 and 8 that there is significant lateral heterogeneity within the crust and lithosphere in our eastern United States study region, even when we compare within groups of stations. Stations that overlie lithosphere likely affected by Grenville orogenesis show little similarity in the timing, character, or geometry of dipping/anisotropic structure across stations; the same is true for the Mountain stations that overlie lithosphere likely affected by Appalachian orogenesis.

A more detailed examination of the crustal interfaces illustrated in Figure 7 yields some clues as to their geometry and possible origin. Beneath ACSO, we identify a dipping interface within the crust, with conversions arriving ∼3.9 s after the direct P arrival. The energy on the transverse component RFs at t = 0 may result from this dipping interface, although the azimuthal variability in its polarity is less clear-cut than the azimuthal variability in the Ps arrival at ∼3.9 s. Lacking a detailed crustal velocity model, we cannot accurately estimate the depth of this interface, but its timing suggests mid-crustal depths. The arrival is modeled well with the k = 1 terms in the harmonic expansion (Fig. 5), so anisotropy need not be invoked. Also beneath ACSO, we see evidence for an anisotropic contrast across the Moho; in the time range (∼5.8 s) associated with the constant (isotropic) Moho arrival, the transverse component exhibits significant azimuthal variability (Fig. 7), with a nonzero sin(2θ) term. Visually, the transverse components associated with the Moho arrival look similar to those for the mid-crustal arrival (upper left panels in Fig. 7), but an additional anisotropic component is suggested by the harmonic stacks. Beneath both BINY and ERPA, and similar to ACSO, there is variability on the transverse component in the time range associated with Moho conversions (Fig. 7) that suggests the presence of a contrast in anisotropy (see nonzero k = 2 terms in this time range in the harmonic stacks in Fig. 5) as well as a dipping component. We see evidence for a contrast in anisotropic properties between the lower crust and the uppermost mantle at all three Grenville stations, suggesting anisotropy in one or both layers.

Our inferences on crustal structure beneath the Mountain stations (Fig. 7) include our identification beneath station MCWV of two intracrustal interfaces (at ∼2.2 s and ∼4.3 s; both times are before the Moho arrival), each of which exhibit significant energy in both the k = 1 and k = 2 stacks (Fig. 6) and thus require some contribution from anisotropic structure. This observation strongly suggests at the presence of multiple anisotropic layers within the middle to lower crust beneath MCWV, likely with dipping anisotropic symmetry axes (and/or with a dip to the interfaces). Similarly, beneath SSPA, we infer the presence of an interface in the lowermost crust (arrival time ∼4.3 s, just before the Moho arrival but with k = 1 and k = 2 peaks that clearly arrive before the k = 0 Moho peak; Fig. 6) that also requires a contrast in anisotropy. The rose diagram for this interface (Fig. 7) clearly shows four polarity changes across the full backazimuthal range, although this pattern is modulated by k = 1 terms as well (Fig. 6). Beneath TZTN, the only identifiable crustal anisotropic contrast is clearly associated with the isotropic Moho, suggesting a contrast between the lower crust and the uppermost mantle, similar to what we observe beneath all three of the Grenville stations.

What are the implications of our findings about anisotropic layering at crustal depths for our understanding of past tectonic deformation? First, our inference of a contrast in anisotropic properties across the Moho at four of our stations (ACSO, BINY, ERPA, and TZTN) suggests significant, “frozen-in” anisotropy in the relatively shallow portions of the lithosphere (lower crust and/or uppermost mantle) due to deformation associated with past tectonic events. The character of the transverse RFs in the time range associated with the Moho does not strictly require anisotropy in the lowermost crust beneath these four stations, but it does imply anisotropy in the deepest crust, the shallowest mantle, or both. In either case, our observations suggest past deformation at depths near the Moho, with enough strain to develop coherent CPO of anisotropic minerals. Second, our observations argue strongly for multiple anisotropic layers within the crust, at middle to lower crustal depths, beneath two of our stations (MCWV and SSPA), both located within the Appalachian Mountains. We interpret this finding as suggesting that (1) the deep crust has undergone significant deformation in the past and has developed anisotropy via the CPO of crustal minerals, most likely associated with Appalachian orogenesis, and (2) the mineralogy, rheology, and/or deformation geometry varied as a function of depth within the crustal column, resulting in multiple layers of crust with contrasting anisotropic geometries. Our inference of layered crustal anisotropy in the ancient Appalachian orogeny is similar to inferences from modern orogens such as Tibet (Liu et al., 2015) and Taiwan (Huang et al., 2015), and provides additional support for the hypothesis that crustal deformation in orogenic systems is complex, varies with depth, and extends to the deep crust.

Regarding the lithospheric mantle, we examine the character of transverse component RFs at time ranges associated with anisotropic and/or dipping contrasts within the lithospheric mantle, shown in Figure 8. Of particular interest are the possible relationships between contrasts in anisotropy and the MLDs inferred from the single-station radial RF stacks shown in Figure 3. The term mid-lithospheric discontinuity, coined by Abt et al. (2010), refers to a sharp decrease in seismic velocity at a depth internal to the continental mantle lithosphere. MLDs have been documented in a number of continental regions (e.g., Abt et al., 2010; Ford et al., 2010; Foster et al., 2014; Hopper et al., 2014; Hopper and Fischer, 2015), although their origin remains debated (e.g., Selway et al., 2015; Karato et al., 2015; Rader et al., 2015). Several workers have suggested a link between anisotropic layering in the continental lithosphere and the observation of MLDs (e.g., Yuan and Romanowicz, 2010; Sodoudi et al., 2013; Wirth and Long, 2014). In particular, our previous work using data from the Granite-Rhyolite Province of North America (Wirth and Long, 2014) found evidence for a clear contrast in anisotropy in the same depth range (∼90 km) as the MLD, although our observations also required a decrease in isotropic seismic velocity with depth and could not be explained solely by anisotropic layering within the lithosphere.

Figure 8 demonstrates the presence of complex dipping and/or anisotropic layering within the lithospheric mantle beneath the eastern United States, with 2–3 clear interfaces that correspond clearly to transverse RF polarity changes with backazimuth, arriving ∼7–12 s after the P arrival, identified at each station. This finding is similar to our previous work using data from the Superior and Wyoming cratons in North America (Ford et al., 2016), where we typically identified multiple layers of anisotropy within the mantle lithosphere. Comparing the rose diagrams in Figure 8 to the single-station radial component RF stacks in Figure 3, we can evaluate whether there is a correspondence (in converted wave arrival timing and thus in interface depth) between isotropic MLDs (corresponding to red or negative pulses in Fig. 3) and contrasts in anisotropy and/or dipping interfaces (as inferred from the harmonic decomposition and illustrated with rose diagrams in Fig. 8). Beneath stations ACSO, ERPA, MCWV, and SSPA, we are able to identify interfaces that apparently exhibit both an isotropic drop in velocity with depth as well as an anisotropic or dipping character. For example, beneath ACSO there is a clear discontinuity that produces a negative pulse on the stacked radial RFs at 9.5 s; the corresponding rose diagram in Figure 8 exhibits clear backazimuthal variations that include an anisotropic component. Similarly, beneath SSPA there are two clear negative pulses on the stacked radial RFs (at t = 6.9 s and 10.3 s) that also require a contrast in anisotropy, likely with a dipping component (nonzero amplitudes in the k = 1 and k = 2 stacks in Fig. 6). However, again similar to the findings of Ford et al. (2016), we do not observe a simple, one to one correspondence between the isotropic MLD arrivals and the anisotropic interfaces. Rather, we observe complex layering within the lithosphere with both anisotropic and isotropic interfaces, and only sporadic correlations between isotropic and anisotropic structure.

Constraints on the geometry of the RF polarity changes due to anisotropic and/or dipping interfaces beneath our stations are visualized in map view in Figures 9–11, which show the phase orientation information for the k = 1 and k = 2 terms in the harmonic expansion for a series of interfaces at different depths. These directions are derived from measurements of the relative amplitudes of the sine and cosine terms, and denote the possible orientations that may correspond to the dip direction (for a dipping interface) or the plunge direction (for a plunging symmetry axis) for the k = 1 terms, and contain information about contrasts across interfaces in symmetry axis orientation (for azimuthal anisotropy) for the k = 2 terms. We caution that these orientations do not uniquely constrain the geometry of anisotropy or of dipping interfaces; in the absence of detailed forward modeling, assumptions must be made about the properties of individual layers in order to estimate the anisotropic geometry. Nevertheless, the phase data contain some information about possible dip orientations and about the geometrical relationships between azimuthally anisotropic symmetry directions across interfaces (e.g., Shiomi and Park, 2008; Schulte-Pelkum and Mahan, 2014a, 2014b; Olugboji and Park, 2016), and regional variations in the possible orientations can be easily evaluated in map view.

Figure 9 shows the phase information for interfaces within the crust and across the Moho for those stations where the Moho interface includes a dipping and/or anisotropic component. This map reveals significant lateral variability in the character of crustal interfaces across the study region, and little in the way of straightforward relationships between the phase orientations and geologic and tectonic indicators. In Figure 10, we show phase orientations for interfaces internal to the lithosphere that do not correspond to strong isotropic contrasts in velocity; rather, these are shallow intralithospheric contrasts in anisotropy, typically with a dipping axis of symmetry. As with the crustal interfaces in Figure 9, we observe a great deal of regional variability. In Figure 11, we show orientations for interfaces in the mid-lithosphere that exhibit both an isotropic drop in velocity with depth, thus corresponding to isotropic MLDs, as well as a dipping and/or anisotropic character. There is some consistency between the k = 1 phase orientations at stations ERPA, MCWV, and SSPA, two of which (ERPA and SSPA) exhibit multiple anisotropic MLD interfaces. All of these have k = 1 orientations directed approximately northwest-southeast, with more variability in the k = 2 phase information.

A direct comparison between our results and those in Wirth and Long (2014) is instructive, particularly for stations ACSO and ERPA, which are located within the Granite-Rhyolite Province (ACSO) or in the adjacent and related Elzevir block (ERPA). Figure 8 shows rose diagrams for arrivals in the time range between ∼7.8 and 10.5 s (corresponding to depths between ∼80–105 km) at Granite-Rhyolite stations; at each of these stations, evidence for layered anisotropy within the mantle lithosphere was identified (Wirth and Long, 2014) with conversions from one or more anisotropic interfaces arriving in the same time range as conversions from the isotropic MLD. Furthermore, at three of the Granite-Rhyolite stations (BWI, BLO, and CCM), the transverse component RFs associated with the main MLD interface exhibited similar behavior, with four clear changes in polarity across the full backazimuthal range. A forward model of the data at station WCI (Wirth and Long, 2014) invokes three lithospheric layers of anisotropy; across the likely main MLD interface, the upper model layer has an approximately north-south anisotropic fast direction and the lower layer is nearly east-west. This model geometry predicts polarity changes similar to those observed at WCI (Fig. 8, lower left), with negative transverse component arrivals in the northeast and southwest quadrants and positive arrivals in the northwest and southeast quadrants. A comparison of this specific geometric pattern with the transverse RF behavior at stations ERPA and ACSO in similar time ranges (Fig. 8, upper right) demonstrates that the backazimuthal patterns are distinctly different beneath the Grenville stations in our study. The fact that station ACSO in particular, which is also located within the Granite-Rhyolite Province (similar to nearby stations BLO and WCI from Wirth and Long, 2014) but to the east of the Grenville Front, does not exhibit a similar geometric pattern implies that the lithosphere beneath ACSO was modified via subsequent deformation. We suggest that this lithosphere was deformed during Grenville orogenesis, thus modifying the preexisting lithospheric structure that is still present elsewhere in the mantle of the Granite-Rhyolite Province.

How do our observations compare to other previous inferences on the anisotropy of the crust and mantle lithosphere beneath eastern North America? In general terms, our results are consistent with previous work that has suggested a significant contribution to seismic observations from anisotropy in the mantle lithosphere, based on surface waves (e.g., Deschamps et al., 2008), SKS splitting patterns (e.g., Long et al., 2016), or on combinations of different types of data (e.g., Yuan and Romanowicz, 2010; Yuan and Levin, 2014). Recent work on Pn velocities and anisotropy beneath the continental United States (Buehler and Shearer, 2017) also suggests significant anisotropy in the uppermost mantle beneath our study region; Buehler and Shearer (2017) further propose that there must be multiple layers of anisotropy within the upper mantle, potentially consistent with our observations of anisotropic layering within the mantle lithosphere. Our finding of significant anisotropic layering with the middle to lower crust beneath at least some of our stations is also generally consistent with results of Schulte-Pelkum and Mahan (2014a), who suggested crustal anisotropy (often with a plunging symmetry axis) beneath a number of USArray Transportable Array (TA) stations in the eastern United States, particularly in the Appalachians (see their Fig. 8). It is interesting that our finding that the layered anisotropic structure of the mantle lithosphere differs across the Grenville deformation front contrasts with the view provided by shear wave splitting observations, which reflect an integrated signal from the entire upper mantle. For example, Sénéchal et al. (1996) observed similar SKS splitting on either side of the Grenville Front in Canada, while in Long et al. (2016) it was noted that there is no obvious correlation between SKS splitting patterns measured at TA stations and the geometry of the Grenville Front in the eastern and central United States.

A major limitation of the anisotropic RF technique is that while we can confidently infer the presence (and roughly estimate the depths) of anisotropic interfaces within the lithosphere, as well as some orientation information from the k = 1 and k = 2 phase values from the harmonic stacking, it is difficult to infer the geometry of anisotropy within each layer without detailed forward modeling. Such forward modeling is computationally intensive and nonunique, with a large number of unknown parameters, so while it can identify plausible models, it often cannot uniquely determine the anisotropy geometry in every layer. These difficulties can be ameliorated by the use of model space search approaches (e.g., Porter et al., 2011; Wirth et al., 2017) or Bayesian inversion schemes (e.g., Bodin et al., 2016), but there are still strong tradeoffs among parameters. In lieu of detailed forward modeling, the harmonic decomposition modeling applied here (or approaches similar to it; e.g., Schulte-Pelkum and Mahan, 2014a) can identify interfaces and provide some primary information about their geometry by examining the relative amplitudes and phases of the k = 0, 1, 2 terms, and provides a natural way of comparing structure among different stations. Furthermore, the RF observations and the harmonic decomposition models presented in this paper can serve as a starting point for future detailed forward modeling, ideally in combination with other data that constrain seismic anisotropy (such as SKS splitting observations, surface wave dispersion, and/or Pn traveltimes). One example of this type of modeling was discussed by Yuan and Levin (2014), who combined different types of anisotropy observations, including backazimuthal RF gathers, into a forward model of multilayered anisotropy beneath three stations in the eastern United States (including SSPA, one of the stations examined in this study).

Despite the limitations of the RF technique, the observations presented in this paper allow us to draw some straightforward inferences about lithospheric deformation beneath the eastern United States and its likely causes. Our identification of extensive anisotropy in the deep crust and mantle lithosphere beneath each of the stations examined in this study implies extensive lithospheric deformation, with enough strain to generate significant CPO. While it is not possible to conclusively identify the deformation events and their timing and geometry, we suggest that the last major orogenic cycle to affect each region is the most plausible tectonic event to have caused widespread lithospheric deformation (e.g., Meissner et al., 2002). Therefore, we propose that the anisotropic structure of the mantle lithosphere beneath ACSO, ERPA, and BINY was shaped by deformation associated with the Grenville orogenic cycle, while Appalachian orogenesis caused lithospheric deformation beneath SSPA, MCWV, and TZTN. Our observations strongly suggest anisotropic layering within the crust beneath two of the three Appalachian stations examined, with an anisotropic geometry that varies with depth. This finding is similar to the documentation of crustal anisotropy in modern mountain belts (e.g., Liu et al., 2015; Huang et al., 2015), and suggests complex crustal deformation accompanying orogenesis. Our main finding of extensive but layered deformation within the lithospheric mantle (at depths to ∼100 km) presents an interesting challenge to the classic idea of vertically coherent deformation (Silver, 1996). Although our observations suggest that the geometry of anisotropy (and thus deformation) changes with depth, they also suggest significant deformation throughout much of the lithosphere associated with orogenesis. This is consistent with the suggestion by Silver (1996) that the lithosphere participates in deformation, leading to seismic anisotropy, but the finding of changes in anisotropic geometry with depth suggests that the geometry of deformation and/or of olivine fabric development is not coherent over the entire lithospheric mantle.

Our identification of lithospheric anisotropy beneath eastern North America is consistent with observations in other continental regions, including those that have undergone (or are currently undergoing) orogenesis (e.g., Meissner et al., 2002; Wüstefeld et al., 2010). Our main conclusions are also consistent with the many previous suggestions that SKS splitting in and around the Appalachians, which is typically parallel to the mountain belt (at least in the central and southern portions), reflects a significant contribution from lithospheric anisotropy (e.g., Barruol et al., 1997; Wagner et al., 2012; Long et al., 2016). We also note that our inference of multiple layers of anisotropy within the mantle lithosphere is generally consistent with the surface wave model of Deschamps et al. (2008) for a region just to the south and west of our study area. Although there is no substantial geographic overlap between our study region and that of Deschamps et al. (2008), we discussed the similarities and differences between the model of Deschamps et al. (2008) and models derived from RFs in our previous work (Wirth and Long, 2014).

We examined data from six long-running broadband seismic stations located within regions of the eastern United States that were affected by the Grenville and Appalachian orogenic cycles. We computed radial and transverse RFs and applied a harmonic stacking method that allows us to identify interfaces within the crust and lithospheric mantle that have a dipping and/or anisotropic character, in addition to interfaces that can be explained in terms of an isotropic velocity contrast. We find evidence for layered anisotropy, often with a likely dipping axis of symmetry, beneath all stations examined in this study. There is often a clear anisotropic contrast associated with the Moho interface, requiring anisotropy in the lowermost crust, the uppermost mantle, or both. Similar to recent findings within Archean regions of North America, we find evidence for complex anisotropic layering with the mantle lithosphere. At two of the stations in the Appalachian Mountains, we observe clear evidence for anisotropic contrasts within the crust, suggesting middle to lower crustal deformation associated with Appalachian orogenesis. A detailed comparison between our observations and previous anisotropic RF analysis at stations located in the Granite-Rhyolite Province but to the west of the Grenville Front suggests a different geometry of anisotropy in the mantle lithosphere. This, in turn, suggests that the lithospheric mantle just to the east of the Grenville Front was deformed during the Grenville orogenic cycle. Our observations can provide a starting point for the future testing of detailed deformation scenarios associated with past tectonic events, ideally in combination with other types of data that constrain seismic anisotropy.

Seismic waveform data from the U.S. National Seismic Network and the Global Seismographic Network were accessed via the Data Management Center of the Incorporated Research Institutions for Seismology (IRIS). Some figures were prepared using Generic Mapping Tools (Wessel and Smith, 1999). This work was carried out in part through the IRIS Summer Internship Program and was funded by the EarthScope program of the National Science Foundation via grant EAR-1358325. We thank Juan Aragon, Maggie Benoit, Vadim Levin, and Jeffrey Park for helpful discussions and two anonymous reviewers for thoughtful and constructive comments that improved the paper.