Spatial Variations in Crustal and Mantle Anisotropy Across the North American‐Caribbean Boundary on Haiti

Haiti, on the island of Hispaniola, is situated across the North American‐Caribbean plate boundary at the transition point between oblique subduction in the east and a transform plate boundary in the west. Here we use shear wave splitting measurements from S waves of local (0–50 km) and intermediate depth (50–150 km) earthquakes as well as SK(K)S phases from teleseismic earthquakes to ascertain good spatial and vertical resolution of the azimuthal anisotropic structure. This allows us to place new constraints on the pattern of deformation in the crust and mantle beneath this transitional region. SK(K)S results are dominated by plate boundary parallel (E‐W) fast directions with ~1.9 s delay times, indicating subslab trench parallel mantle flow is continuing westward along the plate boundary. Intermediate depth earthquakes originating within the subducting North American plate show a mean fast polarization direction of 065° and delay time of 0.46 s, subparallel to the relative plate motion between the Caribbean and North American plates (070°). We suggest a basal shear zone within the lower ductile crust and upper lithospheric mantle as being a potential major source of anisotropy above the subducting slab. Upper crustal anisotropy is isolated using shear wave splitting measurements on local seismicity, which show consistent delay times on the order of 0.2 s. The fast polarization directions indicate that the crustal anisotropy is controlled by the fault networks in close proximity to the major strike‐slip faults, which bisect the north and south of Haiti, and by the regional stress field where faulting is less pervasive.


Introduction
The island of Hispaniola has seen a complex tectonic evolution due to its position along the North American-Caribbean Plate boundary (Figure 1), with its present-day position being at the transition from oblique subduction in the east to oblique collision to the west. Haiti, which forms the western half of Hispaniola, therefore has a high seismogenic potential, with events such as the M w 7.02010 earthquake occurring periodically during its recent geological history (Bakun et al., 2012;Prentice et al., 2010). Observations of the present-day crustal stress state and its link to the underlying mantle dynamics of the plate boundary play an important role in our understanding of the plate tectonics.
Seismic anisotropy, defined as the directional dependence of seismic wave velocities, provides a crucial link between geophysical observations and internal deformation within the Earth's crust and mantle. One of the most common methods to observe seismic anisotropy is shear wave splitting, whereby the shear wave becomes polarized into two orthogonal components as it passes through the anisotropic medium. By measuring the orientation of the fast polarization direction (φ) and the delay time (dt) between the components inferences can be made about the anisotropic system. In the mantle, seismic anisotropy is typically interpreted as the result from strain-induced lattice-preferred-orientation (LPO) of olivine caused by past or present mantle flow (Nicolas & Christensen, 1987). In the crust there are however, many sources of anisotropy which affect measurements of shear wave splitting (Crampin & Peacock, 2008). Two common mechanisms include the alignment of microcracks within the stress field, or the alignment of large-scale structures such as faults.
In this study we have calculated shear wave splitting using data from 31 seismic stations to achieve the highest spatial density of shear wave splitting measurements across Haiti to date. In addition, we perform the analysis on earthquakes from multiple depth ranges to isolate anisotropy in the crust, mantle lithosphere, and subslab mantle. This gives us both spatial and vertical resolution in our measurements. These observations of deformation beneath the island help to constrain both its tectonic history and the geodynamics that underpin the region.

