Lithosphere Structure and Seismic Anisotropy Offshore Eastern North America: Implications for Continental Breakup and Ultra-Slow Spreading Dynamics

The breakup of supercontinent Pangea occurred ∼ 200 Ma forming the Eastern North American Margin (ENAM). Yet, the precise timing and mechanics of breakup and onset of seafloor spreading remain poorly constrained. We investigate the relict lithosphere offshore eastern North America using ambient-noise Rayleigh-wave phase velocity (12–32 s) and azimuthal anisotropy (17–32 s) at the ENAM Community Seismic Experiment (CSE). Incorporating previous constraints on crustal structure, we construct a shear velocity model for the crust and upper ∼ 60 km of the mantle beneath the ENAM-CSE. A low-velocity lid ( V S of 4.4–4.55 km/s) is revealed in the upper 15–20 km of the mantle that extends ∼ 200 km from the margin, terminating at the Blake Spur Magnetic Anomaly (BSMA). East of the BSMA, velocities are fast (>4.6 km/s) and characteristic of typical oceanic mantle lithosphere. We interpret the low-velocity lid as stretched continental mantle lithosphere embedded with up to ∼ 15% retained gabbro. This implies that the BSMA marks successful breakup and onset of seafloor spreading ∼ 170 Ma, consistent with ENAM-CSE active-source studies that argue for breakup ∼ 25 Myr later than previously thought. We observe margin-parallel Rayleigh-wave azimuthal anisotropy (2%–4% peak-to-peak) in ∼ 15% gabbro was crystallized within melt channels. Due to the slow spreading rate, corner flow at the ridge was weak relative to plate motion, resulting in olivine LPO oriented sub-parallel to the plate-motion direction in the underlying asthenosphere, which was subsequently frozen into the cooling oceanic mantle lithosphere. The combination of margin-parallel LPO in the oceanic lithosphere and significantly weaker shape-preferred orientation of the frozen gabbro channels in the lid likely contribute to the observed margin-parallel anisotropy.

and tectonic relationship between CAMP magmatism and rifting is debated (McHone, 2000), it is thought that normal seafloor spreading and opening of the Atlantic basin began sometime between ∼200 and 170 Ma.
Insights into the transition from continental rifting to seafloor spreading are contained in the crust and mantle signature offshore ENAM. The offshore region is characterized by two positive-polarity magnetic anomalies separated by the Inner Magnetic Quite Zone (IMQZ) that lacks well defined magnetic lineations ( Figure 1). The higher amplitude East Coast Magnetic Anomaly (ECMA) occurs just seaward of the continental shelf and has generally been interpreted as marking the transition to oceanic crust (e.g., Kelemen & Holbrook, 1995;Klitgord et al., 1988;Lynner & Porritt, 2017). ECMA emplacement ages range from 175 to 200 Ma (Benson, 2003;Klitgord & Schouten, 1986;Labails et al., 2010), though recent revised estimates of ∼195 Ma have been proposed based on the African conjugate to the ECMA as well as salt deposits off Nova Scotia and Morocco (Labails et al., 2010;Sahabi et al., 2004). Approximately 200 km seaward of the ECMA is the lower amplitude Blake Spur Magnetic Anomaly (BSMA). The age of the BSMA is estimated at ∼170 Ma, but its precise origin and significance are still debated (Greene et al., 2017).
Recent work offshore ENAM has challenged the notion that the ECMA marks the completion of continental breakup and onset of seafloor spreading, and instead, it has been proposed that the BSMA marks this important transition Shuck et al., 2019). Using data collected during The ENAM Community Seismic Experiment (ENAM-CSE) that extends farther offshore than ever before, detailed crustal imaging shows thin proto-oceanic crust with higher lower-crust velocities and rougher basement topography west of the BSMA compared to the east, consistent with deeper melting and ultra-slow opening rates (∼0.65 cm/yr half-rate) Shuck et al., 2019). This interpretation predicts a 15-20 km thick continental mantle lithosphere underlying the oceanic crust west of the BSMA and normal oceanic lithosphere east of the BSMA; however, previous shear-velocity imaging in the region shows little evidence for a distinct change in mantle velocities across the BSMA (Lynner & Porritt, 2017).

of 30
Additional insights into dynamics associated with continental rifting and seafloor spreading may be preserved in the olivine lattice-preferred orientation (LPO) in the lithospheric mantle, which acts as a record of past mantle flow. The fast [100] axis of olivine crystals tends to align parallel to the direction of shearing in the mantle, forming an LPO (Karato et al., 2008;Zhang & Karato, 1995). In models of mid-ocean ridges (MOR), corner flow near the ridge aligns olivine LPO parallel to the spreading direction and is frozen-in as the lithosphere cools (e.g., Blackman & Kendall, 2002;Blackman et al., 1996;Kaminski & Ribe, 2001Ribe, 1989). This frozen-in LPO leads to the azimuthal anisotropy of seismic waves routinely observed in the Pacific lithosphere, with a fast azimuth parallel to the fossil-spreading direction (FSD) (e.g., Forsyth, 1975;Hess, 1964;Mark et al., 2019;Morris et al., 1969;Raitt et al., 1969;Russell et al., 2019). Seismic anisotropy of the deeper asthenosphere reflects present-day mantle deformation and at large length scales broadly aligns sub-parallel to absolute plate motion (APM) beneath the ocean basins (Beghein et al., 2014;Nishimura & Forsyth, 1989;Schaeffer et al., 2016), with deviations associated with smaller-scale convective processes (e.g., Becker et al., 2014;P.-Y. P. Lin et al., 2016).
Previous seismic anisotropy observations at the ENAM-CSE from shear-wave splitting show margin-parallel fast axes, significantly rotated from the FSD and current-day APM and interpreted as present-day asthenospheric flow along the margin (Lynner & Bodmer, 2017). In addition, preliminary comparison of sub-Moho V P along crossing margin-parallel and margin-perpendicular refraction lines in the region suggests a margin-parallel fast direction in the shallow lithosphere that is approximately perpendicular to the FSD (Shuck & Van Avendonk, 2016). However, these observations are limited to two locations in the ENAM-CSE footprint where the refraction lines intersect.
In this study, we use ambient-noise Rayleigh waves to construct a shear velocity model of the offshore ENAM-CSE region that incorporates recent crustal constraints from refraction tomography (Shuck et al., 2019). Our model reveals relatively low-velocity mantle lithosphere extending ∼200 km seaward that we interpret in the context of the detailed crustal architecture, providing further evidence for a prolonged breakup prior to seafloor spreading. We also report margin-parallel Rayleigh-wave anisotropy in the lithosphere, perpendicular to typical expectations for seafloor spreading, and offer an alternative perspective on lithosphere fabric formed at slow-spreading ridges.

Data and Methods
The Eastern North American Margin Community Seismic Experiment (ENAM-CSE) consisted of onshore-offshore active-source reflection and refraction as well as a one year broadband ocean-bottom seismometer (OBS) deployment (Lynner et al., 2020). We use continuous seismic data from 28 broadband OBS deployed during the ENAM-CSE from April 2014 to May 2015. Water depth in the study region ranges from ∼1300 m near the shelf to ∼5200 m on the eastern-most edge of the array. Instrument response is deconvolved to displacement, and seismograms are downsampled to 1 Hz prior to processing.

