LUO ET AL. 10.1029/2023GC011118 2 of 23 of seismic anisotropy can indicate ongoing processes, such as present-day mantle flow in the asthenosphere (e.g., Gung et al., 2003;Lopes et al., 2020), as well as relict structures in the lithosphere from past tectonic events (e.g., Hopper et al., 2017;Long et al., 2017).

The eastern margin of North America has undergone two complete supercontinental cycles, including the Proterozoic Grenville orogenic cycle and the Paleozoic Appalachian orogenic cycle. The Grenville orogenic cycle involved the Mesoproterozoic formation of the supercontinent Rodinia and subsequent Neoproterozoic continental rifting (Bogdanova et al., 2009;Rivers, 1997). The Appalachian orogenic cycle started with diachronous accretion of various Gondwana-derived terranes to the Laurentian margin during the Ordovician Taconic orogeny, the Late Ordovician-Early Silurian Salinic orogeny, the Latest Silurian-Devonian Acadian orogeny, and the Late Devonian Neoacadian orogeny (e.g., Hatcher, 2010;Hibbard et al., 2007Hibbard et al., , 2010;;van Staal et al., 2009). Appalachian orogenesis culminated in the continent-continent collision of composite Laurentia with Gondwana and the formation of the supercontinent Pangea during the Late Mississippian to Permian Alleghanian orogeny (Hatcher, 2002), followed by Mesozoic continental rifting and the breakup of Pangea (Withjack et al., 2012). Multiple episodes of subduction, terrane accretion, continental collision, and rifting involved in these two supercontinent cycles greatly modified both the surface geology and the lithospheric structure of eastern North America. These tectonic processes, while being transient events in the context of Earth's vast geologic timeline, have likely deformed the lithosphere and potentially left a long-lasting anisotropic structure that can be investigated today by seismic techniques.

Shear wave splitting analysis, based on measurements of the splitting or birefringence of shear waves, is a commonly-used and powerful technique to study seismic anisotropy, particularly in the mantle (Long & Silver, 2009;Savage, 1999). Splitting analysis has previously been applied to eastern North America at different scales (e.g., Aragon et al., 2017;Fouch et al., 2000;Link et al., 2022;Long et al., 2016;Lopes et al., 2020;Yang et al., 2017). Long et al. (2016) conducted SKS splitting analysis on USArray Transportable Array stations across eastern North America and found imperfect correlations between the fast axis direction and absolute plate motion (APM). Notably, Long et al. (2016) and White-Gaynor and Nyblade (2017) pointed out an orogen-parallel fast direction pattern beneath the Appalachian orogen, which suggested substantial contributions from crust and lithospheric mantle anisotropy. In contrast, Yang et al. (2017) conducted a shear wave splitting study at a similar scale and proposed a more APM-parallel overall fast direction beneath eastern North America, with little contribution from lithospheric deformation associated with past tectonic events. Applying the SKS splitting technique to a more densely spaced seismic array in the central Appalachians, Aragon et al. (2017) confirmed the presence of an orogen-parallel fast direction pattern beneath the Appalachians, with a sharp lateral transition in splitting behavior along the array, which they attributed to likely contributions from the lithosphere. To the north, Lopes et al. (2020) found generally APM-parallel fast directions beneath southern New England, but found that the splitting delay times decreased drastically from west to east. They suggested that weaker splitting beneath the eastern portion of the study area may result from either destructive interference from contributions from the lithosphere or from a change in the flow direction in the asthenosphere, possibly associated with the presence of the Northern Appalachian Anomaly (NAA), a prominent low velocity anomaly centered beneath New England (e.g., Levin et al., 1995Levin et al., , 2018;;Li et al., 2003;Tao et al., 2020). The ambiguity in the interpretations involved in shear wave splitting studies is rooted in the fact that shear wave splitting is a path-integrated method, so it intrinsically lacks the vertical resolution and it cannot uniquely distinguish contributions from lithosphere and asthenosphere anisotropy (Park & Levin, 2002).

One way to complement shear wave splitting observations and obtain constraints on the vertical distribution of anisotropy is via anisotropic Ps receiver function (RF) analysis, which is based on the observation of directionally dependent variations of Ps converted phases (Frederiksen & Bostock, 2000;Levin & Park, 1997;Park & Levin, 2016a). Anisotropy-aware RF analysis has been applied to both subduction zones (e.g., Bar et al., 2019;Haws et al., 2023;Wirth & Long, 2014) and continental settings (e.g., Ford et al., 2016;Schulte-Pelkum & Mahan, 2014b;Wirth & Long, 2014) to investigate potential multi-layered anisotropy beneath those study regions. This technique has also been applied to eastern North America (e.g., Frothingham et al., 2022;Li et al., 2021;Long et al., 2017;Yuan & Levin, 2014). Yuan and Levin (2014) conducted both SKS splitting and the anisotropic RF analyses on three long-running stations in northeastern US, and found coherent lower lithosphere anisotropy with a fast axis nearly orthogonal to the strike of the Appalachian orogen and APM-consistent asthenosphere anisotropy. Long et al. (2017) compared anisotropic RFs from stations in the Grenville Province and in the Appalachian domains and identified dissimilar anisotropic features in the crust and lithospheric mantle, LUO ET AL.

3 of 23 which suggested that different tectonic processes had modified the lithosphere beneath these two regions. Li et al. (2021) examined RFs of 62 stations in northeastern North America, focusing on discontinuities in the upper mantle. They identified multiple negative velocity gradients accompanied by contrasts in anisotropy in the lithosphere at 60-100 km depth beneath most stations, and they found the lithosphere to be thinnest beneath central New England, where the NAA is centered. Despite the valuable insights on anisotropic layering beneath individual stations gleaned from these studies, their interpretations have been hindered by an inability to establish clear correlations between the features observed beneath each station. As a result, a comprehensive investigation of any given feature becomes exceedingly challenging. This limitation can be partially attributed to the relatively large distance between the seismic stations used in these studies. It is difficult to trace localized structures such as shear zones and faults across stations that are more than 70 km apart (the average spacing of TA stations), especially because the strike and dip of those features are not known, and structures vary abruptly perpendicular and parallel to the orogen. This problem can be mitigated by using a denser seismic array, which potentially allows the identification of individual features across several stations, thus providing constraints on the shape, orientation, and lateral extent of anisotropic features, and aiding the interpretation of their nature and origin.

