Upper Mantle Anisotropic Shear Velocity Structure at the Equatorial Mid‐Atlantic Ridge Constrained by Rayleigh Wave Group Velocity Analysis From the PI‐LAB Experiment

The evolution of ocean lithosphere and asthenosphere are fundamental to plate tectonics, yet high resolution imaging is rare. We present shear wave velocity and azimuthal anisotropy models for the upper mantle from Rayleigh wave group velocities from local earthquake and ambient noise at 15–40‐s period recorded by the Passive Imaging of the Lithosphere Asthenosphere Boundary experiment at the equatorial Mid‐Atlantic Ridge covering 0–80‐Myr‐old seafloor. We find slow velocities along the ridge, with faster velocities beneath older seafloor. We image a fast lid (25–30‐km thick) beneath the ridge that increases to 50–60 km beneath older seafloor. Punctuated, ∼100 km wide low velocity anomalies exist off‐axis. There are multiple layers of azimuthal anisotropy, including (i) a lithosphere (20–40 km depth) characterized by strong anisotropy (4.0%–7.0 %) with fast‐axes that rotate from ridge subparallel toward the absolute plate motion/spreading direction at distances >60 km from the ridge, and (ii) weak anisotropy (1.0%–2.0%) at >40 km depth. Our results are consistent with conductive cooling of lithosphere, although with some complexities. Thickened lithosphere beneath the ridge supports lateral conductive cooling beneath slow‐spreading centers. Undulations in lithospheric thickness and slow asthenospheric velocities are consistent with small scale convection and/or partial melt. Lithospheric anisotropy can be explained by vertical flow and a contribution from either fluid or mineral filled cracks organized melt beneath the ridge and plate motion induced strain off axis. Deep azimuthal anisotropy is suggestive of upwelling beneath the ridge and three‐dimensional flow possibly caused by small scale convection off‐axis.

The predicted lithospheric thickness increases with the square root of the seafloor age for the half space thermal model (Parker & Oldenburg, 1973), while the plate model was empirically constructed to include the apparent flattening of seafloor subsidence and heat flow at seafloor ages >70 Myr (Parsons & Sclater, 1977;Stein & Stein, 1992). The physical process to explain seafloor flattening is commonly hypothesized to be mantle dynamics such small scale convection (Richter, 1973), which would effectively create a constant heat flux to the base of the plate and prevent further thickening. Surface wave derived shear velocities of the oceanic upper mantle can broadly be explained by a thermal plate, and this is particularly apparent when binned as a function of age (Nishimura & Forsyth, 1988;Ritzwoller et al., 2004).
Although it is widely believed that thermal structure plays a dominant role in the evolution of the ocean lithosphere (McKenzie et al., 2005), there is evidence that the oceanic lithosphere may not be purely thermally defined. Scattered wave studies indicate a strong, sharp lithosphere-asthenosphere boundary (LAB), inconsistent with the gradual base predicted by purely thermal models (Kawakatsu et al., 2009;Rychert et al., 2018b;Tharimena et al., 2017). In addition, the depths of the discontinuities reported by scattered waves, and likely related to the LAB, do not necessarily change with respect to the seafloor age according to thermal predictions, indicating greater complexity (Gaherty et al., 1996;Rychert et al., 2018bRychert et al., , 2012Tan & Helmberger, 2007). Similarly, high resolution in situ surface wave imaging frequently finds lithospheric thicknesses that deviate from classic thermal models (Bell et al., 2016;Gao, 2016;Harmon et al., 2009Harmon et al., , 2020Weeraratne et al., 2007). Ridge studies near Juan De Fuca, 17°S on the East Pacific Ridge and the Mid-Atlantic ridge found asymmetric low velocity regions beneath the ridge, that were not consistent with thermal predictions, with velocities that were several percent slower than what could be explained by temperature alone (Bell et al., 2016;Forsyth et al., 1998;Gao, 2016;Harmon et al., 2020). In the asthenosphere beneath older seafloor, low velocity anomalies that cannot be explained by temperature alone were also observed beneath in the Mid-Atlantic . Also, magnetotelluric (MT) studies observe low resistivities (<10 Ohm·m) at asthenospheric depths both near the MOR (Evans et al., 2005;Johansen et al., 2019;Key et al., 2013;Wang et al., 2020) and further off axis (Naif et al., 2013;Wang et al., 2020) that cannot be explained by thermal models. There have been several alternative hypotheses have been proposed to explain observed deviations in lithospheric thickness, slow velocity anomalies, and/or low resistivities. These include near solidus conditions (Priestly & Mckenzie, 2013;Takeuchi et al., 2020), mantle oxidation state (Cline et al., 2018), water content and/or elastically accommodated grain boundary sliding (Karato & Jung, 1998;Karato et al., 2015) or the presence of low melt fraction (Anderson & Spetzler, 1970). However, to date there have only been a handful of high resolution in situ studies of oceanic lithosphere-asthenosphere systems, so the underlying cause of resistivity and seismic anomalies remains the subject of vigorous debate.
Another important aspect of the evolution of the oceanic lithosphere is how the tectonic plate and asthenosphere deform as the mantle vertically upwells beneath the ridge and then moves laterally in the direction of plate motion. One of the primary tools for studying the accumulated shear strain associated with mantle flow is seismic anisotropy. Anisotropy in the mantle is typically thought to be caused by crystallographic alignment of the fast axes of olivine, the primary mineral component of mantle peridotite, which results in lattice preferred orientation (LPO; Mainprice et al., 2000;Nicolas & Chirstensen, 1987). The shear strain could either have occurred in the past, with anisotropy currently frozen-in to the lithosphere (fossil anisotropy) or the strain and associated anisotropy could be currently accumulating, reflecting present day absolute plate motion (APM). The strength of anisotropy owing to olivine LPO beneath the oceanic plate varies depending on the physical conditions at the bottom of the cooling plate and a spreading rate (Michibayashi et al., 2016). However, particularly near ridges melt generation and segregation processes (Holtzman & Kendall, 2010;Phipps Morgan et al., 1987) may also play an important role for controlling the seismic anisotropy through shape preferred orientation (SPO)of melt aligned cracks. Similarly, SPO of fluid filled normal faults may be an important contribution to the anisotropic signals (Dunn & Toomey, 2001).
Global models and in situ studies of anisotropy yield fast directions that generally support the notion that the strain in the upper mantle is dominated by plate motions, although global resolution is broad, high resolution in situ measurements are rare, and important deviations exist. Global surface wave studies provide a good overview of anisotropy in the oceans. These studies often find shallow fast directions consistent with paleo-spreading directions, with deeper fast directions more aligned with APM (Debayle et al., 2005;Maggi et al., 2006), but often have resolutions of >1,000 km, which is too large to capture variations at the ridge segment scale, typically <100 km. In situ surface wave analysis provides information about the depth of anisotropy at higher lateral resolution than global studies, and broadly finds that the lithospheric anisotropy roughly aligns with the spreading or paleo spreading direction, while asthenospheric anisotropy tends to align with APM, although there are exceptions to this generalization, which are discussed in further detail in the Discussion section (Eilon & Forsyth, 2020;Lin et al., 2016;Russel et al., 2019;Takeo et al., 2018Takeo et al., , 2016. SKS studies, sensitive to accumulated upper mantle anisotropy, have measured fast axis directions that are parallel to the plate spreading and/or APM directions (Harmon et al., 2004;Hung & Forsyth, 1999;Martin-Short et al., 2015;Wolfe & Solomon, 1998). Local shear wave splitting and active source P-wave anisotropy can potentially provide even higher resolution of crust and uppermost mantle anisotropy, where source-station geometries allow for it. These methods have been used to constrain ridge parallel fast axes in the crust and lithospheric mantle within 50 km of the ridge (Barclay, 2003;Blackman et al., 1996Blackman et al., , 1995Dunn et al., 2017Dunn et al., , 2005Toomey et al., 2002).
Here, we present results from the Passive Imaging of the Lithosphere Asthenosphere Boundary (PI-LAB) experiment and the Experiment to Unearth the Rheological Oceanic Lithosphere-Asthenosphere Boundary (EURO-LAB) at the equatorial mid-Atlantic Harmon et al., 2018Harmon et al., , 2020Saikia et al., 2020), which was designed to image the base of the tectonic plate (2018a, 2018b; Rychert & Shearer, 2009). In this paper, we present a new shear wave velocity and azimuthal anisotropy model for the upper mantle using a combination of local earthquake events and ambient noise data. We interpret our results to understand the origin and evolution of the lithosphere-asthenosphere system at a slow spreading center.