Daily OBS Tilt Noise Removal
We observe exceptionally strong horizontal noise on the vertical channels (i.e., tilt noise) at ENAM at periods 10 s (Figure 2a), presumably due in part to the strong Gulf Stream current that flows northeastward along the coast. We remove this coherent horizontal energy from the vertical channels for each 24 hr segment of the continuous data using the Automated Tilt and Compliance Removal (ATaCR) software (Janiszewski et al., 2019), which implements the techniques developed by Crawford and Webb (2000). We do not remove pressure coherence from the vertical channel, as this has been shown to degrade the desired fundamental-mode primary microseism (Bowden et al., 2016).
It is not common practice to remove daily tilt noise prior to performing ambient-noise tomography at OBS experiments, but we find that tilt removal improves the overall signal-to-noise ratio (SNR) of vertical component empirical Green's functions (EGFs) (See Section 2.2) by a factor of ∼2 on average, and by an order of magnitude for some station pairs (Figure 2b). The largest SNR improvements occur for station pairs with shallower average water depth (Figure 2d). Similar ambient noise SNR improvements were reported after tilt and compliance corrections at the shallow water stations at the Cascadia Initiative (Tian & Ritzwoller, 2017

Ambient Noise Processing
Ambient noise EGFs are constructed from tilt-removed seismograms following the general procedure of Bensen et al. (2007) (Figure 3); however, we do not apply time-domain normalization. Daily displacement seismograms are split into 15 3-hr segments with 50% overlap. The deployment average coherence cross-spectrum ρ ij is calculated between the vertical channels of two stations i and j as follows (Ekström, 2014): Figure 2. Ambient noise empirical Green's function (EGF) signal-to-noise ratio (SNR) improvement from daily tilt corrections. (a) Example 24 hr smoothed spectra for station X06 on day 47 of the deployment. The raw vertical channel (Z) and tilt-corrected vertical (Z-1-2) are shown. The gray region is bounded below by the Peterson low noise model (PLNM) and above by the high noise model (PHNM) (Peterson, 1993

of 30
where k represents a single time window and N is the total number of windows for the deployment. U(ω) is the vertical component displacement spectra at angular frequency ω, and U*(ω) is its complex conjugate. By utilizing relatively short time windows, any windows containing anomalous signals such as large earthquakes have little influence on the final stacked spectrum, precluding the need for time-domain normalization or other alteration of the original high-quality waveforms. We found that a typical one-bit normalization procedure (e.g., Bensen et al., 2007) degrades signal-to-noise ratio in the 12-35 s period band by more than an order of magnitude ( Figure S1 in Supporting Information S1).

Interstation Phase Velocities
Interstation phase velocities are estimated from the stacked coherence spectra using Aki's spectral formulation, whereby the real part of the cross-spectra takes the form (Aki, 1957): where r is station separation, c is phase velocity, and J 0 is the Bessel function of order zero. Phase velocity dispersion is estimated at each zero crossing following Ekström et al. (2009) and then interpolated to a uniform frequency axis (Figure 4). This process identifies multiple dispersion curves that all satisfy the data, and therefore, we select the one with velocities from 25 to 30 s that are closest to a nominal mantle velocity of 3.9 km/s (Figure 4c), based on the average of our starting shear velocity model predictions at those periods ( Figure 5c). We then discard dispersion curves that are not smooth or do not decrease with increasing frequency. In order to minimize noise in the cross-spectrum prior to the zero-crossing analysis, we apply a cosine-taper window in the time domain defined by a minimum group velocity threshold of 0.3 km/s (that is, energy corresponding to group velocities 0.3 km/s is zeroed). Additional smoothing is applied to for periods 17 s using a moving average with a window length of 1.85 mHz to eliminate the occurrence of spurious zero crossings, particularly for large r. Dispersion measurements for r less than one wavelength are discarded.
Equation 2 is valid only for an evenly distributed noise source (Aki, 1957). It has been shown that an azimuthal bias in noise sources introduces higher order 2θ and 4θ azimuthal terms to Equation 2 that can map into an apparent azimuthal anisotropy of phase velocity (Harmon et al., 2010). This is an important point to consider given that we are solving for azimuthal anisotropy near a coastline, where the noise source could be biased. The EGFs in Figure 3 appear quite symmetric, suggesting an even distribution of noise sources, but we more explicitly evaluate noise source distribution by plotting the real and imaginary part of the cross-spectrum as a function of distance and azimuth for several frequencies ( Figure S2 in Supporting Information S1), following Harmon et al. (2010). We find strong, coherent energy in the real component that is exceptionally uniform with respect to azimuth ( Figure S2a in Supporting Information S1) and nearly no coherent structure in the imaginary component ( Figure  S2b in Supporting Information S1). Taken together, this indicates a uniform distribution of noise sources and validates our use of Equation 2 for extracting interstation phase velocities.

Phase Velocity Inversion
Interstation Rayleigh-wave phase velocities are inverted for 2-D phase velocity maps from 12 to 32 s period and azimuthal anisotropy from 17 to 32 s for the offshore ENAM region. We solve first for phase slowness maps s(x) on a 0.1° × 0.1° grid and take the reciprocal to obtain maps of phase velocity c(x). Anisotropy terms are solved on a coarser 0.5° × 0.5° grid. A perturbation in phase delay time δτ ij between stations i and j due to a perturbation in phase slowness δs(x) is given by where θ ij is propagation azimuth and r ij is the great-circle distance between the stations. Coefficients A c and A s describe the frequency-dependent azimuthal anisotropy within the array footprint with magnitude = √ 2 + 2 and fast azimuth ψ = 0.5 tan −1 (A s /A c ). The 2-D finite-frequency sensitivity kernel K ij = ∂τ/∂s is given by (F. C. Lin & Ritzwoller, 2010) where τ i (x) is the phase delay field due to an impulse at station i, † ( ) is the adjoint phase delay field due to an impulse at station j, and τ i (x j ) is the value of the phase delay field at station j (i.e., the minimum travel time from station i to station j). Reference velocity c 0 is taken as the average of all interstation dispersion measurements at frequency ω. To ensure that the kernel density matches the ray-theoretic value, the kernel is normalized such that This formulation of the kernel, termed the "empirical" kernel by F. C. Lin and Ritzwoller (2010), accounts for both finite-frequency effects and off-great-circle propagation caused by lateral velocity gradients along the ray path. In theory, the quantities τ i (x) and † ( ) can be determined empirically for each station by fitting a smooth surface to the interstation phase delays measured at all other stations across the array. In practice, however, this is challenging due to uneven data distribution and presence of noise in the phase delay measurements.
We take a "semi-empirical" approach and approximate τ(x) numerically via spectral-element method (SEM) simulations by propagating S-waves through a realistic phase velocity map at each frequency of interest for an impulse centered at each station. We estimate τ(x) by cross-correlating waveforms at each point in the model with the reference waveform at the source, τ(0). We then plug these numerical estimates of the phase delay field into Equation 4. The input synthetic phase velocity maps are constructed from the 2-D V P model along Line 1 from Shuck et al. (2019) after converting to V S assuming V P /V S = 1.8 for the crust and mantle and 2.22 for the sediments ( Figure 5) and extrapolating phase velocities to the full ENAM footprint via constant depth contours. These estimates of V P /V S are taken from the V P and V S models along line three of Shuck et al. (2019). Example finite-frequency kernels are shown in Figure 6 for 13 and 20 s period. Equation 4 gives the sensitivity kernel at an instantaneous frequency ω and contains all Fresnel zones, but in practice each phase velocity measurement is made over a finite bandwidth [ω -Δω/2, ω + Δω/2]. We approximate finite-bandwidth kernels by forming a Gaussian weighted average of instantaneous kernels centered on ω with a half width of 10%, effectively limiting the spatial extent of the kernel to the first few Fresnel zones.
Phase slowness maps and anisotropy are inverted via Equation 3 using a linearized iterative least squares approach (Menke, 2012), and the final model is obtained after two iterations. Perturbations to the starting homogeneous, isotropic phase slowness model is regularized using wavelength-weighted second derivative smoothing and norm damping toward the starting model. At periods less than 17 s, anisotropy terms are damped to zero, as variations in phases velocity associated with anisotropy are swamped by the large variations associated with water depth. Sensitivity kernels are updated upon the second iteration to account for off-great-circle propagation associated with lateral gradients in slowness.