The Seismic Experiment for Imaging Structure beneath Connecticut array (SEISConn; Long & Aragon, 2020), a dense broadband seismic array in southern New England, provides an opportunity to use the anisotropic RF technique to investigate anisotropic structures beneath this region in detail. The presence of multiple layers of anisotropy beneath the SEISConn array was proposed in a previous SKS splitting study (Lopes et al., 2020) as well as a previous RF study focusing on radial component seismograms (Luo et al., 2021). In this study, we unravel the anisotropic complexity beneath southern New England in detail by analyzing both radial and transverse components of RFs, performing harmonic decomposition to extract backazimuth-dependent terms, and conducting Bayesian inversion to resolve detailed anisotropic geometries and explore potential interpretations. We identify several anisotropic interfaces in the crust and lithospheric mantle that can be traced beneath adjacent stations and that may represent relict shear zones or fault zones related to past tectonic events. Additionally, two deeper anisotropic interfaces at 80-160 km depths are consistently observed beneath the entire array, which may reflect a change in the dominant mantle flow direction beneath the study region. Findings from this detailed regional study can be extrapolated to other segments of the Appalachian orogen and to other collision zones to understand how orogenesis contributes to lithospheric architecture and how deformation may be recorded by seismic anisotropy. This study also demonstrates how these anisotropic features might be modified or overprinted by younger processes.

The SEISConn array consisted of 15 broadband seismometers, named CS01 to CS15 from west to east, deployed with ∼10 km spacing across northern Connecticut from 2015 to 2019 (Figure 1a; Long & Aragon, 2020). The array traverses the eastern edge of Laurentia near station CS05, which separates the exposed Proterozoic Grenville crust to the west and the Gondwana-derived Moretown terrane crust accreted during the Ordovician Taconic orogeny to the east. Stations CS06-CS09 in the central portion of the array are situated above the Hartford basin, which was opened by the Mesozoic rifting event. To the east, there is likely a buried suture associated with the accretion of Ganderia during the Late Ordovician-Early Silurian Salinic orogeny. Farther east, near station CS14, the array traverses the exposed western edge of Avalonia, another Gondwana-derived terrane that was accreted to the composite Laurentia during the Latest Silurian-Devonian Acadian orogeny.

The main goal of the SEISConn deployment is to understand how the crust and mantle lithosphere beneath southern New England have been modified by past orogenic cycles and younger processes through investigations of the seismic velocity structure, discontinuities, and anisotropy beneath this region. A detailed crustal velocity model has been constructed by Gao et al. (2020) using the full-wave ambient noise tomography technique. They found a low-velocity mid-crust, likely reflecting radial anisotropy produced during the Mesozoic rifting, and a high-velocity lower crust beneath the Hartford Basin, which might be related to magmatic underplating during the Central Atlantic Magmatic Province event. Traditional 1-D RF analysis using the converted Ps phase (Luo et al., 2021) showed a ∼15 km Moho depth offset between station CS03 and CS04, as well as multiple discontinuity features in the crust and upper mantle that might be associated with subduction and terrane accretion events during the Appalachian orogenic cycle. A lithospheric step of similar size is also suggested beneath southern New England by RF analysis using the converted Sp phase (Goldhagen et al., 2022). A follow-up study using a more LUO ET AL. 2022) also further constrained the presence of a relict slab in the lithosphere, possibly subducted during the closure of the Rheic Ocean, and a localized low velocity anomaly that might be associated with the NAA in the upper mantle. The potential presence of the NAA beneath southeastern New England was also implied by SKS splitting analysis (Lopes et al., 2020), which showed a generally APM-parallel fast axis orientation with a gradually decreasing splitting delay time toward the east, near the likely southern margin of the NAA. However, the potentially multi-layered upper mantle anisotropy inferred in this study was not explicitly modeled by applying a grid-search and curve-fitting technique to the splitting results, due to insufficient backazimuthal coverage of SKS phases (Lopes et al., 2020). The scattered wave analyses (Luo et al., 2021(Luo et al., , 2022) ) focused on layers with isotropic velocity contrast, and the SKS split ting study (Lopes et al., 2020) focused on the path-integrated signature of upper mantle anisotropy beneath southern New England. Although multiple layers of anisotropy with varying geometries were suggested by several of those studies, the depths and characteristics of individual anisotropic layers could not be explicitly resolved. This knowledge gap can be filled with the anisotropy-aware RF analysis that we carried out here.

Ps receiver function analysis is based on identifying S-waves converted from incident P-waves at structural discontinuities with changes in material properties (Langston, 1979;Rondenay, 2009). We rotated seismograms to the L-Q-T coordinate system, where L is approximately parallel to the direct P wave, Q is orthogonal to L in the vertical plane, and T is orthogonal to both L and Q. The radial and transverse components of RFs are calculated by deconvolving the L component from the Q and T components of seismograms, respectively, using the multiple-taper-correlation method (Park & Levin, 2000, 2016b). (Note that for simplicity, we use the more conventional term "radial RF" for the RFs calculated based on the Q-component, as in Ford et al. (2016)). A detailed description of our data processing procedures and quality control criteria is in Text S1 of the Supporting Information S1. As shown in Figure 1b, earthquakes used as sources are found at a wide range of backazimuths, resulting in incident waves arriving from a range of directions.

When the seismic structure beneath a receiver is made up of horizontal layers with only isotropic velocity contrasts (Figure 2a), the P-to-S conversion will only produce vertically polarized SV waves whose amplitudes do not vary with backazimuth. However, when the interface is dipping (Figure 2b) or when one of the layers is anisotropic with a plunging (Figure 2c) or horizontal (Figure 2d) axis of symmetry, both SV and horizontally polarized SH waves are produced, and their polarities and/or amplitudes vary systematically with backazimuth (Levin & Park, 1997). Specifically, if a slow layer overlies a fast layer across a west-dipping interface (Figure 2b), the P-to-SV conversion will be stronger for the incident wave coming from the west and weaker for the incident wave coming from the east. If the interface is horizontal but the lower layer is anisotropic with an east-plunging fast axis (Figure 2c), the P-to-SV conversion will be stronger for the incident wave coming from the east and weaker for the incident wave coming from the west. If the fast axis in the lower layer is instead horizontal and oriented E-W (Figure 2d), the P-to-SV conversion will be stronger for the incident wave coming from the east or west and weaker for the incident wave coming from the north or south. For all the cases mentioned above, if the positions of the upper and lower layers are switched, the directional dependence of the P-to-SV conversion will also be inverted. The directionally-dependent P-to-SV conversion is always accompanied by directionally-dependent P-to-SH conversion. The directional dependence of the P-to-SH conversion will be systematically phase shifted by +90° from that of the P-to-SV conversion if the structure is horizontally layered and can be fully modeled by transverse isotropy (TI) with a cylindrical axis of symmetry (Park & Levin, 2016a).

Anisotropic RF analysis can reveal the presence of dipping interfaces and anisotropic layers by exploring the event backazimuth-dependent variations in the radial component RFs, which record the P-to-SV conversion, and the transverse component RFs, which record the P-to-SH conversion. Of course, the Earth structure is usually much more complicated than the two-layer examples shown in Figure 2. There can be multiple anisotropic layers with different anisotropic geometries separated by dipping interfaces, and the anisotropy involved can be more complicated than TI. As a result, it is often difficult to analyze the backazimuth-dependent variations of RFs directly because the signals generated from different features are blended. One way to more clearly characterize dipping interfaces and anisotropic orientations is to apply the harmonic decomposition technique (Bianchi et al., 2010;Olugboji et al., 2016;Park & Levin, 2016a;Shiomi & Park, 2008). This technique uses linear regression to model the behavior of both the radial and transverse RFs, extracting the constant term (the zeroth-order harmonic term) that does not vary with backazimuth as well as the higher-order harmonic terms (e.g., sin(baz), cos(baz), sin(2*baz), and cos(2*baz)) that do vary with backazimuth.