Data and Methodology
In the present study, we use data from the PI-LAB experiment, which includes 39, three-component broadband ocean-bottom seismometers (OBS) each equipped with a differential pressure gauge, deployed from March 2016 to March 2017 ( Figure 1) Harmon et al., 2018). We use a vertical component Rayleigh wave seismogram for both local earthquake events and ambient noise cross-correlation. For local earthquakes, we selected those events having magnitudes greater than or equal to 3 (black star in Figure 1) determined from a local catalog. Although initially 39 stations were installed, 2 (I01D and I36D) were not recovered, and 2 had technical errors caused by lack of recording of one or more channels. The cross-correlations were performed between 34 pairs of stations (triangles, Figure 1). For the earthquake data, 18 stations (marked in blue color triangle, Figure 1) recorded good quality signals.

Estimate of Green's Function
To compute noise cross-correlation Green's functions, we use continuous data from the vertical component seismograms. The data processing is performed following Bensen et al. (2007) and Harmon and Rychert (2016). First, the noise time series is split into 24 h segments and resampled to 1 sample per second. Second, the mean, trend and instruments responses are removed from the waveforms. Third, the waveforms amplitudes are normalized using a running mean of the envelope, whitened, and band pass filtered between 0.005 and 0.14 Hz. Finally, the cross-correlation are performed between all possible pair of stations for each day long record and then the cross-correlations for all days was linearly stacked to increase the signal-tonoise ratio (SNR). Figure 2 shows the cross-correlation functions (CCFs) as a function distance. For further analysis, we use the symmetric (real) part of the cross correlations ( Figure 2 Myr (Muller et al., 2008). The blue triangles represent the OBS which recorded the good quality local earthquake signals, and the black triangles the additional stations that were used in the noise analysis. The filled red diamonds show the seismic locations which were not used due to technical errors. The thick red line shows the location of Mid-Atlantic Ridge plate boundary. OBS, ocean-bottom seismometer.

Group Velocity Measurement
To compute Rayleigh wave group velocity, we apply automated frequency-time analysis (FTAN) (Bensen et al., 2007;Levshin & Ritzwoller, 2001) to each vertical-vertical component stacked noise CCF. Similarly, we apply the same FTAN technique on the local earthquake data set to compute Rayleigh group velocity. As an example, an FTAN image is shown in Figure 3b for the Green's function for a pair of stations as well as a stacked FTAN velocity-period plot for long station paths. There was a total of 595 CCF and around 1,600 local vertical seismograms available for measuring dispersion. The main procedure of FTAN is described by Bensen et al. (2007) and Levshin and Ritzwoller (2001). We select only those dispersion measurements that satisfy the followings two criteria: (i) First, the interstation distance should be greater than the three wavelengths at a given period (Bensen et al., 2007), which is given by where Δ is the interstation distance and λ is the wavelength. (ii) Second, we consider those velocity curves having SNR value of greater than 9. The SNR value is defined here based on the ratio of the maximum absolute value in the signal window to the standard deviation of the noise window. Our SNR cut-off is consistent with SNR cut-offs from previous work (e.g., Bensen et al., 2007), and in practice, it eliminated only a few measurements in our data set ( Figure S1). At short periods (4-9 s), we have limited numbers (60-80) of reliable group velocity measurements for the higher mode and similar numbers (<100) for the fundamental mode from 10-to 15-s period. Although good dispersion measurements are possible at these short periods, there are insufficient raypaths for tomography, especially for azimuthal anisotropy so we only present the one-dimensional (1D) group velocity curves in this period range ( Figure S2). The higher mode is visible mainly on the radial component data as evidenced in the stacked period-velocity plot ( Figure S3). The average group velocity dispersion curve across the region is shown in Figure 4a using FTAN. The error shown in Figure 4a is the standard deviation of all the measurements. The individual group velocity curves for the study region are shown in Figure S4.