Off-Great-Circle Propagation
Water depth across ENAM increases by nearly 4000 m from west to east, leading to a drastic difference in short period Rayleigh-wave sensitivity across the array that is strongly controlled by water depth (Figure 7). A 13-s Rayleigh wave that primarily samples the eastern-most edge of the array where water is deep (4000-5000 m) will travel slower than the same wave sampling shallow water near the shelf. This strong lateral velocity gradient at short periods steers energy toward the faster shallow-water regions leading to off-great-circle propagation (Figure 6). If great-circle paths are assumed, as is routinely done in surface-wave tomography studies, an apparent 2θ azimuthal anisotropy signal is produced, where waves that travel parallel to the margin appear faster on average than waves that travel perpendicular to the margin ( Figure 6e). Left unaccounted for, this off-great-circle propagation can result in (1) apparent 2θ azimuthal anisotropy with a fast axis parallel to the margin and (2) biased fast isotropic velocities, on average. We avoid these biases by using semi-empirical finite-frequency kernels that account for off-great-circle propagation and by limiting anisotropy measurements to periods greater than 17 s, where the anisotropy bias is small (Figure 6e).

Inversion for Shear Velocity, V S
The single-frequency phase velocity maps are resampled to a 0.25° × 0.25° grid and combined to produce a dispersion curve at each grid point in the model. Each dispersion curve is individually inverted for a 1-D depth dependent shear velocity (V S ) profile, and all profiles are ultimately combined to produce the final 3-D V S model. The starting 3-D V S model for the inversion is the same as that used for the semi-emperical kernels in Section 2.4, built up from the 2-D crustal V P model along Line 1 from Shuck et al. (2019) assuming V P /V S = 1.8 in the crust and mantle and 2.22 in the sediments. These estimates of V P /V S are based on averages from V P and V S models along Line 3 from Shuck et al. (2019). The mantle is not well resolved by the P n modeling, so we Figure 6. Finite-frequency sensitivity kernels and importance of off-great-circle propagation. (a) Synthetic phase velocity map at 13 s period after Gaussian smoothing (σ = 0.2°). Gray contours show S-wave phase delay at 25 s intervals due to an impulse at station X10 modeled using the spectral element method (SEM) (i.e., τ(x) in Equation 4). Note the strong wavefield distortion due to large lateral velocity gradients associated with changes in water depth. (b) Finite-frequency sensitivity kernel for station pair X10-X04. The dashed and solid black lines show the great circle (GC) and SEM ray paths, respectively. (c),(d) same as (a), (b) but for 20 s period. The SEM ray path is nearly coincident with the GC. (e) Apparent azimuthal anisotropy of GC phase delays (τ GC ) relative to SEM phase delays (τ SEM ) for all station pairs with r > 100 km at 13 s (red) and 20 s (blue). assume a nominal mantle V P of 8.0 km/s. As seismic structure mostly varies perpendicular to the margin (Shuck et al., 2019) (i.e., with water depth), the 2-D line is extrapolated to the entire 3-D ENAM region along contours of constant water depth.
Each 1-D inversion aims to solve for V S while minimizing dispersion misfit via a standard linearized iterative least-squares approach (Menke, 2012). The perturbation kernels and forward estimates of phase velocity are calculated for a layered Earth model extending to 250 km using the SURF96 software (Herrmann, 2013) (Figure 7). The inversion is regularized with norm damping toward the starting model, second derivative smoothing, and a constraint that seeks to preserve layer gradients in the crust and mantle (Russell et al., 2019). These constraints are broken across prescribed sediment and Moho boundaries to allow for discontinuities. Compressional velocities (V P ) and density are held fixed during a single inversion, but their values are allowed to be perturbed within the Monte Carlo framework described below. A linear increase in damping toward the reference model is applied below 50 km such that the reference model is exactly maintained by ∼100 km depth. A total of 10 iterations are performed: iterations 1-9 maintain the same sensitivity kernels, and dispersion is predicted via perturbation theory using the kernels; upon iteration 10, the kernels and predicted dispersion are recalculated using SURF96.
In order to evaluate starting model dependence, a Monte Carlo approach is used. At each grid point, the starting model sediment, crust, and mantle velocities and layer thicknesses are randomly perturbed by drawing from a zero-mean normal distribution with standard deviation of 10% of the reference values for the sediment, crust, and mantle. This process is repeated 100 times at each grid point. Next, we repeat the entire inversion processes for a 3-D starting model constructed from Line 2, which is parallel to Line 1 but located to the south (Shuck et al., 2019). The median model resulting from the Line 2 inversion is similar to the Line 1 inversion to within uncertainty. We combine the two inversion results into an ensemble of 200 models from which the median and middle 68% of models with χ 2 ≤ 1.75 are taken as the preferred V S and associated 1-σ uncertainties. These model uncertainties capture both uncertainties in the dispersion data as well as uncertainty associated with the starting model and our choice to keep V P and density fixed within each inversion.  Figure 8 shows isotropic phase velocity maps for periods ranging from 13 to 32 s. The dashed white boundary shows the resolving limit defined by a value of 0.1 on the diagonals of the frequency-averaged resolution matrix. At 13 s period (Figure 8a), Rayleigh wave velocities decrease from west to east due to the increasing water depth. At longer periods this general trend reverses, and velocities increase eastward as a result of the transition from thick continental crust in the west to thinner oceanic crust (and increased mantle sensitivity) in the east. This same character is seen in the synthetic phase velocity maps in Figures 6a and 6c (see also Figure S7 in Supporting Information S1). Similar to the synthetic maps, variations parallel to the margin are smooth and relatively minor, suggesting that the measurements are largely 2-D and primarily controlled by structure associated with the continent-to-ocean transition. At the longest periods (27-32 s), the highest velocities (>4.1 km/s) are observed east of the BSMA, although resolution is limited by modest station distribution east of the BSMA.