Figure 3 provides three synthetic examples to demonstrate the capabilities of harmonic decomposition. The first model (Figure 3a) consists of two layers over a half-space, with a 30-km thick isotropic upper crust and a 5-km thick anisotropic lower crust that has the same isotropic velocity as the upper crust, but with an east-plunging +10% fast axis of anisotropy. The Moho is at 35 km depth, and the isotropic mantle half-space has faster seismic velocities than the crustal layers above. We calculate synthetic RFs using the PyRaysum program (Bloch & Audet, 2023), which generates ray-theoretical seismograms for incident plane waves propagating through seismic models with dipping interfaces and layered anisotropy (Frederiksen & Bostock, 2000). Figure 3d displays the calculated synthetic radial and transverse RFs in backazimuth gathers, with trivial random noise added. There are clear backazimuth-dependent variations of both radial and transverse component RFs at arrival times that correspond to 30 and 35 km depth, representing the top and bottom interfaces of the anisotropic lower crust. The variation in arrivals at 30 and 35 km depth on the radial component can be modeled as first-order positive and negative sine functions of backazimuth, respectively, and the patterns in the transverse component are correspondingly +90° phase shifted. These patterns can be extracted using the harmonic decomposition technique and are visible in the sin(baz) harmonic term (Figure 3g). The isotropic velocity increase across the Moho is visible in the constant term at 35 km depth, as expected. The harmonic decomposition also shows that the east-plunging anisotropic axis also produces some minor cos(2*baz) signals in addition to the prominent sin(baz) signals. The synthetic RFs are computed based on a horizontally layered TI model, and therefore no significant energy is present in the unmodeled portion of RFs (Figure S1a in Supporting Information S1). Similarly, Figures 3e and3h show the synthetic RFs in backazimuth gathers and the corresponding harmonic decomposition, respectively, calculated from a seismic model with an anisotropic lower crust with a horizontal +10% fast axis of anisotropy (Figure 3b). A completely isotropic seismic model with a dipping interface (Figure 3c) can also produce backazimuth-dependent variation in RFs (Figure 3f), which can be modeled as harmonic terms (Figure 3i). The presence of a dipping interface can cause backazimuth-dependent variation in the direct P phase, which manifests as a substantial sin(baz) signal at 0 s delay time in Figure 3i. In addition, there are larger unmodeled harmonic signals associated with the dipping interface compared with a model with anisotropy in horizontal layers (Figure S1c in Supporting Information S1).

For the synthetic example shown in Figure 3, we can qualitatively infer the potential anisotropic geometries and dipping interfaces involved in the seismic model (Figures 3a-3c) by visually inspecting the harmonic terms (Figures 3g-3i). In the case of particularly prominent anisotropic features that have significantly larger amplitudes than other signals, it is feasible to derive quantitative constraints on the strength and orientation of the anisotropy. One way to do this is via forward modeling, in which we fit the observed harmonic terms with synthetic RFs calculated for models with parameters based on a combination of prior knowledge and trial-and-error (e.g., Bianchi et al., 2010;Ford et al., 2016;Liu et al., 2015). With many unknown parameters, such forward modeling usually involves a substantial amount of trial-and-error, and it is often difficult to fully explore the tradeoffs among parameters such as the dipping angle of the symmetry axis and the strength of the anisotropy.

These difficulties can be mitigated by applying Bayesian inversion, which infers values of model parameters from a probabilistic perspective, instead of doing extensive trial-and-error forward modeling. A Markov chain Monte Carlo approach has been applied in several previous RF studies (e.g., Bissig et al., 2021;Bodin et al., 2013Bodin et al., , 2016;;Wirth et al., 2016). Here, we use a Bayesian framework with efficient ray-theoretical forward calculations (Link et al., 2020) to invert the anisotropic and dipping interface model that best fits the observed constant and harmonic terms. One advantage of fitting the harmonically decomposed RFs, instead of traditional RFs or raw seismograms, relates to the fact that the harmonic terms are extracted based on the assumption of lateral homogeneity and TI anisotropy with a cylindrical axis of symmetry; therefore, patterns of signals that cannot be produced by dipping interfaces and TI anisotropy (e.g., scattering due to lateral heterogeneity, effects from varied incident wave slowness due to different epicentral distances, anisotropy that cannot be approximated using a cylindrical axis of symmetry) are suppressed. Therefore, applying Bayesian inversion to the harmonically decomposed RFs has the potential to infer the anisotropic seismic model under the TI assumption more faithfully, without attempting to fit other unrelated signals and artifacts.

We first present results from a single station of the SEISConn array, station CS07, both as a case study that illustrates how we apply our analysis to real data, and as an example of a station that exhibits strong, clear signals from crustal anisotropy that can be modeled using our Bayesian inversion technique. Next, we present results from across the SEISConn array, with a focus on identifying signals that are laterally coherent across multiple stations.

CS07 is a station deployed in the Hartford Basin in the central part of the SEISConn array (Figure 4a). The station has good event backazimuth coverage (Figure 4b), with 149 events selected for RF analysis. On the radial component, the Moho conversion is visible as a positive signal at 3.5-4.5 s that is coherent across all backazimuths. Perhaps a more eye-catching feature is a set of signals at 0.5-1.0 s whose amplitudes and polarities vary with backazimuth. There is a negative signal followed by a positive signal at 0.5-1.0 s in the 45°-145° backazimuth range, which is reversed to a positive signal followed by a negative signal at 185°-355° backazimuth range. Signals at the backazimuth range in between (i.e., 145°-185° and 355°-45°) have weaker amplitudes at 0.5-1.0 s. This backazimuthally dependent pattern also manifests in the transverse component with a +90° phase shift, such that the negative followed by positive signals occur at 135°-235° while the positive followed by negative signals occur at 275°-85°. This correlation between the radial and the transverse RFs suggests that this backazimuth-dependent variation pattern can be generated by TI anisotropy. Indeed, this feature is identified with prominent sin(baz) signals in the harmonic decomposition (Figure 4c). We observe a large negative signal at ∼3 km depth followed by another large positive signal at ∼6 km depth in the sin(baz) term; these likely represent the sharp top and bottom interface of a ∼3 km thick anisotropic layer. The pattern in the sin(baz) term can be explained by either a west-plunging fast axis or an east-plunging slow axis of symmetry in the anisotropic layer. This ambiguity is typically difficult to resolve in practice; the difference between the two cases lies in the smaller-amplitude accompanying signals on the cos(2*baz) term. Specifically, a west-plunging fast axis will yield negative followed by positive signals, but an east-plunging slow axis will yield positive followed by negative signals. Uncommonly, the feature we observe beneath CS07 is prominent enough that the accompanying signals in the cos(2*baz) term can be recognized as positive followed by negative at the corresponding depths. Therefore, this anisotropic layer likely has an east-plunging slow axis of symmetry. Above this anisotropic layer, there is a slower layer, indicated by a positive signal at ∼3 km depth in the constant term, which represents an isotropic velocity increase with depth. This slow layer likely represents the sedimentary units of the Hartford Basin, given the location of CS07 and the depth of the interface. The Moho also shows up clearly at ∼26 km depth in the constant term, as expected. Minor signals are present in the unmodeled RFs (Figure S1d in Supporting Information S1), which can arise from effects such as lateral heterogeneities, noise, and irregular event distribution. RF results for all individual SEISConn stations, similar to the example from the CS07 shown in Figure 4, are shown in Supporting Information S1 (Figures S2-S15).