Surface Wave Tomography
We invert group velocity for tomography, which allows us to combine CCF and local earthquake data in the same inversion. Phase velocity inversion is the topic of future work as it would require good estimates SAIKIA ET AL.   for the source mechanism of the local earthquakes. We use the method of Debayle and Sambridge (2004) to perform the tomographic inversion for group velocity maps for 15-40 s periods. The method follows the continuous regionalization algorithm of Montagner (1986) to retrieve the local structure. Both the isotropic component of speed and 2θ azimuthal anisotropy can be estimated by using this algorithm. To stabilize the underdetermined problem, a priori constraints need to be introduced in this regionalization algorithm. In particular, the algorithm uses a horizontal correlation length along the raypath, L corr , that controls the horizontal smoothness of the model at each geographical point (Montagner, 1986). We chose the L corr value to be 100 km. We tested other values of L corr changing by a factor of two and found the results did not change significantly ( Figure S5). Based on the previous studies of regional (Debayle et al., 2001;Nishimura & Forsyth, 1989) and global scales (Debayle & Sambridge, 2004), an a priori standard deviation is chosen equal to 0.05 km/s.

Azimuthal Anisotropy
For the propagation of Rayleigh waves in a medium, the effect of azimuthal anisotropy can be described with two parameters: cos(2θ) and sin(2θ), where θ denote the azimuth along the propagation path (Montagner & Nataf, 1986). The group velocity, U, as a function of period, T, and azimuth can be written as, where U 0 represent the isotropic group velocity, A c and A s are two anisotropic amplitude parameters.
The strength of anisotropy and the fast propagation direction can be estimated by the following formulas for the S v wave.  We use the continuous regionalization algorithm of Debayle and Sambridge (2004) to estimate the abovementioned anisotropic parameters.

Shear Wave Velocity and Azimuthal Anisotropy With Depth Inversion
We use the linearized iterative inversion method (Hermann & Ammon, 2004) to obtain the shear wave velocity structure beneath the region from the group velocity measurement. The Rayleigh wave dispersion data sets are available in the period range of 15-40 s for a fundamental mode at every 50 km distance interval. The initial 1D seismic velocity model during inversion is obtained from the analysis of the long period teleseismic surface wave data set   (Figure 4b). We account for water depth, with a single layer at the top of our model with  1.5 Vp km/s and  0.0 Vs km/s. Beneath the water layer, the model is parameterized with 1-km thick layers down to a depth of 10 km, followed by 2-km thick layers between 10 and 40 km, and 5-km thick layers between 40 and 120-km depth. In the starting model, we assume 5 km of thick crust and a fixed Vp/Vs ratio of 1.8. During inversion, we varied Vp/Vs ratio within the range of 1.79 and 1.81 but testing with these values indicates minimal changes in the shear wave velocity model. We tested several smoothing values (0-10) and used up to 20 iterations, which was sufficient for the models to converge to an acceptable fit (e.g. RMS < 0.024 km/s). To check the influence of the initial velocity model on the inversion results, we also perform the inversion using different initial models, including PREM, half-layer, spline and a model derived from surface wave model ( Figure S6). Testing with different velocity models indicates that velocity model derived from surface wave analysis  produced the best fitting results (RMS < 0.024 km/s) compared to using other velocity models. The crustal velocities in the final model do not vary significantly from our starting model regardless of damping, so we do not show images of the shallow structure <10 km.
We select two areas (marked by yellow boxes, area 1 near the ridge and area two off-axis) from the study region and invert A c and A s for anisotropic parameters as a function of depth for those two areas. As the depth sensitivity of the anisotropic parameters are derived from the depth sensitivity of group and phase velocity (Montagner & Nataf, 1986), we invert for anisotropic shear velocity based on the group velocity partial derivatives with respect to shear velocity described above, which should be valid for <10% anisotropy. The form of the anisotropy as a function of depth is similar to Equation 1 but substituting shear velocity (V sv ) for group velocity (U), and substituting depth (z) for period (T). We report shear velocity anisotropy in terms of Equations 2 and 3 again substituting dV s for dU. We invert for a smooth parameterization as well as a simplified three layer model. We use a damped least squares inversion, and we determine the error in the smooth model via bootstrap resampling of our data within their error bounds and determining the maximum range of acceptable models. In the simplified three-layer models, the depths of the boundaries below the seafloor are free parameters and are optimized in the inversion. The initial guess for the boundary layer depths is based on the smooth inversions, but testing with different starting depths yields similar structures. We also tested models with single layer anisotropy and two-layer anisotropy; however, these models cannot achieve an acceptable fit to the data. So, the approach is essentially a minimum parameterization.