Tectonic Setting
Hispaniola is situated along the plate boundary between the North American and Caribbean plates (Figure 1), which are obliquely converging at 19 mm/year (Demets et al., 2000). Along this plate boundary oblique subduction of the Atlantic oceanic lithosphere occurs under the Greater Antilles, from Puerto Rico through to central Hispaniola. Seismotectonic and geodetic observations both suggest that the subduction interface beneath Puerto Rico is uncoupled, with little stress being partitioned into the upper crust Calais et al., 2016;Symithe et al., 2015). Hispaniola marks the transitions from oblique subduction in the east to oblique collision in the west with significant motion on left-lateral strike-slip faults, the Enriquillo-Plantain Garden fault zone (EPGF) in the south and the Septentrional-Oriente fault zone (SOF) in the north (Leroy et al., 2015;Mann et al., 1995). Compressional structures both on-and offshore, such as the Northern Hispaniola Fault (NHF), also accommodate the collision (Dolan et al., 1998;Mann et al., 2002). The subduction interface along the North Hispaniola Fault is by contrast coupled, which imparts significant stress into the upper crust, as evidenced by the active strike-slip and compressional faulting Calais et al., 2016).
Dense networks of GPS stations in kinematic studies have been used to define three main blocks between the Caribbean and North American plates with major active faults coinciding with the block boundaries Manaker et al., 2008;Symithe et al., 2015). To the west, the Gonâve block is bounded by the SOF in the north and the EPGF in the south. The eastern extent of this block is still debated but has most recently been suggested lie within the Plateau Central-San Juan Valley region across Haiti and the Dominican Republic (Symithe et al., 2015). The main Hispaniola block is bounded by the western part of the Muertos Trough in the south, by the SOF in the north, and the Mona Passage to the east. GPS observations also suggest an additional block in Northern Hispaniola between the SOF and the NHF (Symithe et al., 2015). The FIGURE 1. Regional map of the North American Caribbean plate boundary. Black lines represent the major plate boundary faults. Colored bars represent previous shear wave splitting measurements, scaled by the delay time. Black arrow represents the relative plate motion between the two plates. Colormap and contours show the Caribbean subducting slab model from Slab2 (Hayes et al., 2018). Contours for the shallow subduction interface (<40 km) have been extended westward beyond the model as has been inferred from seismicity (Symithe et al., 2015). EPGF = Enriquillo-Plantain Garden fault zone; SOF = Septentrional-Oriente fault.
Puerto Rico-Virgin Island (PRVI) block in the east is bounded by the eastern part of the Muertos Trough in the south and the Puerto Rico trench in the north.
The collision of Hispaniola with the Bahamas platform has been used to explain the lateral segmentation of the plate boundary (Leroy et al., 2000;Mann et al., 2002). At present, this collision slows the eastward motion of Hispaniola, with respect to the North American plate, due to the compressional forces imparted (Calais et al., 1992;De Zoeten & Mann, 1999). The free PRVI block continues to move eastward at the velocity of the Caribbean plate. This idea is consistent with a coupled subduction interface beneath the Hispaniola block compared to an uncoupled interface beneath the PRVI block.

Crustal Anisotropy
The preferential orientation of cracks in crustal rocks has long been thought to be a common mechanism by which to induce shear wave splitting in the crust (Crampin, 1981). However, the mechanism by which cracks align has been the subject of much debate. Crampin et al. (1984) proposed the idea of extensive-dilatancy anisotropy (EDA), where subcritical crack growth occurs, extending cracks in the direction of maximum horizontal compressive stress (SH max ) (Atkinson, 1984). Thus, the fast polarization direction of the shear wave split should be parallel to the regional SH max and this has been observed in numerous studies (Crampin & Booth, 1985). By contrast, other studies have observed spatial variations in the fast polarization direction over much smaller distances than would be expected by the EDA hypothesis. This has led to alternative origins for observed anisotropy in the crust, which include the following: alignment of fractures proximal to active faults (Tadokoro et al., 1999;Zhang & Schwartz, 1994;Zinke & Zoback, 2000), intrinsic anisotropy from rock fabrics-associated sedimentary bedding planes (Kern & Wenk, 1990), alignment of minerals in metamorphic fabrics (Brocher & Christensen, 1990), and remnant palaeostress features (Aster & Shearer, 1992;Blenkinsop, 1990). In reality, there is likely to be an interplay between two or more of these mechanism beneath any individual station that likely contribute to the observed crustal anisotropy through shear wave splitting.
There is also still some debate over the depth extent of anisotropy which local shear wave phases are sensitive to. Many studies argue that the observed crustal anisotropy is shallow, likely in the upper 2-3 km, to explain sharp variations in splitting parameters only a few kilometers apart (Peacock et al., 1988;Savage et al., 1989;Savage et al., 1990). Shallow crustal anisotropy has also been favored in a number of fault zones such as the San Andreas fault system (Zhang & Schwartz, 1994) and the North Anatolian fault system (Eken et al., 2013;Peng & Ben-Zion, 2004). By contrast, studies have observed a clear increase in delay ties with depth (Keir et al., 2005;Li et al., 2015;Shih & Meyer, 1990), which suggest a more pervasive anisotropy throughout the crust.

Mantle Anisotropy
In the upper mantle anisotropy is typically interpreted as the result of lattice-preferred orientation (LPO) of olivine, which develops due to dislocation creep or recrystallization (Karato et al., 2008;Mainprice et al., 2005;Zhang & Karato, 1995). Under typical upper mantle conditions (A-, C-, or E-olivine fabrics) the fast axis (as sensed by vertically incident shear waves) tends to align parallel to the direction of shear in the mantle. In regions of high stress, high water content, and low temperature, the slow c axis can align with the shear direction (Jung & Karato, 2001), which is known as B-type fabric. If this relationship between the LPO and deformation field is known, observing anisotropy can help us place constraints on the dynamics of the upper mantle. Shearing in the mantle is most commonly associated with mantle flow, but can also be the result of deep shear zones beneath large crustal faults.