Azimuthal Anisotropy
Azimuthal anisotropy of phase velocity is estimated from 17 to 32 s, sensitive primarily to mantle lithosphere depths. Figure 9 shows 2-D maps of azimuthal anisotropy for two different choices of smoothing: (1) A smoothly parameterized inversion that seeks to minimize the second spatial derivatives ∇ 2 A c and ∇ 2 A s . (2) An inversion that effectively solves for a single anisotropy fast azimuth and magnitude on either side of the BSMA by strictly enforcing ∇A c = ∇A s = 0 but breaking this constraint at the BSMA.
Both inversions indicate dominantly margin-parallel anisotropy across the array. In detail, the smooth model (Figure 9a) shows a gradual west-to-east clockwise rotation from margin-parallel. East of the BSMA and for 17-25 s period, both models indicate a ∼25-45° rotation from margin-parallel toward the FSD (Figure 10). Neither inversion suggests FSD-parallel anisotropy anywhere in the model. At periods 27 s, both models indicate a reduction in anisotropy east of the BSMA. The smooth model shows evidence of a southeastward reduction in anisotropy magnitude west of the BSMA, ranging from ∼3% to 4% in the northwest to ∼1%-2% in the southeast; however, these magnitude variations may not be well resolved (see Section 3.2.1; Figure S4a in Supporting Information S1) As the 2-D anisotropy maps suggest a relatively simple pattern dominated by margin-parallel anisotropy, we also solve for a 1-D array-averaged fast azimuth and magnitude for the region. Phase velocity residuals calculated relative to the 2-D path-integrated isotropic velocities are shown as a function of azimuth in Figures 11a-11f. Azimuthal patterns are clearly dominated by a 2θ sinusoid with fast azimuth parallel to the margin and perpendicular to the FSD (Figures 11g and 11h). Anisotropy peak-to-peak magnitudes decrease slightly with increasing period from 3%-4% at 17 s to 2%-2.5% at 32 s. To address potential bias caused by uneven azimuthal sampling, we also estimate anisotropy strength and fast azimuth for data averaged in 20° azimuthal bins. The fast azimuths and magnitudes obtained agree with the unbinned estimates but have uncertainties of a factor ∼2-3 larger.
We explore whether the simpler 1-D anisotropy model can explain the data as well as the 2-D models and find that it can in most cases ( Figure 12). Phase delay residuals are compared for the three anisotropy inversions with varying degrees of freedom-the 1-D array-averaged inversion, 2-D smooth inversion, and 2-D break at the BSMA. Overall, the 2-D smooth inversion yields slightly lower root-mean-square (RMS) phase delay residuals on average, which is unsurprising given that it has the most degrees of freedom. However, for most frequencies the data fits produced by each of the three inversions are nearly indistinguishable (Figure 12g), and therefore we favor the simplest model with 1-D margin-parallel anisotropy, though we cannot rule out a small (25-45°) clockwise rotation east of the BSMA.

Synthetic Recovery Tests: Resolving Anisotropy East of the BSMA
We perform synthetic tests to determine whether our dataset can resolve a 90° change in anisotropy fast direction from margin-parallel west of the BSMA to FSD-parallel east of the BSMA ( Figure S3-S6 in Supporting Information S1). At each frequency, the synthetic dataset is calculated using the synthetic isotropic phase velocity maps extrapolated from Line 1 and the anisotropic model shown in Figure S3 in Supporting Information S1. The input anisotropy model is characterized by a constant 2% zero-to-peak anisotropy magnitude across the array with fast azimuthal parallel to the margin west of the BSMA and parallel to the FSD east of the BSMA. Gaussian noise with standard deviation of 0.6 s is added to the synthetic phase delay times. We test the two different 2-D inversion strategies introduced above to evaluate their ability to recover the change in anisotropy direction across the BSMA.
We find that the smooth anisotropy inversion ( Figure S4 in Supporting Information S1) successfully resolves the margin-parallel anisotropy west of the BSMA and FSD-parallel anisotropy east of the BSMA but has trouble capturing the sharp change directly at the BSMA, where recovered anisotropy magnitudes are small. Anisotropy strength is generally poorly resolved in the smooth model; it is overestimated at the westernmost edge of the array and underestimated east of the BSMA. The inversion containing the break at the BSMA ( Figure S5 in Supporting Information S1) more accurately captures both the anisotropy magnitude and change in direction across the BSMA. This is expected given that this parameterization closely resembles the character of the input anisotropic model.
Overall, both parameterizations are able to capture a 90° change in anisotropy across the BSMA, but the parameterization that explicitly allows a break at the BSMA more accurately recovers the magnitude of anisotropy ( Figure S6 in Supporting Information S1). In both cases, the array-averaged anisotropy is dominated by the margin-parallel anisotropy structure west of the BSMA, due to its larger footprint and the greater number of ray paths on the western side of the array. We conclude that if a change in anisotropy direction occurs from margin-parallel west of the BSMA to FSD-parallel east of the BSMA, the 2-D anisotropy inversions should recover it. It follows that we can confidently state that FSD-parallel fabric is not present within the ENAM-CSE region. Figure 11. (a)-(f) Array-averaged 1-D azimuthal anisotropy ranging from 17 to 32 s period. Gray circles show interstation phase velocity deviations from the pathaveraged isotropic values, and black diamonds show the mean and standard deviation for data within 20° azimuthal bins. The 2θ fit to the individual measurements is shown in red, and the fit to the binned data is in black. The margin-parallel direction and fossil-spreading direction (FSD) are indicated by the vertical blue and brown lines, respectively. (g) Anisotropy peak-to-peak strength and (h) fast propagation azimuth for the unbinned data in red and binned data in black. Error bars show 2-σ uncertainties.

V S Uncertainty
Example 1-D V S profiles from the west and east are shown in Figure 13 with uncertainties estimated via the Monte Carlo inversion described in Section 2.5. Uncertainties are larger in the crust and sediments than in the mantle, particularly for the western profile (Figures 13a and 13b) where the rapid change in water depth makes resolving crustal structure more difficult. In the western profile, crustal uncertainties generally range from 0.5 to 0.55 km/s compared to 0.25-0.3 km/s in the east. Mantle velocities for both the western and eastern profiles appear well resolved with uncertainties ranging from 0.1 to 0.3 km/s in the upper 10 km of the mantle and 0.03-0.1 km/s below (see also Figure S9 in Supporting Information S1). The relatively large uncertainties in the lower crust and uppermost mantle along the western edge of the array indicate an unresolved tradeoff (Figures S9b and S9bc in Supporting Information S1). For this reason, we do not interpret structure in that part of the model.

Crustal Variations
Horizontal slices through the 3-D V S model are shown in Figure 14. In the upper 4 km of the crust, velocities are slow (2-2.5 km/s) in the west corresponding to sediments and transition sharply to faster velocities (3-3.6 km/s) as the sediments thin to 4 km in the east (Shuck et al., 2019). The sharpness of this basement transition is a product of the layer discontinuities built into the starting reference model ( Figure S8 in Supporting Information S1), which are preserved but otherwise poorly resolved by Rayleigh waves. In the lower crust (8 km), a slow anomaly is observed just east of the ECMA.  In the 4-15 km depth slices, a band of anomalously slow velocities coincides with the shelf, correlating with large gradients in seafloor depth (>2 km/°) and large uncertainties. This rapid decrease in water depth west of the shelf break is just outside the array and poorly resolved by our relatively smooth phase velocity maps. This results in slower than predicted phase velocities along the shelf that map directly into the crustal velocities. The region of the model most affected by this bias is indicated by a semi-transparent mask in Figure 14, and we do not interpret structure within that part of the model.