Resolution Testing
We examine the resolution of our tomography using checker board tests and by examining our raypath coverage. The station and earthquake source distributions for a 20-and 40-s period, is shown in Figures 5a and 5c. In the checkerboard test, we first create square cells of size 1° from period 15-40 s having 6% velocity perturbations about a reference background velocity. We compute the synthetic data set by considering the same model distribution, frequency content and source-receiver distribution as in the observed data set. In addition, to test the stability of the velocity images, we added Gaussian noise with standard deviation of average travel time uncertainty for the calculation of synthetic data (Mottaghi et al., 2013). The same procedures that were applied to the synthetics to generate the checkerboard output, and the synthetic travel times were then used for the tomographic inversion. We use fast marching method (Rawlinson & Sambridge, 2004a, 2004b to compute the travel time residual of synthetic group velocity using the same source-receiver pairs as in the observed data set. The checkerboard test for selective periods (20 and 40 s) indicates that we can resolve our anomalies as small as 1° in majority part of our study region (Figures 5b  and 5c).

Vertical Resolution
The sensitivity of Rayleigh wave group velocity to our average shear velocity structure in our region is shown in Figure 6. For the period range used here (15-40 s), our result has sensitivity to the upper 80 km, with peak sensitivities at our longest periods at ∼40-km depth. We also tested the sensitivity kernel at different periods using modified ak135 model (Hermann & Ammon, 2004). Both the inverted models converge to almost the same value, but we present the case using the best-fit regional model. In addition, to evaluate the vertical resolution, we plot the diagonal components of the formal resolution matrix in Figure   our shear velocity inversion. The average value is around ∼0.1, indicated that we can uniquely resolve the average shear velocity over 10 layers. Specifically, we can resolve the average velocity over the uppermost 8 km and then have a unique piece of information from the inversion every 10-50 km as the layer thickness increases with depth.

Results
Our average 1D group velocity dispersion curve for the study region with the standard deviation of all measurements is shown in Figure 4a for both ambient noise and earthquake data. Both data sets produce similar group velocity values and are within error of each other (Figure 4a). In general, the group velocity increases from 3.66 ± 0.09 km/s at 16 s to 3.95 ± 0.12 km/s at 40-s periods. A posteriori formal error of the tomography model from the final iteration is shown in Figure 7 and the area where the error is close to 0.05 km/s can be considered the poorly resolved area.
The group velocity tomography maps at different periods (16-40 s) are shown in Figure 8. Group velocity varies more smoothly from 3.60 to 3.80 km/s at 16-s period. At 20-s period, group velocity ranges from 3.60 to 3.85 km/s with the fastest region (3.85 ± 0.05 km/s) near −1.00°S latitude and 13.00°W longitude. A pronounced slow group velocity anomaly that varies from 3.60 to 3.80 km/s is consistently found at all periods from 16 to 40 s at the ridge segment between the Chain and Romanche Fracture Zone. The anomaly is approximately ∼150-200 km wide (indicated by yellow box, Figure 8). Another slow group velocity anomaly is visible at the ridge segment south of Chain Fracture Zone, but the anomaly is 0.05-0.10 km/s faster than the segment to the north from 24-36-s period. At 36-40-s periods, there are two off-axis low group velocity anomalies located east of the ridge (3.80-3.90 km/s) ( Figure 8). Otherwise, in the remaining study area, the observed group velocities vary from 3.80 to 4.10 km/s. There are ∼100 km wide areas east of the ridge, where we found slightly faster group velocities that vary from 4.00 to 4.10 km/s. The 3D shear wave velocity model yields an image of a fast lid with variable thickness across the region and a low velocity zone beneath, with several anomalous regions. The shear wave velocity variation at 30km depth is shown in Figure 9a, which provides information about the upper most mantle in our study region. At 30-km depth, a slow shear wave velocity anomaly is visible beneath the ridge segment between SAIKIA ET AL.     16,20,24,28,32,36, and 40 s periods. Colored background shows isotropic group velocities. The red bars indicate the magnitude and fast direction of azimuthal anisotropy. Thick black and thin black (both side) arrows represent the absolute plate motion direction and spreading direction, respectively. Data from yellow boxes are used for the inversion of peak-to-peak azimuthal anisotropy for areas 1 and 2 as labeled, as shown in Figure 11. this profile (A1-A2) away from the ridge anomaly, the velocity variation at 10-60-km depth ranges from 4.75 to 4.50 km/s. There is a maximum velocity at the depth of around 10-25 km between 600 and 1,000 km distance along the profile, where velocity varies from 4.60 to 4.75 km/s. Similarly, along the B1-B2 profile, the velocity variation is in the range of 4.40-4.75 km/s in the upper 50 km of the mantle. There are slow velocities visible near the ridge, but they are off set to the west of the ridge by ∼100 km, and more muted, with a minimum value of 4.45 km/s. Lower velocities are visible at depths >50 km at 800 km distance along the profile, with average values 4.25 km/s. Around ∼60-km depth, velocity drops from 4.40 to 4.20-4.30 km/s over <10 km in both profiles. There are also distinct minima that range from 4.15 to 4.20 km/s near the ridge and also off-axis.
There is a notable change in orientation and strength of anisotropy near the ridge, relative to older seafloor, and similarly between lower and higher periods that is visible in the group velocity azimuthal result (Figure 8). At 16-18-s periods, near the ridge axis the magnitude of the anisotropy is 1%-3% with a fast direction of 100° to 135° for areas with >2% anisotropy, approximately in the strike direction of the ridge. At 16-18-s period further off-axis (>50 km from the ridge), the anisotropy is weak <0.5% across the region. At 20-s period and higher, within 50 km both ridge segments near Chain, the magnitude of anisotropy is relatively weak (<1%), and varies between the spreading direction (79°) and the strike of the ridge segments (−10°). At 20-s period and higher, at distances >50 km from the ridge axis, anisotropy is typically 1.0%-3.0%, although with some weaker areas, and the fast directions rotate toward the spreading direction 78° (around 50% data) and/ or APM 54° (around 25% data). At >20 s, the weaker (0.5%) anisotropy is observed from 20 to 28 s at 8-10°W longitude, but increases to a maximum of 1.5% at longer periods. Also the anisotropy north of Romanche remains consistently weak across all periods at ∼0.5%. The strength of anisotropy along with the deviation of anisotropy orientation from APM, spreading direction and strike of the ridge are shown in Figures S7-S10.
We tested the robustness of the requirement of azimuthal anisotropy, possible trade-offs between isotropic and anisotropic structure, and overall effects on isotropic structure. We present the azimuthal variation of SAIKIA ET AL.
10.1029/2020GC009495 10 of 22 group velocity residual, which indicates a 2θ variability, suggestive of predictions for azimuthal anisotropy ( Figure 10). We also carried out the group velocity inversion both including and excluding anisotropy (2 terms) ( Figure S11). The isotropic structure is very similar between the inversions; however, the magnitude of the velocity variations is stronger in the isotropic case. The data misfit is smaller for the anisotropic case than the isotropic case ( Figure S12).
The azimuthal anisotropy depth inversion for areas 1 and 2 (yellow boxes in Figure 8) are shown together in Figure 11 for comparison. The azimuthal distribution of data for region two is shown in Figure S13. For area 1 near the ridge axis (Figures 11a-11c), the group velocity peak-to-peak anisotropy is up to 1.0% on average for the whole region with fast directions that rotate from 110°-116° at 16-18-s period to 20°-30° at longer periods. We present the magnitude of anisotropy with depth after Equation 2, giving peak-to-peak shear velocity dV s in percent. The resulting smooth model (Figures 11b and 11c, red line) requires a peak of 1.4% ± 2.0% anisotropy in the upper 10 km of the mantle with a fast direction of ∼128°, while the layer beneath requires a peak anisotropy of 1.5 ± 2.0% down to 55-km depth with a fast direction of ∼40°. At greater depths the peak strength of anisotropy is 1.0 ± 1.0% and not significant from zero at depths greater than 60 km, with the fast direction rotating toward the ridge strike direction. We also present a minimum parameterization of a three-layer model (Figures 11b and 11c, blue lines), which has similar fast directions 127°, 32°, and 8°, respectively, and lower anisotropy, 1.9%, 1.5%, and 1.5% with increasing depth, and distributed over similar depth ranges. In area two the average group velocity peak-to-peak anisotropy ranges from 0.5% to 2.0%, with the strongest anisotropy at 24 s period, with fast directions that range from 128° to ∼57° at the longest periods. The smooth model (red line, Figures 11e and 11f) again requires anisotropy in the upper 10 km of the mantle, up to 2.0 ± 1.0% with a fast direction of 175°. From 24 to 53-km depth the peak strength of anisotropy is 3.7 ± 1.0% with a fast direction of ∼72°. From 53 to 120-km depth, the peak strength of anisotropy is 2.0% ± 1.0% with a fast direction centered ∼0° although it wraps around 180° over the depth range. The minimum parameterization (blue line, Figures 11e and 11f) again requires moderately greater strength of anisotropy 3.6% in the 10-22-km depth layer and 5.3% and 1.4% in the 22-42 and 42-112km depth layer, respectively, with fast directions of 159°, 72°, and 0°. At 100 km, where our model suggests very weak anisotropy in both regions, there is some resolution, although it is poor. Forward modeling tests indicate that anisotropy could also be stronger than that of our model in this depth range. In summary, to SAIKIA ET AL.  We have chosen only high quality data with SNR > 15 and interstation distance greater than three wavelengths. The datasets are chosen for the selective area (−16 < longitude < −11; −2 < latitude < 1). The vertical green line defines fast direction at 80°. SNR, signal-to-noise ratio. explain the changes in strength and fast direction in our group velocity peak-to-peak anisotropy data, multiple layers of anisotropy are required and that young seafloor of area 1 has a significantly different structure to older seafloor at area 2.

Shear Wave Velocity Variation
We evaluate our observed lithospheric thicknesses in comparison to the half space cooling model and a geodynamic model which incorporates lateral conductive cooling assuming half-spreading rate of 20 mm/ yr Jha et al., 1994) and mantle temperature of 1,350°C (Figure 12). We estimate the seismic velocity from temperatures using the relationships of Jackson and Faul (2010) assuming a 5.0-mm grain size and accounting for attenuation. We added 5 km to the depths in the models presented in Figure 12 to make them comparable to the depths in our surface wave velocity results, which are with respect to the sea surface. We consider the depth to the gradient in our model at around the 4.50 km/s shear wave velocity contour as the transition of lithosphere and asthenosphere which also roughly coincides with the 1,000°C SAIKIA ET AL.  isotherm. Although the exact choice of contour is arbitrary, it allows us to make comparisons among the observed and synthetic models.
Our observed lithospheric thickness beneath the ridge (25-30 km) (Figures 9 and 12), is thicker than the zero thickness near the ridge of the half-space cooling model and the 14-km thickness in the 2D geodynamic model. Thicker lithosphere is predicted for slow-spreading ridges when geodynamic models account for lateral heat conduction (Parmentier & Morgan, 1990). Therefore, our result suggests that lateral conductive cooling is important near both ridge segments.
Although the lithospheric thickness is similar (around 25-30-km depth) beneath the ridge axis of both the profiles, the two segments have different subridge characters (Figures 9b and 9c). Beneath profile B1-B2 the 4.50 km/s contour remains relatively flat at 25-30-km depth over a broad (±20 Myr, i.e., 600-700 km) SAIKIA ET AL.  (e and f) The difference of temperature and shear wave velocity between these two geodynamic models. Shear velocities are estimated using the parameterizations of Jackson and Faul (2010) that accounts for temperature, pressure and anelastic effects, assuming a 5 mm grain size. 2D, two-dimensional. subridge region (Figure 9c). For comparison, in profile A1-A2 the 4.50 km/s contour is at 25-30-km depth from 10 Myr on the eastern side, extending at least to 10 Myr on the western side (Figure 9b). The model loses resolution on the western side of A1-A2, but if the thickening of the lithosphere is symmetric it would correspond to a 400 km wide region over which the 4.50 km/s contour is shallow. The observed shallow, subridge contours show better agreement with the 2D geodynamic model than the half space cooling model ( Figure 12). Including upwelling over a broader region in the 2D model would likely broaden the region of thin lithosphere beneath the ridge and enhance the match with our observations. Beneath older seafloor, >10 Myr, away from the ridge again we find some similarity, but also important differences in comparison to the model predictions. At 10-20 Myr seafloor age, the half-space cooling and geodynamic models predict the lithosphere, based on the 4.50 km/s contour, should increase from 20 to 40 km thickness, whereas our model increases from 30 to 50 km based on the same contour. Our result is similar or a little thicker than the thermal models, suggesting equal to slightly more cooling and thickening in these regions than predicted by thermal models. At the older seafloor ages (20-40 Myr) the half space cooling and 2D geodynamic model predict thickening to 40-50-km depth. However, we do not observe monotonic lithospheric thickness variation; rather the 4.50 km/s contour undulates within the range of 50-60 km. At the eastern edges of our profiles the lithosphere appears to thin to 30 km based on the 4.5 km/s contour; however, these regions are outside the raypath coverage, and this thinning may be an artifact of the inversion.
When compared to the lithospheric thickness beneath young seafloor (<10 Myr) at other ridges, our observed lithospheric thickness 25-30 km is 10-21 km thicker than the results of Rayleigh wave tomography beneath fast (EPR: ∼6 km) and intermediate (Juan da Fuca: 20 km) spreading ridges (Bell et al., 2016;Harmon et al., 2009). But the result is roughly comparable with the receiver function study that revealed around 20-45-km thick lithosphere near the Juan de Fuca Ridge beneath 0-10 Myr seafloor (Rychert et al., 2018a). Taken together, these results support general thickening of the subridge lithosphere with slower spreading rates caused by conductive cooling (Parmentier & Morgan, 1990). The estimated lithospheric thickness (50-55 km) beneath 10-15-Myr-old seafloor is also consistent that observed at the fast spreading EPR (40-50 km) from Rayleigh wave tomography (Harmon et al., 2009).
We observe several low velocity anomalies (4.15-4.20 km/s) at >55-km depth in both of the profiles onand off-axis. Velocities <4.18 km/s are difficult to explain without the presence of partial melt (Jackson & Faul, 2010;Stixrude & Lithgow-Bertollini, 2005), for example, our 2D geodynamic model predicts minimum velocities of 4.18 km/s. Temperature differences of 200°C might explain these low velocities based on the parameterization of Jackson and Faul (2010). However, geochemical studies suggest that our region is ∼50°C cooler compared to the other sites along the Atlantic which are thought to have an average mantle temperature of 1,350°C-1,400°C (Schilling, 1995;Schilling et al., 1994). A change in grain size might partially explain the result (Behn et al., 2009), but it is hard to imagine why such a difference would exist off-axis and between the different ridge segments. Increased water content in the mantle has been suggested (Karato & Jung, 1998) as a means of reducing the velocity through attenuation effects, although it has been recently suggested that water may not reduce seismic velocities, except in regions where oxygen fugacity is high such as subduction zones (Cline et al., 2018). Near-melt conditions could have an enhanced effect on seismic velocity (Yamauchi & Takei, 2016). However, such conditions do not likely exist over the large ∼100 km swaths of the mantle required to explain our result given predicted geotherms and solidi. Enhanced elastically accommodated grain boundary sliding, beyond that present in the model of Jackson and Faul (2010) used here could also reduce velocity via attenuations effects (Karato et al., 2015). However, the presence of partial melt is the most likely explanation as the partial melt is expected beneath the ridge, due to the upwelling of mantle material. Off axis, upwelling and partial melt may be the result of small-scale convection (Richter, 1973) as discussed in detail in Harmon et al. (2020). We also prefer the partial melt explanation due to evidence from other techniques such as the colocated MT study which also found evidence for partial melt in the mantle . In other words, since seismic and MT are both sensitive to melt, but MT lacks sensitivity to most of the subsolidus mechanisms, melt is the most likely explanation.
There is also a marked difference in the asthenospheric anomalies between the two ridge segments within our study region, which also might indicate differences in the amount of partial melt present or thermal structure in the shallow asthenosphere at 30-50-km depth. The ridge segment north of Chain has a minimum velocity 4.27 km/s versus the segment south of Chain which has a minimum velocity of 4.45 km/s, ∼5% faster, and resolvable outside of error (Figure 9b). This would require a substantial difference in temperature between the two ridges as discussed above. The presence of a greater amount of partial melt fraction in the northern segment, ∼0.5%, based on a 7.9% velocity reduction for 1.0% melt fraction (Hammond & Humphreys, 2000) or 2.0% melt fraction using less conservative estimates of the effect of melt on shear velocity (Clark & Lesher, 2017;Schmelling, 1985) can easily explain our result.
Velocities in the asthenosphere in our region (50-80 km depth) are faster (4.15 km/s) compared to other MOR systems such as Eastern Lau Back Arc spreading center (∼3.6 km/s; Wei et al., 2015), Southern East Pacific Rise (∼3.9 km/s; Harmon et al., 2009), and the Juan De Fuca/Gorda Ridges (∼3.8 km/s; Gao, 2016), although Lau may be complicated the presence of volatiles from the Tonga subduction system. Differences between our observation and the observations in the Pacific can likely be explained by the arguments we have made previously about the relative amount of partial melt for the difference between ridge segments in our study region.