Shear Wave Splitting in the Caribbean
There have been a number of studies that have focused on the anisotropic structure of the Caribbean-North American plate boundary, mostly focusing on mantle anisotropy through observations of shear wave splitting from teleseismic earthquakes (Figure 1). In general, fast polarization directions have been found to be parallel to the plate boundary (Benford et al., 2012;Hodges & Miller, 2015;Lynner & Long, 2013;Meighan & Pulliam, 2013;Piñero-Feliciangeli & Kendall, 2008;Schlaphorst et al., 2017). However, given the lateral segmentation of the plate boundary, differing explanations have been used to explain these observations. Beneath the PRVI block and to the east where oblique subduction is occurring, trench-parallel mantle flow is suggested. Moving west, as indicated in GPS studies , Hispaniola marks a transition from oblique subduction to strike-slip faulting in the western segment of the plate boundary. As a result of this change plate boundary parallel fast polarization directions in the Dominican Republic have been associated to large shear zones beneath the major strike-slip faults (Benford et al., 2012).
At present, however, no measurements of shear wave splitting have been made on Haiti. Using data from the Trans-Haiti network and permanent seismic stations, we extend the lateral extent of shear wave splitting observations along the Caribbean-North American plate boundary ( Figure 2). This dense network of seismic stations deployed across the island allows us to observe variations in the anisotropic structure of the highly segmented plate boundary in this region. Using a local catalog of earthquakes, recorded by the Trans-Haiti Network stations (Possee et al., 2019), we are also able to observe in detail spatial variation of anisotropy in the crust and examine the results in relation to detailed GPS models of crustal stress and known fault structures, which will help to constrain crustal and upper mantle motions in the region.