Mantle Variations
Velocities in the mantle are generally well resolved in comparison to the crust as the mantle-sensitive Rayleigh waves are not affected by changes in water depth. Within the array footprint, mantle velocities range from ∼4.4 km/s in the west to ∼4.7 km/s in the east (Figures 14d and 14e). East of the BSMA and within the array footprint, velocities peak at 4.6-4.7 km/s in the upper 10-25 km of the mantle and decrease with depth, suggestive of thermal ocean lithosphere. At 55-60 km depth, velocities are nearly homogeneous (∼4.5 km/s) across the region, likely due in part to our choice of a linear increase in damping toward the starting model below 50 km depth.

2-D Margin-Perpendicular Structure
The average 2-D margin-perpendicular structure of our model shows more clearly the contrast in shallow mantle velocities that coincides approximately with the BSMA (Figure 15). In the upper 15-20 km of the mantle, relatively slow shear velocities (4.4-4.55 km/s) extend east of the margin ∼200 km terminating at the BSMA. East of the BSMA, velocities are elevated (>4.6 km/s) and characteristic of depleted, cold oceanic lithosphere (e.g., Nishimura & Forsyth, 1989  across the region, and therefore, we interpret that the low-velocity lid within the IMQZ is underlain by mantle lithosphere with temperature and composition consistent with the adjacent oceanic lithosphere. The average margin-perpendicular V S structure in our model is broadly similar to the offshore structure in the previous shear velocity model of Lynner and Porritt (2017). However, upon closer inspection important differences do exist. Their inferred offshore crust is ∼10 km thicker than what we observe. They also do not observe low velocities in the uppermost mantle west of the BSMA, nor do they observe a change in mantle velocities across the BSMA. These differences are likely attributed primarily to our differing starting model assumptions, and in particular, our direct incorporation of the crustal constraints from Shuck et al. (2019). In addition, our phase velocity maps explicitly account for azimuthal anisotropy and therefore should more accurately capture the true variations in isotropic phase velocity. Finally, that our modeling is focused exclusively on the offshore region rather than the full margin-crossing ENAM region both simplifies the modeling and avoids potential smearing of continental structure into the oceanic domain.

Low Velocity Mantle Lid West of the BSMA
We explore two hypotheses for explaining the 15-20 km thick low velocity mantle lid west of the BSMA. First, velocities of the lid are similar to continental mantle velocities previously imaged just onshore (Shen & Ritzwoller, 2016) (Figure 16a) and therefore may correspond to stretched and/or thermo-chemically modified continental lithosphere associated with rifting (e.g., Hopper et al., 2020). This interpretation of the mantle structure is consistent with the detailed crustal structure recently imaged at the ENAM-CSE Shuck et al., 2019). Shuck et al. (2019) observed relatively thin crust (6-8 km) with fast lower crustal velocities (V P ∼ 7.5 km/s) on Lines 1 and 2 between the ECMA and BSMA that is best explained by deeper than usual mantle melting. This led them to predict the presence of a 15-20 km thick continental lithosphere that truncated the upper melting regime producing less voluminous, more mafic melts. The slower shear velocities (V S < 4.55 km/s) that we observe in the upper mantle west of the BSMA are consistent with this prediction.
Second, the slow lithospheric velocities could be explained by gabbroic inclusions trapped within oceanic mantle lithosphere during the initiation of ultra-slow spreading (e.g., Gaherty & Dunn, 2007;Lizarralde et al., 2004). Efficient conductive cooling at the MOR during periods of very slow spreading leads to a thicker mantle lid that could crystallize gabbroic melts and trap them in the mantle, reducing shear velocity and producing thinner crust overall. This process may be particularly present in systems where the melt extraction process is inefficient relative to the rate of melt production, as could have been the case during the earliest stages of breakup. West of the BSMA, basement roughness is characteristic of typical ultra-slow spreading crust with a half-spreading rate of ∼0.65 cm/yr . This offers an alternative mechanism for truncating the upper melting regime without invoking the emplacement of proto-oceanic crust on top of continental lithosphere. However, one possible caveat is that the notion of a thicker thermal lid may be difficult to reconcile with the relatively hot mantle potential temperatures of 1395-1420°C that have been inferred based on the crustal structure in the region (Shuck et al., 2019), but detailed modeling of thermal evolution of ultra-slow spreading centers in the presence of elevated mantle temperatures is needed to evaluate temperature tradeoffs. Lithospheric lids with normal (V S > 4.5 km/s) upper-mantle velocities are observed along several stretches of the slow-spreading Mid-Atlantic Ridge (Gaherty & Dunn, 2007;Grevemeyer, 2020;Harmon et al., 2020;Ryberg et al., 2017;Saikia et al., 2021). This suggests that only ultra-slow spreading environments and/or spreading systems with relatively inefficient melt extraction, such as during the earliest stages of seafloor spreading, are conducive to trapping gabbro in the mantle.
We evaluate whether the addition of gabbro to nominal oceanic mantle (second hypothesis) can explain the slow V S west of the BSMA (Figure 16). We start by defining a reference oceanic V S profile as the average east of the BSMA. A nominal gabbro V P profile is estimated following Carlson and Miller (2004): where gabbro 0 = 6.85 km/s is the reference velocity of gabbro at the Moho (z moho ) and the change in V P due to pressure and temperature effects at depth are (∂V P /∂z)| T = 0.02 s −1 and (∂V P /∂z)| P = − 0.01 s −1 . This profile is converted to V S assuming V P /V S = 1.8. Next, using the Voigt average we solve for the proportion of gabbro required to reduce velocities of the oceanic reference model to match velocities observed west of the BSMA. The resulting depth distribution of gabbro in Figure 16b shows that up to ∼35% gabbro by volume would be required to explain velocities west of the BSMA for an oceanic reference model. Integrating this gabbro distribution implies crustal thinning west of the BSMA of ∼4 km (assuming crust that is pure gabbro), which is more thinning than observed (0.7-1.8 km) in the reflection and refraction imaging Shuck et al., 2019) (Figure 16c). In contrast, for a continental reference profile calculated from the onshore structure of Shen and Ritzwoller (2016), only up to ∼15% gabbro and ∼1 km of thinning are required, which is consistent with the crustal observation. Though simple, this back-of-the-envelope calculation suggests a combination of the two proposed hypotheses: slow mantle velocities west of the BSMA are due to extended continental lithosphere that contains up to ∼15% gabbro retained in the mantle during the volcanic process that built the overlying basaltic crust.