Significance of Anisotropy
Overall, the depth dependent shear velocity anisotropy models suggest highly variable structure, likely owing to a combination of past shear strain and current deformation resulting in olivine LPO, and also reflecting both subridge upwelling and plate motion related shearing (Figures 9 and 11). There is a strong difference in the fast direction beneath the ridge and older, more distant lithosphere in the group velocity anisotropy, and also a general increase in strength from ∼0.5% near the ridge to 2.0%-3.0% off-axis. Changes in anisotropic fast direction and strength with period show variability that reflects at least two and possibly three distinct layers, likely related to one to two layers of frozen in lithospheric anisotropy, and also an asthenospheric layer (Figure 11). Testing indicates that these layers, characterized by strong variations in anisotropic fast direction with depth, are required to accommodate the subtle changes in fast directions in the group velocity anisotropy due to the broad peaks in the group velocity depth sensitivity kernels (Figures 4c and S14). For instance, a single layer of anisotropy of 4% anisotropy with the fast direction in the spreading direction at 20-40-km depth can explain the magnitude of the peak-to-peak anisotropy in area 2 reasonably well but does not predict the subtle rotation in fast direction (blue lines, Figure S14). Shallower or thicker layers do not produce a similar pattern in the magnitude of area 1 or 2 (green and red dashed lines Figure S14). On the other hand, multiple layers allow for rotation of the fast direction with period (black dotted, Figure S14) as is observed in the inversions. Despite the heterogeneous model, there are also large regions of general agreement that can be well-represented by a subridge region (area 1) and older seafloor (area 2), which we interpret in greater detail.
Near the MOR at asthenospheric depths >25 km, for example, in our area 1 (Figures 11a-11c), vertical shear is predicted to align the olivine a-fast axis in the vertical direction, and therefore weak azimuthal anisotropy may exist owing to preferential alignment of the second fastest axis, the c-axis, in the ridge-strike direction, in comparison to the slowest axis, the b-axis, in the ridge perpendicular direction (Blackman & Kendall, 2002;Blackman et al., 2017Blackman et al., , 1996. For buoyant flow, predicted at slow ridges, this effect is predicted to be more pronounced than at fast ridges in 2D geodynamic models (Blackman et al., 1996). This mechanism could explain our anisotropy at 25-110 km depth in area 1 beneath the ridge, where the required anisotropy is between 1.0% and 2.0% on average with fast directions oriented within 30° of the N-S direction, that is, closer to ridge parallel.
At the shallowest depths (10-20 km) in area 1 near the ridge axis, our result that has a fast direction that is partially between the ridge strike and spreading direction could be explained in a variety of ways. One possible physical explanation is that the surface wave result reflects recent (<8 Myr) shearing in the upper most mantle and or crust, possibly related to readjustment of the tectonics near the ridge and the Romanche transform (Bonatti et al., 1994;Harmon et al., 2018). Alternatively, it could represent the effects of complex 3D mantle flow, which are predicted to be important at slow spreading centers (Parmentier & Morgan, 1990). Another explanation is that the result is crustal, related to normal faulting. Finally, it could reflect averaging of on-axis and off-axis structures. The map-view image at 16-s period suggests this represents an average of very small, ridge parallel anisotropy beneath the ridge axis and relatively strong anisotropy (2.0%-3.0%) just off-axis rotated slightly off the ridge axis direction at ∼120°, which is the dominant direction we observe.
The idea of strong changes in fast directions within 30 km of the axis is similarly supported by local and teleseismic shear wave splitting results. Evidence from the local S wave splitting study from the PI-LAB data suggests that on the ridge axis, there are of ridge parallel fast directions, and within 30 km (within area 1) a few fast direction results that are oriented at ∼120° (Kendall et al., 2019). SKS splitting results near the ridge also find evidence for up to 3 s splitting in the ridge parallel direction (Kendall et al., 2019), which could reflect SPO from melt bands (Kendall, 1994;Kendall et al., 2019). Source side splitting measurements also suggest both strong (>2 s) ridge subparallel fast direction splitting near the ridge axis, as well as spreading direction fast axes and nulls, from different back azimuths, also suggesting complexity possibly due to multiple layers of anisotropy caused by 3D flow (Eakin et al., 2018;Nowacki et al., 2012). In summary, considering map-view and the independent constraints from shear-wave splitting, the just-off ridge-parallel surface wave fast directions can be explained by an average of ridge parallel directions on-axis and the influence of one of the other mechanisms off-axis and averaged with the SPO from melt bands by the broad surface wave sensitivity.
Beneath older seafloor, that is, area 2 (Figures 13d-13f) the 1.0%-2.0% anisotropy oriented roughly N-S at >42-km depth may be caused by upwelling or downwelling from small scale convection . Its muted magnitude could be related to a lack of organized flow or because it reflects the difference between the b-and c-axes, similar to the deep anisotropy in area 1. The deformation history of the mantle rocks may also be important, as the mantle flow transitions from upwelling near the ridge to more horizontal shearing away from the ridge. Laboratory studies suggest a weakening of anisotropy and potentially complex rotations of fast directions as the olivine fabric evolves (Boneh & Skemer, 2014).
In area 2, the ∼5% anisotropy with a fast direction of 76° at 22-42 km depth (Figures 13e and 13f), is consistent with olivine a-axis alignment in the spreading direction and/or APM and therefore, likely due to shear strain caused by plate motions. The anisotropy is likely frozen-in to the lithosphere given that this anisotropy lies within the fast cold lithosphere and is likely to be fairly strong. Our observation would translate to ∼0.3 s vertical S-wave split time, at least partially explaining nearby local S wave splitting and SKS splitting results which find ∼0.5 s splitting times and fast directions in the APM/spreading direction (Kendall et al., 2019). Therefore, the magnitude or depth extent of the anisotropy in our model may be damped owing to surface wave resolution, or broader lateral surface wave sensitivity could be averaging more complex structures. In addition, the shallower and deeper ridge parallel fast directions would need to be canceled by deeper perpendicular anisotropy to fully reconcile the result with SKS, a topic of future work.
At the shallowest depths (<20 km) in area 2, the 3.0%-4.0% anisotropy with a 160°-170° fast direction could be explained by shearing in the mantle just beneath the Moho related to a frozen-in signature of complex 3D flow just beneath the Moho, potentially different than that occurring today, given the different fast direction observed at similar depths beneath the ridge today, for example, 120° in area 1. Alternatively, SPO of fluid or mineral e.g. serpentine, filled cracks in the crust can explain the ridge parallel fast directions (∼160°-180°) we observe. Previous active source studies on the Mid-Atlantic ridge have observed ridge parallel fast directions in the upper 4-6 km of up to 3.0% P wave anisotropy (Dunn et al., 2017(Dunn et al., , 2005, although admittedly weaker than the 2.0%-3.0% S wave anisotropy we observe in the upper 10 km of the mantle. These results were attributed to fluid filled fractures near the ridge. However, at older ages, some of the fractures may be cemented by alteration products such as anisotropic serpentine and talc minerals, which may enhance the anisotropy in the upper 20 km if these fractures are altered to great depth (Faccenda et al., 2009). If anisotropy at <20 km is explained by crustal cracking phenomena such as this, it also suggests that tectonic stresses operating by the ridge are much different today than what we observe in area 1, likely because of current readjustment, in order to give significantly different anisotropic fast directions. Our depth inversion result is also generally consistent with the local S wave splitting study results, which are confined to the upper 20-30 km based on earthquake depths, in area 2 where there are two populations of fast directions, in either the spreading direction (∼79°) or a smaller population in the ridge strike direction (169°), with 0.1-0.3 s split times (Kendall at el., 2019).
In both areas, where we performed the anisotropy inversion for depth, the strongest required anisotropy was located in the lithosphere (0-42 km) rather than in the asthenosphere, which is contrary to predictions for strain development of olivine from state of the art geodynamic models (Blackman & Kendall, 2002;Blackman et al., 2017Blackman et al., , 1996. However, several active source studies (Dunn et al., 2017(Dunn et al., , 2005Gaherty et al., 2004;Kodaira et al., 2014) including the classic Raitt et al. (1969) study have observed strong 3.0%-9.8% P-wave azimuthal anisotropy in the uppermost mantle beneath Moho, suggesting that strong olivine alignment occurs at very shallow depth, which agrees with our findings. One possibility is that shear zone may form at the base of the crust, which aligns olivine a-axes into the spreading direction (Kodaira et al., 2014). However, this process is not well understood as precise quantification of low temperature plasticity appropriate SAIKIA ET AL.
10.1029/2020GC009495 17 of 22 Figure 13. Comparison of lithospheric and asthenospheric anisotropy from in situ surface wave studies beneath oceans. Red colored symbols indicate asthenospheric results (>60 km depth) and blue symbols indicate lithospheric results (<60 km). We report the maximum strength within each layer and its associated direction. We translate various anisotropic parameterizations in the literature to dV s . Symbol shape indicates results from studies as shown in the key. (a) Peak to peak magnitude of anisotropy as a function of seafloor age. (b) Fast direction (ϕ) relative to spreading/ paleospreading rate for the lithospheric results or relative to APM for the asthenospheric results as a function of seafloor age. (c) Same as (a) but as a function of spreading/paleospreading rate derived from Muller et al. (2008). (d) Same as (b) but as a function of spreading/paleospreading rate. (e) Same as (a) but as a function of APM rate with both lithosphere and asthenosphere measurements differenced from APM, and (f) same as (b) but as a function of APM rate with both lithosphere and asthenosphere measurements differenced from APM. APM, absolute plate motion.
for sub-Moho conditions and its effects on crystallographic alignment and olivine deformation are still an actively evolving area of research (Mei et al., 2010).
The fast axis rotation from ridge parallel at the ridge to spreading direction from 0 to 12 Myr at ∼40-km depth allows us to examine the rheology of the mantle lithosphere (Figure 9). We can estimate approximate mantle strain rates and viscosities for the base of the lithosphere from these observations. Our result suggests that it takes 12 Myr for the anisotropy to align in the spreading direction from the ridge (area 1) to area 2. Assuming it takes a strain of roughly four to five to reorient an aggregate into a strong anisotropic fabric aligned in the shear direction (Skemer et al., 2012), gives a strain rate of 1.1-1.3 × 10 −14 s −1 that yields a lower lithospheric viscosity of ∼8 × 10 20 to 8 × 10 21 for tectonic stresses ranging from 10 to 100 MPa.
The anisotropic structure we observe in the slow spreading Mid Atlantic Ridge provides a complementary measurement for our global understanding of the structure of the oceanic lithosphere. We compare our results to several results using Rayleigh waves measured by OBS from the Pacific (Eilon & Forsyth, 2020;Lin et al., 2016;Russel et al., 2019;Takeo et al., 2014Takeo et al., , 2018Takeo et al., , 2016) as a functions of seafloor age, spreading/paleo spreading rate and APM rate in a no net rotation frame of reference (Figure 13). For each study, we take the maximum value of the magnitude of anisotropy in the lithosphere and asthenosphere and the corresponding fast direction. For comparison purposes, we convert anisotropy estimates in a variety of parameterizations the literature to dV s in % as presented for our result. We subtract the spreading directions from the lithospheric fast directions and the APM directions from the asthenospheric fast directions in the seafloor age and spreading rate comparisons, and we subtract the APM directions from both layers in the APM rate comparison.
We find that in general that magnitude of anisotropy in the lithosphere and the asthenosphere is similar across all studies, with the lithosphere ranging from 3.0% to 7.0% with a mean value of 4.1%, while the asthenosphere ranges from 1.0% to 5.0%, with a mean value of 2.2%. The lithospheric anisotropy is greater than the asthenospheric anisotropy except in one case the Juan De Fuca Ridge (Eilon & Forsyth, 2020). There does not appear to be any relationship between the magnitude of APM and the strength of lithospheric anisotropy (correlation coefficient 0.47; Figure 13e). There is a weak positive correlation between the strength of lithospheric anisotropy and spreading rate and seafloor age (correlation coefficient 0.61 and 0.57, respectively, Figures 13c and 13a), although this correlation is mostly due to the strong anisotropy reported for the ultrafast paleo spreading rate for the NW Pacific (Muller et al., 2008;Takeo et al., 2014). There is no clear relationship between the strength of the asthenospheric anisotropy in comparison to seafloor age, spreading rate, or APM rate (correlation coefficients < 0.3; Figures 13a-13e). Our result is in general agreement with the lack of a trend in the magnitude of azimuthal anisotropy with age from global dispersion curves from Rayleigh Wave phase velocities, at least outside of error in one study (Smith et al., 2004) and also the relatively constant magnitudes from 0 to 80 Myr for short periods (25-75 s) in another (Eddy et al., 2019). A weak decrease in the magnitude of anisotropy beneath seafloor >80 Myr from these studies is not observed in the few points in our compilation. However, lithospheric properties from regional studies typically exhibit more variability than trends gleaned from age averages of global models . Comparisons between P wave anisotropy from refraction work in the western Atlantic suggested weak lithospheric anisotropy, interpreted as the result of slower spreading, for instance, in comparison to the Pacific (e.g., Gaherty et al., 2004). This is in general agreement with the weak correlation with spreading rate in our compilation, primarily reflecting the strong anisotropy beneath the oldest seafloor in the western Pacific (Takeo et al., 2014). Therefore, spreading rate may exert a weak control on the development of lithospheric anisotropy as it directly impacts strain development, while asthenospheric anisotropy is relatively insensitive to the plate motion rate.
The fast directions in the lithosphere are within 30° of the spreading/paleospreading directions, which is also in line with the ∼20°-30° fast direction error reported in many studies, with a few deviations (gray region, Figures 13b-13f). The deviations include our study near the ridge, the Juan De Fuca Ridge result near the ridge (Eilon & Forsyth, 2020) a result near Tahiti (Takeo et al., 2016) and another near the Shatsky Rise (Takeo et al., 2018). The near ridge environment may be complicated by upwelling and a very young lithosphere, while the regions near Shatsky and Tahiti may have been affected by hotspot interaction. The general agreement suggests the lithospheric fast direction is typically determined by plate spreading.
Similarly, many of the asthenospheric fast directions relative to APM in a no net rotation reference frame are also are within 30° of each other, again within the ∼20°-30° error reported in many studies (gray region, Figures 13b-13f). In the asthenosphere our fast directions and those from the Shatsky Rise are outside the nominal error. This may reflect upwellings due to small scale convection (Takeo et al., 2018). Overall, we conclude that in general the spreading direction appears to be commonly frozen into the lithosphere, and the APM is common in the asthenosphere, although there do appear to be exceptions. The exceptions hint at an important role for processes such as hotspot volcanism, small scale convection and 3D flow near ridges.

Conclusions
We find shear wave velocities that are slower along the ridge and tend to increase with sea floor age. Thickened lithosphere (25-30 km) beneath the ridge supports the influence of lateral heat conduction at slow spreading centers. Slower velocities beneath the northern ridge segment in comparison to the southern segment are likely explained by recent upwelling and potentially more melt. The lithosphere generally increases in thickness with age but also has significant thinning and thickening relative to simple thermal models. There are several punctuated low velocity regions in the asthenosphere that likely caused by partial melt. Group velocity azimuthal anisotropy maps generally support relatively weak (1.0%) ridge-subparallel fast directions within 50 km of the ridge that rotate to plate spreading/APM directions further off axis and increase in strength to (1.0%-3.0%). However, depth inversion of anisotropic structure suggests multiple layers of anisotropy are required, with a weaker asthenospheric anisotropy (∼1.0%-2.0%) likely related to complex convection and/or vertical upwelling. Stronger lithospheric anisotropy (2.0%-5.3%) could be caused by a combination of plate shearing off-axis and shallow fluid or anisotropic mineral filled fractures. Whereas, the similarly strong lithospheric anisotropy near the ridge could be influenced by SPO caused by melt or fluid filled fractures or LPO due to 3D upwelling or complex tectonic readjustment. The deviations from simple thermal models highlight the role of lateral conductive cooling, small scale convection, and melt in the evolution of the oceanic lithosphere.