We take advantage of the large amplitude of the conversions associated with the shallow crustal anisotropic layer beneath CS07 and apply Bayesian inversion to obtain more quantitative constraints on its anisotropic geometry as well as its dipping top and bottom interfaces (Figure 5). The inversion is conducted on the harmonic terms derived from RFs that are band-passed filtered between 0.02 and 1.5 Hz because signals from the shallow crustal S1 in Supporting Information S1; Kennett & Engdahl, 1991). Unmodeled terms are shown in Figure S1d of the Supporting Information S1.

interfaces can be well isolated with this frequency band, as shown in Figure 4. An overview of the parameter space setup is shown in Table 1, and more details of the inversion procedure are available in Supporting Information S1. The mean values of the Bayesian statistics (Figure 5a) can reconstruct synthetic RFs that capture all major crustal signals observed in the raw RFs of CS07, as evident in the resemblance between Figures 4b S1 of the Supporting Information S1 and assuming vertical incidence.

Lower/upper bounds Note. The parameters for P-wave velocity and density (v P and ρ) are kept fixed in all inversion runs due to limited sensitivity for these parameters. Thickness (H), P-wave to S-wave velocity ratio (v P /v S ), anisotropic strength (a), the trend of slow axis from north (ϕ), its plunge from the horizontal (θ), the interface dip (δ), and the strike of the layer interfaces (γ) are parameters which are searched for in the inversion and vary between the boundaries separated by the dashed lines. Bayesian inversion is performed for several different initial starting models, where the strength of anisotropy and fast axis direction are chosen randomly from the given interval, while the remaining free parameters are set to the parameters given in the lower half of the table. The inversion resolves an east-plunging slow axis in the second layer with >20% anisotropy situated between a slower top layer and a (nearly) isotropic deeper, third layer with similar velocities, consistent with our qualitative interpretation based on visual inspection. Furthermore, the inversion sheds light on the potential existence of a west-plunging slow axis with ∼10% anisotropy in the top layer, associated with the sedimentary units of the Hartford Basin, which was not as obvious from visual assessment. The more general application of such quantitative analysis to other, less prominent features is hindered by the ambiguity and non-uniqueness involved in the harmonic analysis. For example, the difference between an east-plunging slow axis and a west-plunging fast axis lies in very minor signals in the cos(2*baz) term, which are likely buried in noise and undistinguishable in typical cases. In the following analysis of the entire array, in which most features do not have as significant amplitudes as the shallow crustal features observed beneath CS07, we apply a more qualitative approach to interpretation.

Figure 6 shows the compiled results of harmonic analysis applied to data from the entire SEISConn array using a bandpass filter between 0.02 and 0.1 Hz. The constant RF profile is generally consistent with results from previous isotropic analyses (Luo et al., 2021(Luo et al., , 2022)), showing a large Moho depth offset between stations CS03 and CS04. In the profiles for the higher-order terms of the harmonic analysis, we can observe a plethora of signals with varied amplitudes and polarities, indicating an abundance of sharp anisotropic contrasts and/or dipping interfaces beneath this region. In general, harmonic signals in the crust have larger amplitudes than those at greater depths, and harmonic signals beneath the western portion of the array have larger amplitudes than those beneath the eastern portion of the array. We can identify several features, both in the crust and in the upper mantle, that have coherent harmonic signals across several adjacent stations.

In order to first identify signals within the crust and shallowest mantle, we show in Figure 7 a set of profiles similar to those in Figure 6 but focused on the upper 80-km depth range beneath the study region and with major anisotropic and/or dipping features labeled. Data used to construct these images were bandpass filtered between 0.02 and 1.5 Hz to emphasize shallower, smaller-scale features. The largest harmonic signals are associated with Feature #1 in the sin(baz) profile (Figure 7). This feature is identifiable beneath stations CS05 to CS09, with larger amplitudes beneath CS05, CS06, and CS07 and smaller amplitudes beneath CS08 and CS09. As discussed in Section 3.1, this negative followed by a positive signal pattern in the sin(baz) term likely denotes an anisotropic layer in the shallow crust with east-plunging slow axis anisotropy (Figure 5). This feature displays a zig-zag shape in cross section, and the top interface of this anisotropic layer aligns with the bottom of the Harford basin, except for the down-going branch beneath CS05. Feature #2 consists of positive signals near or before 0 s time and negative signals at shallow depths in the sin(baz) terms beneath CS01, CS02, and CS03. The signal around or even before 0 s delay time in harmonic RFs indicates backazimuth-dependent variation of the direct P arrival amplitude, which suggests the presence of dipping interfaces (Figure 3i; Schulte-Pelkum & Mahan, 2014a). In addition to the signal at or near 0 s delay time, a dipping interface will also result in the variation of P-to-S conversion amplitude at the depth of the interface, with opposite polarity to the variation of direct P. Therefore, Feature #2 likely represents a series of east-or west-dipping interfaces at 0-5 km depth beneath CS01, CS02, and CS03. Feature #3 is an elongated west-dipping interface that expresses itself in the cos(baz) term, which suggests a generally north-or south-plunging anisotropic axis. There is no accompanying interface with opposite polarity in the cos(baz) term, which implies that this interface can be either the top or the bottom interface of an anisotropic layer, with the other interface being diffuse and thus not observable by RF analysis. Features #4 and #5 in the sin(baz) profile align with the Moho geometry shown in the constant RF profile. The negative and positive sin(baz) signals are likely produced by the west-(in the west) and east-(in the east) dipping Moho discontinuities.

In the cos(baz) term, we also observe positive signals associated with Feature #4 and negative signals associated with Feature #5; these provide information about the Moho geometry in the north-south direction, outside the profile plane. The character of these signals suggests that the Moho beneath the western portion of the profile is northwest-dipping, whereas the Moho beneath the eastern portion of the profile is southeast-dipping. Feature #6 is a west-dipping interface with a positive velocity gradient in the upper mantle that only shows up in the constant RF as a series of minor positive signals from ∼56 km depth beneath CS15 to ∼77 km depth beneath CS06. This isotropic feature is also observed in the previous isotropic RF and scattered wavefield migration studies (Luo et al., 2021(Luo et al., , 2022)).