Seismic Data
In this study, we use seismic data from 31 broadband seismic stations deployed on Hispaniola, 30 of these deployed on Haiti and 1 deployed on Dominican Republic ( Figure 2). Twenty-seven of the stations are part of the Trans-Haiti project, a temporary deployment active from April 2013 to June 2014 (https://www.fdsn. org/networks/detail/ZU_2013/). The remainder of the station are permanent seismic station deployed on the island, three from the Canadian Seismograph Network (https://doi.org/10.7914/SN/CN) and one from the United States Geological Survey (USGS)(https://doi.org/10.7914/SN/CU). For these four stations, we use seismic data in the time frame from January 2010 to December 2016. All data are available from the Incorporated Research Institutions for Seismology Data Management Center (IRIS-DMC).

Local Earthquakes
For the analysis of crustal anisotropy, we use a local catalog of 1,058 earthquakes recorded using the Trans-Haiti Network seismic stations from 2013-2014 (Possee et al., 2019). The earthquakes were manually picked and hypocenter locations determined using a 1-D velocity model for the region ( Figure 2) (Possee et al., 2019). At each station we limit our analyses to direct S phase arrivals with an incidence angle of less than 35°, calculated using ray tracing through the 1-D velocity model. In this way, we ensure all arrivals are within the shear wave window and avoid any phase conversions at the free surface that could cause error in measurements (Crampin & Booth, 1985). This resulted in a total of 4,134 direct S-phase arrivals on which shear wave splitting analyses were performed.
In addition, to observe any possible temporal changes in crustal anisotropy we also analyze local earthquakes from the aftershock sequence of the 2010 M w 7.0 earthquake . This analysis is limited to station LGNH, which was the only station in close proximity active during this time period (January to March 2010).

Teleseismic Earthquakes
To observe mantle anisotropy, we analyze shear wave splitting on SK(K)S phases from teleseismic earthquakes at a distance of 85-140°( Figure 2). The P-to-S conversion this phase undergoes at the core mantle boundary polarizes this phase in the plane of the raypath, which mean any shear wave splitting observed is limited to the receiver-side and any source side splitting effects can be ignored. Due to the increased noise associated with island seismic stations we limit our search to earthquakes with a M w > 6.0. Events were then manually screened to remove noisy events where no clear phases could be identified. For stations from the Trans-Haiti network this search criteria and screening resulted in 27 possible earthquakes on which to perform the shear wave splitting analysis. The extended active time period of the permanent stations allowed us to increase this to 117 earthquakes for those stations.

Intermediate Depth Earthquakes
In order to improve the depth resolution of our study we also analyze intermediate depth earthquakes in the range of 50-150 km depth from the NE Caribbean subduction zone. The maximum crustal thickness across Hispaniola is~40 km as seen in receiver function studies (Corbeau et al., 2017), therefore most of these intermediate earthquakes occur within the subducting North American slab. They therefore allow us to discriminate anisotropy from the mantle wedge versus subslab effects. Earthquakes are selected with a M w > 4.0 from the Advanced National Seismic System Comprehensive Catalog (ANSS ComCat), which contains earthquake both located by the USGS NEIC and regional networks (https://doi.org/10.5066/ F7MS3QZH). As with the local earthquakes, the events selected are all within the shear wave window to avoid any phase conversions at the surface.
To ensure a high signal-to-noise (SNR) ratio we use earthquakes with a M w > 4.0 which are taken from the ComCat earthquake catalog. As with the local earthquakes the events must be within the shear wave window to avoid any phase conversions at the surface.

Station Corrections
In shear wave splitting analysis it is critical to know the orientation of the sensor correctly. Therefore, we re-evaluate the station orientations using the software package Doran-Laske-Orientation-Python (DLOPy), which isolates Raleigh waves in a ray-based coordinate system on the radial component to find the correct station orientation (Doran & Laske, 2017). This method was primarily developed for ocean bottom seismometers (OBS) where station orientation is unknown but provides an excellent method to check for misaligned land stations.
For permanent island stations such as LGNH and SDDR the analysis of surface wave through DLOPy finds orientations ±2°of north for the north component ( Figure S1a in the supporting information). This both confirms the orientation of these stations but also validates the method for use on this data. When applied to stations of the Trans-Haiti network we find that the majority of stations were orientated correctly, though some were misaligned by up to 30°with respect to north (Table S1, Figure S1d). We therefore applied a rotation correction to any station where the misalignment was greater than the 4σ error calculated using bootstrap analysis in DLOPy. A full list of results can be found in Table S1.

Shear Wave Splitting Analysis
All shear wave splitting analyses were performed using the MFAST code (Savage et al., 2010) The code extends the shear wave splitting procedure of Teanby et al. (2004) to include an automated grading criterion to assess the quality of measurement. MFAST allows multiple filters to be tested on a window surrounding the S-phase pick to determine which produces the highest signal-to-noise ratio (SNR). The SNR is calculated as the ratio of the S-phase amplitude in the signal window, to the amplitude of noise in a window proceeding the S-phase arrival. The three band-pass filters with the highest SNR ratio are selected for splitting analyses. . Examples of good shear wave splitting measurements made using MFAST for (a) local S phase recorded at station JACM (b) SKS arrival recorded at station SDDR. Left-hand panels show waveforms rotated to the incoming polarization direction (p) and its perpendicular direction (p ⊥ ), which is equivalent to the radial (R) and transverse (T) component for SKS phases. The original filtered waveforms are on top and the corrected waveforms beneath. The shaded gray region represents the chosen window to perform the shear wave splitting measurement on. The central panels show the linearization of particle motion before (left column) and after (right column) correction for the chosen measurement window. The right panel shows the contours of the smallest eigenvalue of the covariance matrix for the preferred measurement, blue cross shows the final value for the two splitting parameters.
The final filter was selected as the one which produced the highest-quality splitting measurement. Events are rejected if there is no filter producing an SNR > 3. The final filters used for each event are recorded alongside each splitting measurement in the supporting information (Data Sets S1-S3). The splitting procedure implemented is based on the method of Silver and Chan (1991), which minimizes the eigenvalue of the covariance matrix of horizontal particle motion by using a grid search over a range of fast polarization directions (φ) and delay times (dt). For S phases from local events we search over delay times from 0-0.4 s, the range is increased to 0-1 s for intermediate events, and for SK(K)S phases we search over 0-4 s. This is performed on multiple time windows to avoid the subjective decision of manually selecting a time window and the cluster analysis procedure from Teanby et al. (2004) is used to select the final optimum shear wave splitting parameters.
For direct S phases from local earthquakes we use the automated grading criteria of the MFAST code to select high-quality measurements, which avoids any potential bias that could be introduced in manual procedures (Figure 3a). This procedure gives a final overall grade based on the cluster analysis, measurements relationship to incoming polarization direction, cycle skipping, and contours from the eigenvalues of the covariance matrix for the best measurement. For more details on this grading procedure, see Savage et al. (2010). For teleseismic phases we acknowledge that the MFAST automated grading procedure was primarily designed for use on local direct phases, therefore we apply additional manual quality control to the measurements to select high-quality measurements and also to identify null measurements. A high-quality measurement should therefore have (1) a clearly visible shear wave pulse with an amplitude larger than the noise, (2) removal of significant energy on the corrected transverse component, (3) elliptical particle motion which is linearized after the correction, and (4) a well-defined minimum in the map of smallest eigenvalues of the covariance matrix from a grid search over splitting parameters (Figure 3b).

SKS and SKKS Results
We obtained 150 good shear wave splitting measurements across 21 of the seismic stations analyzed in this study (Data Set S1). Seventy-eight of these measurements were obtained from two permanent stations (LGNH, SDDR), the remainder were obtained from the Trans-Haiti temporary network. To first order all measurements show a consistent fast orientation with a mean φ of 110°± 28, which show little variation moving from north to south across the plate boundary (red bars, Figure 4). Delay times, which provide an indication of the magnitude of anisotropy are also consistent across the region with a mean value of 1.9 s ± 0.9. We did not observe any significant number of null splitting measurements at any of the stations analyzed in this study.
At station SDDR due to the increased number of measurements, 54, we were able to observe possible backazimuthal variation in the splitting parameters, although as can be seen in Figure 2 the events tend to cluster at backazimuths of~50°, 250°, and 330°and do not sample a continuous range. This variation was particularly evident in the fast orientations, but also in the delay times despite being more scattered ( Figure S2). Similar patterns could be seen at station LGNH, however, observations were less clear at this station due to the reduced number of measurements.

Local S Results
We obtained over 600 shear wave splitting measurements on S phases from local earthquakes within the crust (Data Set S2). Spatially, these measurements are concentrated in the south as detected seismicity rates are much higher in the region beneath each station that defines the shear wave window (Figure 2) (Possee et al., 2019). However, we have still made at least 10 shear wave splitting measurements at 23 of the 31 seismic stations used in this study. One hundred fifty-one of these measurements belong to the aftershock sequence of the 2010 M w 7.0 earthquake at station LGNH. Compared to SK(K)S results measurements of the fast orientation from local S phases were more variable with many changes observed across the network (Figure 4).
Stations in southern Haiti show the most variation in fast directions with large differences even among stations that are in close proximity to each other ( Figure 5). Broadly fast directions can be described as being NW-SE or NE-SW. The exception to this general rule is station BELC, which is located just south of the main EPGF and the 2010 M w 7.0 earthquake hypocenter. The delay times at all stations are consistent with mean values all ranging from 0.11-0.16 s.
Stations across central Haiti show variations in splitting parameters; however, compared to either southern or northern Haiti we observe a dominance in fast orientations with a NNE-SSW trending direction ( Figure 6). The exception to this trend is station PIGN, which shows a dominant NW-SE fast direction. Events for station PIGN are all sourced from north of the station, so the anisotropic signature being sampled is likely from northern Haiti as opposed to central Haiti. Mean delay times for these stations are consistent with those in the south ranging from 0.12-0.20 s.
In northern Haiti we were only able to achieve 10 or more shear wave splitting measurements at two stations, BOIS and MEND. Both show different results, BOIS shows a single dominant fast orientation at NW-SE, while MEND shows two dominant fast directions NW-SE and NE-SW. Mean delay times are still consistent though ranging from 0.16-0.17 s.

Intermediate S Results
From intermediate earthquakes ranging in depth from 50-150 km we were able to make 17 good quality shear wave splitting measurements (Figure 4; Data Set S3). All these were measured from earthquakes to the east of the station SDDR, likely from the subducting slab. The mean fast polarization direction from these measurements was 065°with a mean delay time of 0.46 ± 0.17 s.

Discussion
It is clear from our splitting analyses of local, intermediate, and teleseismic S waves that distinct zones of anisotropy must exist within the crust, supraslab mantle, and subslab mantle beneath Hispaniola indicating that multiple anisotropic domains are present. This is suggested by the step-wise increase in delay time with depth (Figure 7), as well as the variation in fast polarization directions between local, intermediate, and SKS measurements (Figure 4). This immediately rules out a single uniform source of anisotropy commonly sampled by all raypaths. Instead, combining shear wave splitting measurements from multiple event depth ranges, we are able to build a large-scale model of seismic anisotropy consisting of several anisotropic domains from the subslab mantle through to the crust. Hispaniola marks the transition from oblique subduction under the PRVI block in the east, to oblique collision across the Hispaniola and Gonâve blocks on the North American-Caribbean plate boundary. As a result, seismic anisotropy in the crust and lithosphere beneath the Hispaniola block reflects the changing tectonic forces along the plate boundary, such as the influence of the Caribbean plate's collision with the Bahamas Platform. Anisotropy in the subslab mantle, by comparison, reflects the underlying geodynamics of the regional mantle flow in the NE Caribbean. We discuss the characteristics and implications of the anisotropic patterns within each zone in the sections that follow.

Mantle Anisotropy
The shear wave splitting results from both SK(K)S phases and intermediate depth earthquakes associated with the subduction system, described in the results, show that the anisotropic nature of the mantle cannot be described simply by a one-layer model. This is seen in both the delay time measurements from each event  Figure 4). Rose diagrams of fast polarization directions are plotted at each station. Individual measurements are plotted as bars oriented by their fast direction and at a distance half way between the event and the station. Black lines represent surface faults mapped according to geological observations (Lambert et al., 1987). SH max for this region is estimated at 010°from Calais et al. (2016). type (Figure 7) and the backazimuthal variations we observe in the SK(K)S splitting measurements at station SDDR ( Figure S2), which typically indicate a multilayered or dipping anisotropic fabric (Savage, 1999). Benford et al. (2012) proposed a simple model in which mantle anisotropy was dominated by large-scale lithospheric fabrics associated with the large transform fault zones that bisect the north and south of the island. However, further to the east other studies prefer subslab trench-parallel mantle flow with an isotropic mantle above (Lynner & Long, 2013;Meighan & Pulliam, 2013;Piñero-Feliciangeli & Kendall, 2008;Schlaphorst et al., 2017). Both these mechanisms would generate seismic anisotropy with an east to west trending fast polarization direction. It is therefore hard to distinguish between these mechanisms using only SK(K)S phases, since they are sensitive to anisotropy throughout the entire mantle and crust beneath the seismic station.
The intermediate earthquakes from the subducting slab allow us to make some additional constraints. If the primary source of anisotropy is located in the lithospheric upper plate, then intermediate earthquakes should have delay times approaching those of the SK(K)S phases. Figure 7 clearly shows that beneath station SDDR intermediate earthquakes primarily sampling the mantle wedge have delays times of~0.5 s, much less than the delay times observed from SK(K)S phases~1.6 s. This indicates the bulk of the delay time is accumulated in the deeper part of the mantle (below 150 km), consistent with a subslab mantle anisotropic source that exhibits trench-parallel mantle flow. Studies to the east, in the region of Puerto Rico, suggest the mantle wedge above the slab does not provide a significant anisotropic contribution compared to the subslab mantle (Piñero-Feliciangeli & Kendall, 2008;Schlaphorst et al., 2017). Schlaphorst et al. (2017) observed a mean dt measurment of 0.21 s from seismicity associated with the slab. However, the mean dt we observe at station SDDR for these intermediate earthquakes is 0.46 ± 0.17 s, compared to 0.19 ± 0.09 s for earthquakes originating within the crust (Figure 7). This implies that there is 0.27 s of dt being accumulated between the base of the seismogenic zone and the top of the slab. In Puerto Rico, geodetic studies have suggested the subducting slab is de-coupled from the overriding crust, whereas in Hispaniola it is much more coupled leading to the large transpressive fault zones, such as the EPGF, SOF and NHF Calais et al., 2016;Symithe et al., 2015). This change in coupling and segmentation along the North American-Caribbean plate boundary has been related to the collision of the Caribbean plate with the Bahamas Platform (Leroy et al., 2000;Mann et al., 2002;Symithe et al., 2015). The increase in dt we observe from the intermediate earthquakes is therefore indicative of this change in kinematics moving along the plate boundary.
When interpreting the depth of intermediate earthquakes, we assume that the majority will have originated within the subducting slab. Since there is no obvious increase in dt with depth (Figure 7), it suggests that anisotropy is not distributed uniformly throughout the mantle wedge, but instead more likely resides in a relatively thin layer that all rays sample equally. We propose two possible explanations for this layer, the first being a layer above the slab with a fabric induced by the simple shear deformation associated within the entrained mantle above the slab. This mechanism has been modeled to produce up 10% anisotropy with a fast axis parallel to the down-dip direction of the slab (Faccenda & Capitanio, 2013).
The second hypothesis is a shallower layer, with a fabric formed by lithospheric shear zones associated with plate boundary faulting, as suggested by Benford et al. (2012). Assuming the magnitude of anisotropy is consistent with an upper mantle composition and fabric (5%) (Mainprice & Silver, 1993), this would require these shear zones to be approximately 25 km thick. The deepest local earthquakes, with depths >40 km, also show shear wave splitting with an increase in dt compared with those in the shallow crust (Figure 7). This suggests that the shear zone may also be present in the lowermost part of the crust where deformation is dominantly ductile as opposed to brittle. A lithospheric shear zone beneath Hispaniola demonstrates the stress transfer, resulting from the collision of the Caribbean plate with the Bahamas Platform, across the plate boundary and transpressional faulting observed on the island.
The fast directions of these intermediate events have a mean of 065° (Figure 4), which is close to the relative plate motion between the Caribbean and the downgoing North American plate,~070° (Demets et al., 2000). A fast orientation subparallel to relative plate motion is likely consistent with either of the two proposed layers of anisotropy above the slab and the base of the lithosphere. One difference that would be testable between the two layers is in the dip of the anisotropic layers. The lithospheric shear zone is likely to have a fast axis that is horizontal, whereas any fabric induced by the subducting slab should be dipping with the slab, likely 20-30°. However, the difference in splitting parameters between a flat versus shallowly dipping layer would be expected to be small (e.g., Eakin et al., 2018) and would require a wide range of backazimuths, from 0-360°, to be able to distinguish between the two. With a limited number of results from intermediate earthquakes, we therefore lack the range of measurements required to explore the parameter space necessary to discriminate between a flat anisotropic layer at the base of the lithosphere versus a dipping layer of