Implications for the Onset of Normal Seafloor Spreading
The 15-20 km thick slow lithospheric lid west of the BSMA that we attribute to continental mantle has implications for the timing and mechanics of breakup and onset of normal seafloor spreading in the North Atlantic. In keeping with Shuck et al. (2019) and Bécel et al. (2020), our results suggest that complete breakup and the onset  Bécel et al. (2020). Crustal thinning is consistent with up to ∼15% gabbro added to continental mantle but is inconsistent with oceanic mantle. of normal seafloor spreading occurred at ∼170 Ma, approximately 10-25 Myr after the initiation of rifting and formation of the ECMA. Crust that was formed following emplacement of the BSMA is consistent with normal seafloor spreading conditions (Shuck et al., 2019), and mantle structure is characteristic of typical oceanic lithosphere (Figure 16a). The BSMA also marks an increase in spreading rate from ultra slow (∼0.65 cm/yr half-rate) to slow (∼1.3 cm/yr half-rate) as inferred by the decrease in basement roughness stepping eastward across the BSMA .
Previously, complete breakup and the onset of normal seafloor spreading was thought to have taken place over a relatively short-lived period that immediately followed formation of the ECMA (Holbrook et al., 1994;Holbrook & Kelemen, 1993;Kelemen & Holbrook, 1995). Rather, the detailed lithospheric architecture at the ENAM-CSE indicates that after formation of the ECMA and prior to the BSMA, a proto-seafloor spreading mode was active with mantle-derived melts migrating vertically through the continental lid forming thin proto-oceanic crust Shuck et al., 2019). Mantle velocities suggest that up to ∼15% gabbro was crystallized and retained in the lid during melt migration, likely in channelized structures (e.g., Katz & Weatherley, 2012). The potential for melts to remain trapped in the mantle has implications for our understanding of melt extraction, crustal accretion, and thermal evolution at ultra-slow spreading centers.
Numerical models have shown that mantle can remain intact even after the crust has separated at rifted margins if the crust fails through brittle faulting while the underlying mantle lithosphere deforms through ductile necking (Huismans & Beaumont, 2011). Such ductile deformation of the lithosphere is conceivable considering the elevated mantle potential temperatures (1430-1480°C) inferred from CAMP lavas (Callegaro et al., 2013) and the inferred ultra-slow extension rates (<2 cm/yr half-rate) Davis et al., 2018).

Comparison to Previous Anisotropy Observations
The Rayleigh-wave anisotropy that we observe in the lithosphere at the ENAM-CSE has a fast direction sub-parallel to the margin. Previous shear-wave splitting observations at the ENAM-CSE also showed margin-parallel anisotropy, but it was attributed to deeper present-day asthenospheric flow along the margin associated either large-scale density or pressure driven flow or edge-driven convection due to the large lithospheric root of the continent (Lynner & Bodmer, 2017). Our observations do not preclude such flow in the asthenosphere but show that strong margin-parallel anisotropy is present in the lithosphere, potentially contributing at least in part to the previous shear-wave splitting observations. A preliminary comparison of shallow mantle V P along the margin-parallel and -perpendicular refraction lines at ENAM also indicated (∼8%) faster velocities parallel to the margin (Shuck & Van Avendonk, 2016) in agreement with our observations, though spatially limited to the crossing points of the profiles.
At the Far-offset Active-source Imaging of the Mantle (FAIM) seismic refraction experiment (115-130 Ma) ∼800 km southeast of the ENAM-CSE, Gaherty et al. (2004) inferred a lithospheric olivine LPO sub-parallel to the FSD (to within ∼15°). This implies that a rotation in lithosphere LPO from spreading-perpendicular to spreading-parallel occurred between 165 and 130 Ma. Observations of Love-Rayleigh scattering offshore ENAM offers evidence for such lateral gradients in anisotropy (Servali et al., 2020), though the frequencies considered (∼100 s) have significant asthenospheric sensitivity, and the inferred scattering points are widely distributed throughout the western Atlantic including within the ENAM-CSE footprint.
Recent surface-wave imaging at the Passive Imaging of the Lithosphere Asthenosphere Boundary (PI-LAB) experiment shows complex patterns of lithospheric anisotropy near the Mid-Atlantic Ridge (Saikia et al., 2021). Generally, at the ridge axis anisotropy is sub-parallel to the ridge and off-axis rotates to the FSD, but there is significant lateral variability. The Mid-Atlantic Ridge INtegrated Experiments at Rainbow (MARINER) experiment found anisotropy that is perpendicular to spreading in the upper 4 km of the lithosphere, interpreted as aligned cracks due to extensional stresses during rifting (Dunn et al., 2017). In addition, global models show variable anisotropy at lithospheric depths in the North Atlantic that is often highly rotated from the FSD (e.g., Becker et al., 2014;Debayle & Ricard, 2013;Schaeffer et al., 2016;Yuan & Beghein, 2013). As structure in the Atlantic varies significantly between different global models, our high-resolution estimate provides a useful benchmark for the northwest Atlantic and supports the notion of variable anisotropy in the Atlantic lithosphere.

Plate-Motion Controlled LPO
The surface-wave anisotropy that we observe is similar in strength and overall character to that observed in other oceanic regions (e.g., Forsyth, 1975;Nishimura & Forsyth, 1989;Takeo et al., 2018;P.-Y. P. Lin et al., 2016;Eddy et al., 2019;Russell et al., 2019), suggesting that olivine LPO is a likely mechanism for explaining the anisotropy. However, the margin-parallel orientation is ∼90° rotated from what is expected for a typical seafloor spreading environment, where corner flow near the ridge generates spreading-parallel olivine LPO that is locked into the lithosphere (e.g., Blackman & Kendall, 2002). Instead, our results suggest margin-parallel shear deformation during continental breakup and initial seafloor spreading. Based on recent seismic anisotropy observations in Pacific and Juan de Fuca lithosphere that deviate from the FSD by 10-70° (e.g., Shinohara et al., 2008;Shintaku et al., 2014;Takeo et al., 2016;Takeo et al., 2018;Vanderbeek & Toomey, 2017), it has been suggested that shear deformation associated with APM at the time of plate formation may dominate the spreading-related deformation if the absolute plate velocities outpace the spreading rate at the MOR (Vanderbeek & Toomey, 2017). This may explain the overall poor correlation between lithosphere anisotropy and FSD in the slow-spreading Atlantic (Becker et al., 2014). At ENAM, estimates of spreading rate prior to the BSMA formation are ultra-slow (∼0.65 cm/yr half-rate), and increase to ∼1.3 cm/yr half-rate just after BSMA formation .
We explore whether our observations of margin-parallel anisotropy can be explained by fast margin-parallel plate velocities at the time of ENAM formation (165-200 Ma) using four recent plate reconstruction models via the GPlates software (Boyden et al., 2011): S12-ESR (Seton et al., 2012); M16-AREPS ; M16-GPC ; M19-T (Müller et al., 2019) (Figure 17). These studies utilize continuously closing topological plate polygon networks that account for inception and cessation of plate boundaries in order to reconstruct plate motions from present day back ∼200 Ma. Relative plate motions in these models are well constrained primarily by seafloor magnetic anomaly picks and fracture zones, particularly in the Atlantic where seafloor is preserved on both conjugate flanks of the MOR. However, absolute plate motions are less well constrained and depend strongly on the choice of absolute reference frame (e.g., Shephard et al., 2012;Torsvik et al., 2008), which varies between studies. Evidence for this is seen in Figure 17, where North America plate motion varies strongly between the four models, particularly for ages 160 Ma. In general, the four models use hybrid reference frames with present day to 70-100 Ma described by a global moving hotspot reference frame and 100-230 Ma constrained by true-polar wander corrected paleomagnetic data. The exception is M19-T, for which authors inverted for a reference frame from 0 to 80 Ma that minimized trench migration velocities and global net lithospheric rotation, and this model is the only one, which includes diffuse deformation along plate boundaries.
Though differences in plate motion do exist between the four models, important similarities are observed. Absolute plate velocities during the ENAM breakup were variable but often fast (1.5-9 cm/yr) relative to spreading (0.5-1.5 cm/yr half-rate), especially early in the breakup . From 170 to 200 Ma, APM directions were significantly different from the spreading direction and often similar in azimuth to ENAM anisotropy (up to ±90° rotated from spreading) for models S12-ESR, M16-AREPS, and M16-GPC. A ∼180° reversal in plate motion direction from approximately north along the margin to south is accompanied by a drop (and shortly followed by an increase) in plate velocity in all models except M19-T. This abrupt change in plate direction occurs at ages ranging from 170 to 190 Ma, depending on the model. Although its origin is not well understood, it roughly correlates in time with far-field plate reorganization processes such as opening of the Gulf of Mexico and the subduction polarity reversal from west-dipping to east-dipping across the Wrangellia Superterrane prior to its collision with western North America (Shephard et al., 2013). As there is a 180° ambiguity in the interpretation of flow direction for a given orientation of seismic anisotropy, both orientations of margin-parallel plate motion are consistent with our observations. At around the time of BSMA formation (∼170 Ma), the spreading rate increased and plate velocities rotated closer to the FSD, on average ( Figure 17). Therefore, fabric east of the BSMA likely represents LPO influenced both by spreading and APM related shearing. This is consistent with the ∼25-45° clockwise rotation in anisotropy east of the BSMA suggested by our 2-D inversions (Figure 9). Although the data do not strictly require this rotation of anisotropy east of the BSMA (Figure 12), it cannot be ruled out. The fast azimuths east of the BSMA and within the ENAM-CSE footprint (165-170 Ma) are not FSD-parallel, indicating that the fabric was still modified by absolute plate motion even when spreading rates were on the order of APM velocities. East of the ENAM-CSE at the younger FAIM experiment (115-130 Ma), APM velocities were fast (2.5-6.5 cm/yr) relative to spreading (1-2.5 cm/yr half-rate) and approximately parallel to the FSD, in agreement with the observed FSD-parallel anisotropy from Gaherty et al. (2004). Given that both the ENAM and FAIM anisotropy observations correlate with APM at the time of spreading, our preferred interpretation is that shear associated with APM dominates LPO formation when plate velocities outpace the spreading rate.