In order to identify features deeper in the mantle, we show in Figure 8 profiles similar to those in Figure 6 with data filtered at lower frequencies (0.02-0.5 Hz); in this view, deeper and broader structures come into focus. Feature #7 represents an extensive anisotropic layer beneath the Moho with a generally NE-SW fast or NW-SE slow horizontal axis of symmetry, indicated by a series of positive signals followed by negative signals in the sin(2*baz) term. Also, the presence of positive cos(baz) signals related to Feature #7 may imply varied anisotropic geometries within the inferred anisotropic layer. Deeper in the upper mantle, features #8 and #9 can be coherently observed beneath almost the entire array, and they are not artifacts associated with the arrival of multiple Moho reflection signals (magenta dashed lines in Figure 8). Feature #8 is a flat-lying feature observed at 80-100 km depth in the cos(baz) term, suggesting that it is either a north-or south-dipping interface, or an anisotropic interface associated with a north-or south-plunging anisotropic axis. Feature #9 is a west-dipping feature observed at 80-160 km depth in the sin(2*baz) term, associated with either a NE-SW fast or NW-SE slow horizontal anisotropic axis. The fact that these deep features are not as well imaged in higher frequencies (Figure 6) indicates that they represent more diffuse, gradational boundaries compared to S1 in Supporting Information S1; Kennett & Engdahl, 1991). For each trace, the black curve represents the bootstrap mean, and the pale green band shows ±1 standard deviation range of the trace derived from the bootstrap test. Positive and negative regions under the curve that are beyond one standard deviation are filled with red and blue, respectively. The top left panel shows the constant RF profile, and the panels at the middle and to the right show harmonic RF profiles for higher-order terms, as indicated by labels at the top. The spacing between traces reflects the physical distances between stations along the profile in kilometers. The aspect ratio of the plots in this figure is 1:1. Unmodeled terms are shown in Figure S24 of the Supporting Information S1.

the sharper boundaries observed in the crust. Also, the converted signals associated with these deep features are likely further modified and defocused when traveling through the shallower heterogeneous and anisotropic structures.

Figure 9 shows the major features observed in the harmonic RFs along with their preferred interpretations. The signals appearing in the higher-order terms of harmonically decomposed RFs usually indicate the presence of a dipping interface, an interface with anisotropic contrast, or both (Park & Levin, 2016a). Because the extents and shapes of these features beneath the SEISConn array are constrained by the dense station spacing, we can use the RF variation patterns extracted by harmonic decomposition along with a priori knowledge about the geologic/ tectonic background to interpret specific crust and upper mantle structures (Figure 9). The set of features we document, with drastically different anisotropic geometries, is highly unlikely to have been formed by a single tectonic event. Instead, they likely reflect the varied deformation patterns of several different past tectonic events and/or present-day processes. In the following section, we will discuss in detail the potential origins and implications of the anisotropic and/or dipping interfaces beneath southern New England, and how the insights from this study can have wider applications and broader implications beyond this specific region.

Beneath stations CS01, CS02, and CS03 in the westernmost portion of the array, we observe a series of dipping interfaces at very shallow depths (0-5 km), identified as Feature #2 in the sin(baz) term (Figure 7). They can be either west-or east-dipping. Resolving this ambiguity solely based on the observed RFs proves to be challenging, given that these interfaces are so close to the surface and that the signals in the constant term can be severely affected by artifacts such as side-lobes of the incident P phase. The geographic location of this feature is on the Laurentian margin, to the west of the suture between Laurentia and the Gondwana-derived Moretown terrane accreted during the Ordovician Taconic orogeny (Karabinos et al., 2017;Macdonald et al., 2014). To the north of our study area, a previous seismic reflection study traversing eastern New York, Vermont, and New Hampshire (the COCORP line; Ando et al., 1984) revealed numerous east-dipping interfaces beneath the Laurentian margin in the New England Appalachians that may be interpreted as a system of Taconic thrust faults (Stanley & Ratcliffe, 1985). Feature #2 may then signify the presence of this east-dipping fault system beneath southern New England.