10.1029/2019JB018438
Journal of Geophysical Research: Solid Earth anisotropy just above the slab. It therefore remains to be investigated as more data becomes available in the future.
Since Haiti lies at the transition from subduction to a transform plate boundary it is useful to compare our observations to other similar systems. One of the most studied examples is the transition from the Hikurangi subduction zone into the oblique Alpine fault system in New Zealand (Wallace et al., 2012). Beneath northern and central New Zealand, where subduction is present, teleseismic SK(K)S and local S phases have been used to infer trench-parallel mantle flow both above and beneath the subducting slab (Audoine et al., 2000;Marson-Pidgeon et al., 1999). We observe that the same subslab mantle flow, however, do not see evidence for mantle flow above the slab. The majority of the mantle wedge beneath Hispaniola appears to be isotropic, with discrete anisotropic layers caused by lithospheric shear or shear just above the slab. As in our study, Marson-Pidgeon and Savage (2004) observe a wide range of shear wave-splitting parameters at individual stations above the transitional region and infer complex layered anisotropic related to absolute plate motion in the lower layer and geological features/shear zones in the upper layer. However, despite improved backazimuthal coverage and incorporating dipping layers, two-layered models were still shown to be nonunique.
On the South Island of New Zealand, SKS fast polarization directions have shown the shear zone associated with the strike-slip fault zone need to be 100-200 km wide (Zietlow et al., 2014). Local and regional earthquakes were then used to confirm the shear zone is likely located in the upper lithosphere (Karalliyadda et al., 2015). This is consistent with the lithospheric shear zone we interpret beneath Hispaniola and suggests it may be present across the whole island and possibly connects the EPGF in the south and the SOF in the north, which are approximately 200 km apart. In other transform plate boundaries, such as the Northern Anatolian Fault zone (NAFZ) and San Andreas Fault (SAF), anisotropy has been attributed to deeper asthenospheric mantle flow. In the SAF a two-layer system is favored with fast polarization directions in the lower layer subparallel to absolute plate motion and the upper layer subparallel to the surface trace of the SAF (Hartog & Schwartz, 2001). In the NAFZ a single-layer model is favored with a 150 km thick anisotropic layer present within the asthenosphere with the mantle lithosphere not thought to be a major contributor (Biryol et al., 2010;Sandvol et al., 2003).