Contribution of the Continental Lid
The inferred gabbroic inclusions frozen within the thin continental lid, if preferentially oriented, could also contribute to the observed margin-parallel anisotropy through a shape-preferred fabric (e.g., Holtzman et al., 2003;Kawakatsu et al., 2009;Schlue & Knopoff, 1976). Shape-preferred models have been invoked to explain margin-parallel structures in active rift environments (e.g., Kendall et al., 2005), where the large velocity contrast between melt-rich dikes and surrounding country rock can produce a significant apparent anisotropy. For the thermal state of the ancient ENAM margin, the available velocity contrast between retained gabbro and surrounding peridotite is much smaller. We calculate the equivalent shear-wave anisotropy of vertical gabbro sheets embedded within a background mantle model via Backus averaging valid for seismic wavelengths much greater than the gabbro channel width (Backus, 1962). We calculate shear-wave azimuthal anisotropy (G) for a vertically Figure 17. Comparison of anisotropy observations with paleo plate motions from 100 to 200 Ma from plate reconstruction models. Solid lines indicate absolute plate motion of North American (NA), while dashed lines indicate relative spreading motion of NA with respect to Africa. (a) Azimuth of NA plate motion relative to the spreading direction. Line colors correspond to four different plate reconstruction models [S12-ESR (Seton et al., 2012); M16-AREPS ; M16-GPC ; M19-T (Müller et al., 2019)]. Gray and tan regions mark the anisotropy fast azimuths at ENAM from this study and the Far-offset Active-source Imaging of the Mantle (FAIM) experiment , respectively. Approximate timing of the ECMA and BSMA emplacements are indicated along the top. (b) NA plate speed (solid) and half-spreading rate of NA with respect to Africa (dashed; M16-GPC, M16-AREPS, and M19-T overlap one another). Gray and tan regions indicate the approximate range of seafloor ages at ENAM and FAIM, respectively. (c) Map of ENAM (yellow) and FAIM (red) OBS with shading that indicates approximate seafloor locations of anisotropy observations. Seafloor age contours from Müller et al. (2008) are shown in black in increments of 5 Myr up to 170 Ma, beyond which their ages are less certain.
oriented alternating layered structure with the relative proportions of gabbro from Figure 16b and find anisotropy of 1% for the continental reference model. Anisotropy peaks directly beneath the Moho because that is where gabbro is most abundant in the model. The Rayleigh-wave anisotropy predicted from the depth-dependent G anisotropy model is 0.5% for the continental reference model, which is far too weak to explain the 2%-4% anisotropy we observe (Figure 18b). We also carry out the calculation for gabbro in the oceanic reference model, which results in slightly stronger G anisotropy of up to ∼2%, but the predicted Rayleigh-wave anisotropy is still much too weak (<1%).
We conclude that while shape-preferred anisotropy due to aligned gabbro channels within the lithospheric lid may slightly contribute to the margin-parallel anisotropy we observe, it is not strong enough to explain the entire signal. Additionally, our inversion tests evaluating 2-D variations in anisotropy suggest a continuity in structure across the BSMA with at most a subtle rotation, arguing against a purely continental-lithosphere origin.
Given the strong relict margin-parallel LPO inferred from shear-wave splitting studies of the eastern U.S. with splitting that largely parallels the structural grain of the Appalachians (e.g., Barruol, Helffrich, & Vauchez, 1997;Barruol, Silver, & Vauchez, 1997;Long et al., 2016;Wagner et al., 2012), it may be tempting to attribute margin-parallel anisotropy west of the BSMA to a preserved accreted terrain. However, predominantly null shearwave splitting observed in North Carolina directly onshore the ENAM-CSE can be interpreted as negligible (or vertical) lithospheric LPO (Lynner & Bodmer, 2017). Furthermore, Rayleigh-wave anisotropy of the onshore region shows weak and laterally variable anisotropy in the lithosphere (Wagner et al., 2018), in contrast to the regionally consistent anisotropy we observe. Olivine LPO formed during collisional orogenic processes would likely be disrupted and altered by later extensional deformation during rifting (Barruol, Silver, & Vauchez, 1997), and therefore, it is unlikely that such a strong, regionally coherent LPO would remain intact at the ENAM-CSE.  Figure 16b. The calculation uses Backus (1962) to estimate the horizontal layered long-wavelength equivalent elastic tensor, which is then rotated 90° to simulate vertical layers. G anisotropy is calculated from the elastic tensor at each depth following Montagner and Nataf (1986). Colors are as in Figure 16. Estimated V P anisotropy is shown in Figure S10 in Supporting Information S1.