We find that the Moho beneath the SEISConn array is shallowest beneath the Hartford Basin in the central portion of the array, and it gradually deepens to the west and to the east, consistent with our previous results (Luo et al., 2021(Luo et al., , 2022)). Because of the large isotropic velocity contrast between the crust and the upper mantle, the Moho manifests clearly in the constant RFs as positive signals (Figures 678). Furthermore, the west-and east-dipping portions of the Moho discontinuity result in negative and positive signals beneath the western and eastern portions of the array, respectively (Features #4 and #5). Feature #4, associated with the west-dipping Moho, may extend all the way to beneath CS02, beyond the prominent Moho depth offset located between CS03 and CS04 (Luo et al., 2021). The deeper Grenville Moho to the west of the offset is flat lying, depicted as a horizontal black line in Figure 9. Thus, unlike the west-dipping shallower Moho, the deeper Moho does not generate any coherent higher-order harmonic signals. In this way, the two distinct Moho interfaces can be clearly identified and distinguished beneath CS03 and CS02, which was not possible using our previous traditional RF analysis (Luo et al., 2021). This observation is another piece of evidence supporting the idea of a doubled Moho in the vicinity of the Moho depth offset, with the shallower Moho in the east extending above the deeper Moho in the west, as suggested by Luo et al. (2022). The observation of a doubled Moho has important implications on the origin of the Moho depth offset in southern New England. It favors a mechanism involving the overthrusting of the rifted Grenville crust and the Moretown terrane from the east onto undeformed Grenville in the west, along a reactivated Proterozoic rift fault (Luo et al., 2022). At the other end of the profile, the sin(baz) and cos(baz) signals associated with Feature #5 suggest that the Moho beneath the eastern portion of the array dips not only to the east but also to the south. This observation points to locally thicker crust to the southwest of the SEISConn array in central Rhode Island, which was also hinted at in previous larger-scale RF studies (e.g., Li et al., 2018).

Using the scattered wavefield migration technique, Luo et al. (2022) also identified a prominent west-dipping interface in the upper mantle cutting across the entire profile that represents an isotropic velocity increase with depth. This interface was interpreted as the Moho of a relict slab from a Paleozoic subduction event that was retained in the lithosphere after the deeper part had broken off. This feature is substantially weaker and less coherent in the traditional RF analysis (Luo et al., 2021), which can be attributed to a combination of the destructive interference between the top of the slab and the slab Moho as well as the possible partial eclogitization of the oceanic crust. Here, in the harmonically decomposed RFs, we can only glean the presence of this feature (Feature #6) in the constant RF profile (particularly in the higher frequency image in Figure 7). The higher-order harmonic signals associated with dipping interfaces are usually much smaller than the constant signals (e.g., Features #4 and #5). The velocity contrast across this dipping interface is likely too small for its harmonic signals to stand out and be distinguishable from other signals and background noise. This dipping interface may be associated with the formation of an anisotropic zone above it (Feature #7), which will be further discussed in Section 4.3. 7 and8. Black solid lines denote flat lying isotropic interfaces; cyan solid lines denote dipping isotropic interfaces; cyan dash-dotted lines denote interfaces that have both an isotropic seismic velocity contrast and an anisotropic contrast; cyan dotted lines denote interfaces that only have an anisotropic contrast. Arrows denote the inferred anisotropic fabric orientations or mantle flow directions. Crosses and circles denote in-and-out of the plane directions (north and south in this case). The sketch is not scaled to allow for the visibility of all features, but approximate depths of key interfaces are noted to the right.

Our Bayesian inversion results (Figure 5) suggest that anisotropy with a steeply west-plunging slow axis may be present in the strata of the Hartford Basin, one of the abandoned rift basins created during the Mesozoic rifting (Hubert et al., 1992;Withjack et al., 2012). The Hartford Basin is a half-graben with a well-developed eastern border fault, and the layers within the basin are generally tilted to the east (Wise, 1992). The observed anisotropy can be explained by layering in the Hartford Basin. Fine layering in sedimentary basins can result in TI anisotropy with a layer-perpendicular axis of symmetry when layer thicknesses are smaller than the wavelength (Backus, 1962;Yan et al., 2016). The Hartford Basin consists of a series of layers with different seismic velocities including several basaltic flows (Philpotts & Martello, 1986). The thicknesses of those layers are much less than the wavelengths used in the RF analysis, which are on the order of kilometers, so they can result in anisotropic behaviors in RFs with a symmetry axis perpendicular to the layers. The observation of a generally west-plunging axis is consistent with the east-dipping layers in the Hartford Basin. The plunge has a mean value of 56° from the Bayesian inversion (Figure 5), which is in striking agreement with the ∼55° dip of the eastern border fault in Connecticut (Digman, 1950) and complements the 15-35° dip of the strata (Wise, 1992). Furthermore, a mean axis trend of 110° is also consistent with the NNE strike of the eastern border fault and the basin layers (Wise, 1992).

Beneath the sedimentary strata of the Hartford Basin, there is a prominent layer with strong anisotropy as well as sharp top and bottom interfaces. This feature, observed beneath stations CS05-CS09, has a thickness of 3-4 km and a lateral extent of ∼50 km (Feature #1 in Figures 7 and8). The layer likely has an east-plunging slow axis of anisotropy, with >50° plunge and >20% anisotropy, as investigated in detail beneath CS07 (Figures 4 and5). Anisotropy observed in the shallow upper crust with low confining pressure may be dominated by aligned cracks whose orientations are controlled by the stress field (Almqvist & Mainprice, 2017). Cracks orthogonal to the maximum compressive axis will be preferentially closed, and the seismic waves will propagate faster along this direction (Crampin & Chastin, 2003;Crampin & Lovell, 1991). Consistent with this mechanism, the observed fast directions of upper crustal azimuthal anisotropy across the contiguous U.S. are generally aligned with the maximum horizontal compressive stress direction (Lin & Schmandt, 2014). However, with a lateral extent of only ∼50 km, we suggest that Feature #1 is likely not associated with the present-day regional stress field in the upper crust beneath southern New England; furthermore, its plunging geometry is not what would be expected for an aligned crack mechanism.

We instead favor an alternative explanation for this feature, which involves foliated metamorphic rocks and the CPO of crustal minerals. The CPO of mica and/or amphibole is usually invoked as the major contribution to substantial seismic anisotropy observed in the mid to lower crust (Lloyd et al., 2009;Rudnick & Fountain, 1995;Shapiro et al., 2004). Feature #1 is at 0-12 km depth today, but was likely at greater depths in the past. Insights into the history of crustal thickening and subsequent erosion and exhumation beneath New England come from geological constraints. Specifically, Hillenbrand and Williams (2021) proposed the existence of a long-lived orogenic plateau in southern New England during the Acadian and Neoacadian orogeny, which collapsed due to reduced compressive stresses after the Neoacadian orogeny. At present, surface rocks in southern New England were at larger depths in the past and were exhumed to the surface due to erosion and crustal thinning during and after the collapse of the Acadian orogenic plateau. Near the locality of Feature #1, the paleodepths of surface rocks are around 16-22 km (Hillenbrand & Williams, 2021), so this anisotropic layer would have been at 16-34 km depth before the proposed plateau collapse. This layer likely consists of metamorphic rocks such as schist and gneiss that are widespread in the nearby Gneiss Dome belts, Connecticut Valley basin, and the Bronson Hill arc (Connecticut Geological Survey, 2013). These metamorphic rocks contain minerals such as mica, amphibole and quartz that can form CPO under ductile deformation in the mid and lower crust (Almqvist & Mainprice, 2017). A hot and weak, perhaps partially molten, mid-crust susceptible to ductile flow was proposed for the Acadian orogenic plateau (Hillenbrand et al., 2022), similar to the mid-crust today beneath the Tibetan Plateau (e.g., Nelson et al., 1996). The anisotropy caused by CPO of mica can be approximated by a TI with a slow axis of symmetry, but deviations from TI are observed for amphibole, sillimanite, and quartz (Ji et al., 2015;Ward et al., 2012). The observation of strong slow-axis anisotropy with minor unmodeled harmonic signals (Figures S1d and S26 in Supporting Information S1) suggests that Feature #1 can be almost entirely described by TI, which may imply that the anisotropy is dominated by CPO of mica. Furthermore, the cleavage plane of mica is usually aligned with the foliation, which can result in seismic anisotropy with a slow axis of symmetry that is perpendicular to the foliation plane and a fast plane parallel to the foliation plane, while amphibole and sillimanite typically have crystallographic axes aligned with the lineation, which instead result in a seismic fast axis of symmetry that is subparallel to the lineation and a slow plane perpendicular to the lineation (Sherrington et al., 2004). Feature #1's well-constrained slow axis of symmetry is therefore consistent with the anisotropy caused by CPO of mica.

The observed east-plunging slow axis of Feature #1 suggests generally west-dipping foliation planes in this anisotropic layer. This west-dipping fabric might be associated with the shear zone resulting from the terrane accretion episodes during the Salinic and/or Acadian orogenies, which likely had west-dipping polarity (Robinson et al., 1998;van Staal et al., 2009). In the case of simple shear near suture faults, metamorphic rocks can develop S-C fabric, with the S-plane forming at ∼45° to the shear direction, and the C-plane parallel to the shear direction (Lister & Snoke, 1984;Lloyd et al., 2009). The anisotropic behavior of Feature #1 is likely controlled by west-dipping foliation planes, which are parallel to the shear direction during west-dipping terrane accretion. Therefore, this layer likely experienced high shear strain during the Paleozoic terrane accretion event and is thus dominated by the shear motion-parallel C-fabric. A similar structure has been observed in the southern Appalachians at the Alleghanian Laurentia-Gondwana suture (Hopper et al., 2017). The Paleozoic fossil fabric associated with Feature #1 has been preserved until today, but the shape of this originally west-dipping layer was likely distorted by Mesozoic rifting and the opening of the Hartford Basin, resulting in the observed zig-zag shape beneath stations CS05-CS09.

Beneath Feature #1, there is a more elongated anisotropic interface identified as Feature #3, which may correspond to the anisotropic mid-crustal layer inferred by Gao et al. (2020). This feature extends from ∼5 km depth beneath CS13 to ∼20 km depth beneath CS03, which corresponds to an Acadian paleodepth of 21-42 km (Hillenbrand & Williams, 2021). This anisotropic feature is also likely associated with fossil CPO of crustal minerals in foliated metamorphic rocks because microcracks will generally be closed at pressures above 200-250 MPa (Ji et al., 2013) or about 10 km depth. Different from Feature #1, Feature #3 has a north-or south-plunging axis, a much smaller amplitude (suggesting weaker anisotropy), and a missing paired top or bottom interface that accompanies the observed interface, which makes the identification of its exact anisotropic geometry and strength impossible. The generally north-south fabric can be plausibly explained by a pure shear scenario with orogen-parallel minimum compressive axis, associated with the orogen-parallel crustal flow in the mid or lower crust during Acadian orogenic plateau collapse (Hillenbrand et al., 2022). These types of fabrics are evident in exposed metamorphic rocks with orogen-parallel fabrics to the north of our study area in southern Vermont and central Massachusetts (Karabinos et al., 2010;Massey et al., 2017). We speculate that the observed interface might be the relatively sharp top interface of the crustal flow zone, perhaps representing the brittle-ductile transition in the thickened Acadian crust (Kohlstedt et al., 1995;Mainprice & Nicolas, 1989), while the bottom interface is more diffuse and thus not observed in the RFs. The observation that Feature #3 is generally west-dipping may then imply either higher paleoelevation and/or more erosion or exhumation in eastern Connecticut than to the west.

An alternative explanation of Feature #3 may involve channel flow and ductile extrusion of the Putnam-Nashoba terrane in southeastern New England (Severson, 2020). Channel flow and extrusion of partially molten weak mid-crustal materials have been proposed for a number of convergent tectonic settings, including the Himalayas, the Canadian Cordillera, and the U.S. Appalachians (e.g., Beaumont et al., 2001;Kuiper et al., 2006;Merschat et al., 2005). In southeastern New England, the Putnam-Nashoba terrane might have extruded to the surface via channel flow during the Acadian orogeny, supported by geological evidence including the presence of partially melted rocks between lower grade metamorphic rocks above and below, and a transition of shear movements from normal to increasingly reverse toward the southeast (Severson, 2020). Feature #3 can be traced to the surface at the station CS13, which is near the western edge of the Putnam-Nashoba terrane. It is possible that Feature #3's anisotropy arises from the CPO of crustal minerals in the flow channel of the Putnam-Nashoba terrane. However, it is difficult to reconcile the eastward flow/extrusion with the observed north-or south-plunging anisotropic axis of Feature #3.

Continental scale shear wave splitting studies have reported generally E-W or NE-SW oriented fast axes in the northeastern US, which suggest major contributions from plate motion parallel shearing of the asthenospheric mantle (Long et al., 2016;Yang et al., 2017), with the fast axis of olivine aligned with the flow direction (Zhang & Karato, 1995). Nevertheless, regional scale variation in splitting parameters is observed beneath southern New England, with consistent plate motion parallel ENE-WSW fast directions but decreasing delay times toward the eastern end of the SEISConn array (Lopes et al., 2020). From our RF harmonics, we can identify two anisotropic layers with NE-SW trending symmetry in the upper mantle whose characteristics can potentially account for the decreasing delay times observed in the splitting analysis. One is a well-defined layer in the lithosphere associated with Feature #7, and the other is a zone in the asthenospheric mantle sandwiched between Feature #8 and Feature #9.

Feature #7 corresponds to a ∼20 km thick anisotropic layer in the lithospheric mantle with well-defined top and bottom interfaces in the sin(2*baz) term (Figure 8). The harmonic signals observed in the sin(2*baz) term can be explained by either a NE-SW fast axis or a NW-SE slow axis. Given that seismic anisotropy in the upper mantle is usually attributed to olivine CPO, which under the dry conditions associated with continental lithospheric mantle (Peslier et al., 2017) has a fast axis roughly aligned with the shear direction (Karato, 2008;Long & Silver, 2009), Feature #7 likely corresponds to a generally NE-SW fast axis, reflecting fossil olivine CPO induced by a past episode of NE-SW directed shear. This shear zone overlies the relict slab lodged in the lithosphere (Feature #6; Luo et al., 2022) and may be associated with the deformation during the subduction of this slab in the Paleozoic.

As discussed in Luo et al. (2022), the timing of this subduction event is not well constrained; it could have been during the Late Devonian Neoacadian orogeny, which does not have a clear record in southern New England (Kuiper et al., 2017;van Staal et al., 2009), or during the Carboniferous-Permian Alleghanian closure of the Rheic Ocean, whose subduction polarity is controversial (Domeier & Torsvik, 2014;Hermes & Murray, 1988;Michard et al., 2010Michard et al., , 2023;;Nance & Linnemann, 2008;Nance et al., 2012). The observation that the shear zone, which lies only ∼15 km deeper than the continental Moho, smoothly undercuts the prominent Moho depth offset between CS03 and CS04 (Figure 7) may suggest that the formation of the shear zone, and thus the subduction event, postdated the formation of the Moho depth offset. Hillenbrand et al. (2021) proposed that the Moho offset was mostly formed during the collapse of the Acadian orogenic plateau at ca. 330-310 Ma; this timing may favor an Alleghanian subduction origin for the relict slab. Regardless of the timing of subduction, the NE-SW fast axis in the shear zone implies that subduction was considerably oblique to the margin, given the general NNE strike of the Appalachian orogen at this latitude. One potential complication, pointed out by van Staal and Zagorevski (2023), is that there are apparently no signs of modification of the apparent relict slab proposed by Luo et al. (2022) by Central Atlantic Magmatic Province (CAMP) volcanism (e.g., Marzoli et al., 2018). Our interpretation of the mantle lithospheric shear zone, which relies on this interpretation of the dipping slab-like structure, thus having some uncertainties and warrants further investigation. Signals in cos(baz) term inside this anisotropic layer suggest complicated variations of anisotropic geometry within Feature #7, which may indeed reflect modifications by younger processes such as CAMP volcanism.

Features #8 and #9 lie at asthenospheric mantle depths and represent two interfaces that can be observed beneath nearly the entire SEISConn array. Because of the dipping nature of Feature #9, there is a triangular "wedge" of the mantle sandwiched between these two interfaces (Figure 9). These two features manifest in different harmonic terms, which implies that they are not just the top and bottom interfaces of a single anisotropic layer; furthermore, there may be additional anisotropy above and/or below this wedge-shaped layer. Feature #9, with negative sin(2*baz) signals, could be generated by an anisotropic layer with horizontal NE-SW fast axis above the interface, while the region beneath the interface can be either isotropic or have a vertical axis of anisotropic symmetry; these possibilities cannot be distinguished by anisotropic RF analysis (or shear wave splitting analysis). Feature #8 has dominant signals in the cos(baz) term, so the region above Feature #8 is likely also anisotropic with a north-or south-plunging axis of symmetry.

One possible explanation for Feature #8 is that it corresponds to the lithosphere-asthenosphere-boundary (LAB). An 80-100 km deep LAB beneath southern New England is somewhat deeper than the estimates from several previous RF studies (e.g., Goldhagen et al., 2022;Hopper & Fischer, 2018) but more consistent with others (e.g., Abt et al., 2010;Rychert et al., 2007). The nature of the LAB, and how it manifests in seismic observations, is still imperfectly understood (e.g., Karato & Park, 2019). A transition from a plunging fast axis above a horizontal fast axis below effectively represents a decrease in radial anisotropy with depth, which can be a plausible origin of seismically observed LAB beneath continents (Rychert & Shearer, 2009). Therefore, we propose that Feature #8 represents the top of an asthenospheric flow zone, which may also correspond to the LAB beneath southern New England at 80-100 km depth. In this interpretation, the anisotropic aspects of Feature #8 reflect the contrast between anisotropy near the base of the lithospheric mantle with a plunging symmetry axis and anisotropy in the asthenosphere with a nearly horizontal symmetry axis parallel to the plate motion.

We further hypothesize that the plate motion-parallel asthenospheric flow beneath New England is modified by the presence of a vertical upwelling flow near the eastern end of the profile. In this interpretation, Feature #9 represents the sharp transition of the mantle flow direction from vertical (below the boundary) to horizontal (above the boundary), such that the plate motion parallel flow zone becomes thinner to the east. The sharp transition of the mantle flow direction suggested by the anisotropy contrast can be achieved by a small-scale or edge-driven convection cell (e.g., King & Anderson, 1998;King & Ritsema, 2000). A thinner zone of plate motion parallel flow would result in smaller path-integrated shear wave splitting due to olivine CPO in the zone, which explains the observed weaker SKS splitting beneath the eastern portion of the SEISConn array (Lopes et al., 2020). One possibility is that the proposed mantle upwelling is associated with the NAA centered beneath central New England, thought to be a region of buoyant present-day upwelling (e.g., Levin et al., 1995;Menke et al., 2016). Our inference of mantle upwelling associated with the NAA beneath southern New England does not necessarily negate a possible contribution from the Great Meteor hotspot (e.g., Eaton & Frederiksen, 2007;Tao et al., 2020) in the formation of the NAA.

The anisotropic RF analysis for the SEISConn array presented here reveals the presence of multiple dipping interfaces and anisotropic layers in the crust and upper mantle beneath southern New England. The harmonic decomposition technique separates RF signals according to their backazimuthal dependence, which can then be further explored either by Bayesian inversion, to obtain quantitative constraints on the seismic model, or by tracing signals with the same harmonic patterns across adjacent stations to resolve the shape and lateral extent of individual features. We are able to identify relatively shallow anisotropic and dipping features in the crust with various potential origins, including fine layering in the Hartford Basin (Hubert et al., 1992), east-dipping faults associated with the Taconic thrust belt (Stanley & Ratcliffe, 1985), as well as foliation of metamorphic rocks and frozen-in CPO of crustal minerals developed during Salinic or Acadian terrane accretion events (van Staal et al., 2009). A more elongated anisotropic interface at mid-crustal depths may record orogen-parallel crustal flow associated with the collapse of the Acadian orogenic plateau (Hillenbrand et al., 2022) or channel flow and ductile extrusion of the Putnam-Nashoba terrane (Severson, 2020). The backazimuth-dependent signals produced by the dipping Moho discontinuity are also well captured by the harmonic analysis. The extension of harmonic signals associated with the west-dipping shallow Moho across the Moho depth offset represents additional evidence for a previously inferred doubled Moho, supporting an overthrust model for the formation of the Moho depth offset (Luo et al., 2022). There is an anisotropic layer in the lithospheric mantle that may represent a shear zone above a proposed relict slab lodged in the lithosphere due to an episode of oblique, west-dipping subduction during the Alleghanian closure of the Rheic Ocean. This anisotropic layer might also reflect modifications by younger processes such as CAMP volcanism (Marzoli et al., 2018). In the asthenosphere, our results suggest plate motion parallel flow in the asthenosphere. This asthenospheric flow zone may be modified by mantle upwelling associated with the NAA beneath the eastern end of the array; this model can explain the previously observed eastward decrease of SKS splitting delay times (Lopes et al., 2020).

This study provides perhaps the most comprehensive and in-depth regional analysis of crustal anisotropy beneath eastern North America to date, which has generally been lacking compared to the extensive studies of mantle anisotropy (Fouch & Rondenay, 2006). Quantitative analysis of crustal anisotropy using Bayesian inversion constrains the detailed geometry of prominent anisotropic structures beneath a single station, while the more qualitative analysis based on the dense array data can help identify robust features with smaller amplitudes. The constraints that we obtained on layered mantle anisotropy complement, and provide further insights into previously obtained results from shear wave splitting. The depth resolution provided by RF analysis is particularly valuable in regions with complicated mantle flow patterns that cannot be fully resolved using splitting analyses, such as the upwelling flow we infer at the eastern end of the SEISConn line.

The workflow that we have developed for this study also has great potential to be applied to other dense arrays that cross different segments of the Appalachian orogen (or to other orogenic regions). The general alignment of fast splitting directions with the strike of the Appalachian Mountains in the central and southern Appalachians, reported in previous SKS splitting studies (Aragon et al., 2017;Long et al., 2016), suggests strong contributions from lithospheric deformation associated with the Appalachian orogenesis. Furthermore, Appalachian orogenesis is known for substantial along-strike variation, which is manifested in not only the surface geology (Hatcher, 2010;Hibbard & Karabinos, 2013) but also in the Moho geometry (Luo et al., 2023). Detailed information about layered lithospheric anisotropy along the Appalachians can further our understanding of deformation processes and their lateral variations during Appalachian orogenesis, and can potentially be extrapolated to shed light on other orogenic belts, including those that are active today. that helped improve the manuscript. Collection and analysis of SEISConn data were funded by Yale University and by the National Science Foundation via Grants EAR-1150722 and EAR-1800923. The authors thank the station hosts and the field personnel, including the participants in the Field Experiences for Science Teachers (FEST) program (Long, 2017), for the successful SEISConn experiment. The analysis and interpretation of the results benefit from the discussions with colleagues, including Jay Ague, Shunichiro Karato, Jeffrey Park, and the Yale seismology group. The authors also thank Jeffrey Park for providing the original versions of the MTC and harmonic decomposition codes.