Crustal Anisotropy
Shear waves recorded from local earthquakes within the crust clearly show that the crust is pervasively anisotropic across Haiti and the plate boundary. There is also clearly heterogeneity in the splitting parameters spatially across Haiti, suggesting either competing mechanisms for anisotropy or that one mechanism is systematically varying across the island. The delay times, however, show little spatial variation or any significant correlation with depth at individual seismic stations ( Figure S3). To first order, this indicates that anisotropy is not uniformly distributed throughout the crust across Haiti as delay times do not systematically increase with depth.

Fault-Controlled Anisotropy
To assess the relationship between fault structures and shear wave splitting parameters we compare all known mapped surface structures to our results (Figures 5 and 6). Along Haiti's southern peninsula, in close proximity to the EPGF, two primary sets of fault structures can be observed, one trending at~060°and the other at~130°. At seismic stations in this region there is good agreement between the bulk of local S splitting directions and the structures which are near to the seismic station ( Figure 5), indicating these structures may be the dominant control on the observed anisotropy.
Station LGNH is situated directly over the most seismically active faults therefore making it the ideal place to test this hypothesis. When the local S splitting results are separated out by their backazimuth/fault zone origin, there is a good agreement with the orientation of different fault systems (Figure 8b). Just south of LGNH seismicity is associated with the steeply north dipping Léogâne fault which trends 75° (Calais et al., 2010;Douilly et al., 2013), seismicity from this fault produces fast polarization directions that parallel this fault orientation. Seismicity to the west of LGNH, likely originate from the Trois Baies fault system, which is southward dipping and trends at 130° (Momplaisir, 1986). Shear wave splits from these events have a mean fast polarization direction of 150°, close to that of the fault system. These two fault systems result in a bimodal appearance to the final rose diagram of fast directions for the station LGNH (Figure 8b). We therefore suggest that much of the heterogeneity in southern Haiti could be explained by the interaction of the S waves with multiple fault structures, which are extensive throughout the crust in southern Haiti.
We also attempt to look for any temporal changes in anisotropy, comparing shear wave splitting results from the 2010 aftershock sequence of the M w 7.0 earthquake  to those made on events during the Trans-Haiti network window of 2013-2014 (Possee et al., 2019). Events recorded in the 2010 aftershock sequence have a mean fast polarization direction of 132°with a delay time of 0.09 ± 0.09 s. In 2013-2014 there was little change in the mean fast polarization direction which is 125°and a delay time of 0.11 ± 0.08 s. The small changes observed are more likely related to the differences in event-receiver pair geometries as opposed to any changes in stress or strain during the time period.
In northern Haiti, stations BOIS and MEND seem to show a similar pattern to those stations in southern Haiti, with fast directions that reflect the complexity of the mapped fault structures near to each station ( Figure 6). The measurements recorded at station PIGN, which dominantly sample the anisotropic signature of northern Haiti, also show a fast direction consistent with the fault structures that they sample ( Figure 6). The SOF bisects the region just north of island (Figure 4), we would therefore expect to see this return to structurally controlled anisotropy as the distance to this major strike-slip fault zone decreases.