Mantle Serpentinization Hypothesis
Serpentinization of the shallow lithospheric mantle via deep cutting normal faults could produce slow velocities and shape-preferred anisotropy (Cai et al., 2018), analogous to gabbro channels. Such alteration has been interpreted from seismic imaging at the non-volcanic Galicia Margin northwest of Iberia (Bayrakci et al., 2016) and has been observed in peridotite samples from various locations in the Atlantic and other slow and ultra-slow spreading regions (Früh-Green et al., 2004). However, serpentinization of the mantle is not typically associated with volcanic rifted margins such as the ENAM with little normal faulting and high mantle temperatures Shuck et al., 2019). Regardless, we explore this hypothesis for explaining the observed low upper mantle velocities and margin-parallel anisotropy.
We perform an analysis similar to that of the gabbro channels above, whereby we mix serpentinite with continental and oceanic endmembers in order to match velocities west of the BSMA. We use constant reference velocities of serpentinite of V S = 2.32 km/s and V P = 5.1 km/s (Cai et al., 2018). As serpentinite is significantly slower than gabbro, only up to ∼5% is required in the continental model and ∼12% in the oceanic model ( Figure S11a in Supporting Information S1). These are equivalent to shear-wave anisotropy (G) values in the continental and oceanic models of ∼1.5% and ∼4%, respectively, assuming 45° dipping faults ( Figure S11b in Supporting Information S1). The oceanic model with ∼4% anisotropy produces surface-wave anisotropy of 0.5%-1.5%, underpredicting our observed values of 2%-4% (Figure S11c in Supporting Information S1).
Normal faulting that extends into the mantle is not observed in the high-resolution reflection images for the ENAM-CSE region . In addition, this mechanism alone is unable explain the high lower-crustal velocities (V P > 7.5 km/s) and thin crust observed west of the BSMA (Shuck et al., 2019). It is also at odds with the relatively bright Moho observed in the region. Taken altogether, we consider this hypothesis unlikely.

Alternative Olivine LPO Fabric Type
Anisotropy that is oriented perpendicular to the spreading direction could instead be indicative of an alternative olivine LPO fabric type. It is typically assumed that olivine with an A-type LPO dominates in the upper mantle with fast [100] aligning parallel to the shear direction and slow [010] oriented perpendicular to the shear plane (Zhang & Karato, 1995). However, recent laboratory experiments on olivine have documented several alternative fabric types with various orientations of the crystallographic axes with respect to the shear plane (Bystricky et al., 2000;Jung et al., 2006;Jung & Karato, 2001;Karato et al., 2008;Katayama et al., 2004).
A commonly invoked mechanism for explaining anisotropy perpendicular to shear is the so-called "a-c switch" (Holtzman et al., 2003;Qi et al., 2018). Experiments on melt-bearing olivine show segregation of melt into bands that leads to strain partitioning that produces LPO with [100] primarily oriented perpendicular to the direction of shear and in the shear plane (Holtzman et al., 2003). Such a fabric could conceivably be responsible for the margin-parallel anisotropy at ENAM, particularly if large volumes of melt were present. However, in their more recent experiments on partially molten samples, Qi et al. (2018) emphasize that melt tends to produce LPO with girdled [100] and [001] (i.e., fast axis dispersed within the shear plane), which has only a slightly stronger [100] alignment perpendicular to shear but is significantly weaker than melt-free olivine LPO. Instead, they predict that the a-c switch would most likely be observed seismically as a transversely isotropic structure with very weak azimuthal anisotropy. As the Rayleigh-wave anisotropy we observe is quite strong (2%-4%), this mechanism is unlikely to explain our observations. A second possibility is the presence of B-type LPO with fast [100] aligning perpendicular to shear and in the shear plane (Jung & Karato, 2001). B-type LPO has been observed in natural samples from a variety of pressure/ temperature conditions (e.g., Bernard et al., 2019) and produced in the laboratory at low temperatures (high stress) and moderate water content (Jung & Karato, 2001). It has been invoked at subduction zones to explain trench-parallel anisotropy observed in the cold mantle wedge (Karato et al., 2008;Nakajima & Hasegawa, 2004) but has not previously been proposed in extensional environments. Given the relatively hot mantle potential temperatures estimated during ENAM formation (1395-1420°C) (Shuck et al., 2019), B-type LPO appears unlikely.

Summary of Interpretation
We favor the interpretation that the lithospheric olivine LPO retains a record of strong margin-parallel APM during early ultra-slow spreading of the Atlantic (Figure 19). This explanation satisfies anisotropy observations at both the ENAM-CSE and FAIM experiments separated in age by ∼50 Myr. It may also explain why we do not observe a significant rotation in anisotropy to FSD-parallel immediately east of the BSMA. It is also likely to produce anisotropy down to at least 60 km depth, which corresponds to the full depth sensitivity of the 17-32 s observations. We cannot rule out an additional small contribution to anisotropy from the thin continental lid between the ECMA and BSMA produced either by preferentially oriented gabbro channels or a relict LPO, but that alone is unable to account for the magnitude and coherence of Rayleigh-wave anisotropy that we observe across the region. Instead, LPO fabric frozen into the lithosphere at slow spreading environments, such as in the Atlantic, is likely to retain a complex deformation signal that records the relative balance between absolute plate motion and seafloor spreading.

Conclusions
We present a shear velocity model of the uppermost mantle lithosphere together with observations of Rayleigh-wave anisotropy offshore the Eastern U.S. that provide new constraints on the late stages of rifting and onset of seafloor spreading at the ENAM (Figure 19). The shear velocity model contains a proto-oceanic domain defined by oceanic crust overlying a 15-20 km thick slow (4.4-4.55 km/s) lithospheric lid interpreted as continental mantle with up to ∼15% gabbro extending from the margin ∼200 km east to the BSMA. East of the BSMA, shallow mantle velocities are fast (>4.6 km/s) and indicative of more typical oceanic lithosphere, suggesting that the BSMA marks the final breakup of the ENAM from West Africa and onset of normal seafloor spreading ∼170 Ma, as previously hypothesized Shuck et al., 2019). Figure 19. Schematic demonstrating our preferred interpretation of the structure and processes that occurred ∼180 Ma prior to the completion of Pangea breakup. The 15-20 km thick continental lid truncated the top of the melting region at the time. As melt percolated through the lid forming the proto-oceanic crust, up to ∼15% gabbro was crystallized within melt channels. Due to the slow spreading rate, corner flow at the ridge was weak relative to plate motion, resulting in olivine LPO oriented sub-parallel to the plate-motion direction in the underlying asthenosphere, which was subsequently frozen into the cooling oceanic mantle lithosphere. The combination of margin-parallel LPO in the oceanic lithosphere and significantly weaker shape-preferred orientation of the frozen gabbro channels in the lid likely contribute to the observed margin-parallel anisotropy.
We observe Rayleigh-wave anisotropy in the lithosphere with a fast direction parallel to the margin, correlating approximately with the plate motion direction at the time, rather than the direction of ultra-slow spreading (≲2 cm/yr half-rate). Nearly 800 km southeast of ENAM at ∼50 Myr younger seafloor, the FAIM experiment showed FSD-parallel anisotropy  that also correlates with the fossil-APM direction. We propose that lithosphere LPO formed at slow-spreading MORs, such as the Atlantic, primarily records mantle shear imparted by absolute plate motion rather than by classic corner flow.

Data Availability Statement
All broadband waveform data was accessed through the IRIS Data Management Center under network code YO using the open-source Python package ObsPy (https://docs.obspy.org/). The script used for retrieving daily waveform data can be found here: https://github.com/jbrussell/fetch_NOISE. Codes used for calculating ambient-noise cross-correlations and extracting phase velocities can be found here: https://github.com/jbrussell/MATnoise. The codes used for 2-D wave propagation can be found here: https://github.com/jbrussell/SEM2D_waveprop. The final shear velocity model and uncertainties are included in Supporting Information S1.