Stress-Controlled Anisotropy
Across central Haiti there is less agreement between the NNE-SSW preferred local S fast polarization directions at stations in relation to the mapped fault structures (Figure 6). This may be due to the apparent lower density of faults in the region that is visible in Figure 6. Another cause of crustal anisotropy is EDA, where microcracks in the crust align parallel to SH max (Atkinson, 1984;Crampin et al., 1984). The fast polarization direction should therefore also be parallel to SH max . Geodetic and seismotectonic studies show that the SH max has a NNE-SSW orientation across much of central Haiti (Calais et al., 2016;Corbeau et al., 2019;Rodriguez et al., 2018), which is generally more consistent with our measurements of shear wave splitting in the crust (Figure 6).
The stations MGOA and PETG (Figures 2 and 4), which lie further west along the southern peninsula, also show NE-SW fast polarization directions. The anisotropic mechanism for these observations would therefore be better explained by the EDA hypothesis than by structurally controlled aligned faults/fractures. This is consistent with the observation that seismic activity increases toward the east along the EPGF with present-day activity peaking at the intersection of the EPGF with the NW-SE trending compressional fault systems (Possee et al., 2019).

Conclusions
In this study we have made shear wave splitting measurements at 31 seismic stations dominantly on Haiti but also in the Dominican Republic, from teleseismic, intermediate depth, and local earthquakes. The measurements provide new insights into the anisotropic nature of the crust and upper mantle beneath the island of Hispaniola along the North American-Caribbean plate boundary. Figure 9 shows in schematic form, our interpretation of the anisotropic structure with depth based on the observed shear wave splitting measurements. We identify multiple distinct source regions of anisotropy beneath Hispaniola. The deepest of these lies beneath the slab, in the subslab mantle, for which we infer a continuation of trench parallel mantle fabrics, consistent with what has been observed along much of the Lesser Antilles subduction zone to the east. Moving upward, the mantle above the slab appears to be mostly isotropic with the exception of an anisotropic layer, which we suggest could be explained by either an entrained shear fabric just above the subducting slab, or a lithospheric shear zone associated with plate boundary faulting. In either scenario, the observation of an anisotropic layer between the slab and the seismogenic zone beneath Hispaniola is consistent with previous modeling of geodetic data Calais et al., 2016), which has suggested that the subduction system is highly coupled in this location.
Overall, the crust itself is widely anisotropic and while observations from crustal earthquakes are heterogeneous, to first order, pervasive crustal faulting appears to be the dominant source of anisotropy in the south and north of the island when close to the active fault zones. Particularly on the southern peninsula of Haiti, the anisotropic signature recorded is dominated by the primary faults associated with the 2010 M w 7.0 Haitian earthquake (Figure 8). There is no temporal change though between the immediate aftershocks of the 2010 earthquake and those deployed during the Trans-Haiti deployment window (2013)(2014). Across central Haiti we observe an overall rotation of fast directions to parallel the maximum horizontal stress orientation suggesting that anisotropy is primarily stress controlled in this region. (1) Shallow anisotropy in the crust related to either the orientation of crustal faults or the maximum horizontal stress (SH max ). (2) A possible basal shear zone representing a ductile mantle lithosphere where the dominant source of anisotropy is the lattice preferred orientation (LPO) of olivine. (3) Possible thin shear zone above the subducting slab orientated parallel to the relative plate motion between the Caribbean and North American plate. (4) Subslab trench parallel mantle flow with anisotropy being caused by the LPO of olivine. 19 mm/year arrow represents the relative plate motion between the Caribbean and North American plates.