Seismic Discontinuities Across the North American Caribbean Plate Boundary From S‐to‐P Receiver Functions

The evolution of the Caribbean plate has resulted in the formation of volcanic arcs, the Caribbean Large Igneous Province (CLIP) and micro‐plates across the plate boundary zones. The northern plate boundary with the North American plate has been particularly segmented with the transition from oblique subduction to oblique collision moving from east to west. However, there are few constraints on the seismic structure of the upper mantle across the plate boundary. Here we use S‐to‐P receiver functions to map seismic velocity discontinuities across the plate boundary, placing constraints on crustal and lithospheric thicknesses, as well as the structures associated with subduction and collision. We image a velocity increase with depth, consistently seen at 28–34 ± 4 km along the plate boundary, which corresponds to the Moho. A second strong velocity increase with depth is observed at depths of 64–66 ± 5 km, which is related to the presence of subducting slabs and anisotropic effects. We image a velocity decrease with depth at 95–135 ± 7 km, which reflects a lithosphere‐asthenosphere boundary that varies in depth across the plate boundary. The deepest negative discontinuity spatially maps to the CLIP. We suggest that a deep melting depth at 135 km, associated with an elevated potential mantle temperature of 1585 ± 20°C during CLIP formation, caused a depleted and dehydrated root to the base of melting, thus thickening the lithosphere.

• We image a crustal thickness of 28-34 km on islands across the North American-Caribbean plate boundary • We image a thicker lithosphere beneath the Caribbean Large Igneous Province (CLIP) compared to the remainder of the plate boundary • Elevated mantle temperatures and deeper melting during formation of the CLIP resulted in a thicker depleted and dehydrated lithospheric root in seismic velocities (Jackson & Faul, 2010). However, the depth and strength of observed seismic discontinuities beneath the oceans and continents has led to suggestions of a change in composition and/or the presence of partial melt, across which, the strength of the mantle changes (Kawakatsu et al., 2009;Rychert et al., 2005). In particular, velocity decreases with depth observed beneath both large igneous provinces, which are typically thought to be caused by previous hotspot volcanism, and current day hot spots are too deep and sharp to be explained thermally (Byrnes et al., 2015;Rychert, Harmon, & Armitage, 2018;Tharimena et al., 2017). This has led to suggestions of a compositional boundary that reflects the depth extent of depletion during melting (Gaherty et al., 1999;Hirth & Kohlstedt, 1996;Yamamoto & Morgan, 2009).
The formation of the CLIP is widely associated with the Galapagos hot spot (Escuder-Viruete et al., 2016;Hastie & Kerr, 2010;Herzberg et al., 2000; J. L. Pindell & Kennan, 2009). High potential mantle temperatures of 1560°C-1620°C associated with the CLIP (Herzberg et al., 2000) and its eruption onto old unaltered Pacific oceanic lithosphere, make the CLIP a likely place to observe a thickened lithosphere associated with a depleted mantle root at its base. The CLIP comprises a large proportion of the Caribbean plate (Mauffret & Leroy, 1997), therefore a thicker lithosphere beneath the CLIP is likely to have important implications for the behavior of the Caribbean plate in response to subsequent tectonic events, such as the collision with the Bahamas platform.
In this study, we use passive seismic data from stations across the North American-Caribbean plate boundary to determine S-to-P (Sp) receiver functions to investigate the structure of seismic discontinuities beneath the plate boundary, such as the Moho and LAB. We migrate the receiver functions to depth in 3-D, which allows us to map the spatial extent of discontinuities and relate them to the tectonic terranes and structures that are present across the plate boundary zone.

Geological Setting
The Caribbean plate is moving NE relative to the North American plate at a rate of 19 mm/yr (DeMets et al., 2000). In the east, this is accommodated by the oblique subduction of Atlantic oceanic lithosphere from Puerto Rico through to eastern Hispaniola. This subduction transitions to oblique collision in central Hispaniola with significant left-lateral strike-slip motion on the Enriquillo-Plantain Garden fault zone (EPGF) and Septentrional-Orient fault zone (SOFZ). West of Jamaica, these faults bound the Cayman Trough, a ∼1,600 km long and ∼150 km wide marine depression formed as a pull-apart basin in the early Eocene (50 Ma), which eventually resulted in a full oceanic spreading center along the mid-Cayman Ridge.
The bulk of the Caribbean plate offshore is comprised of the CLIP with a crustal thickness of 10-20 km, imaged by seismic refraction surveys (Granja-Bruña et al., 2010;Mauffret & Leroy, 1997;Mauffret et al., 2001;Núñez et al., 2016). Thin crust is observed in the Venezuelan basin (10-15 km), whereas thick crust is observed under the Beata Ridge and the Lower Nicaraguan Rise (15-20 km). It can also be seen outcropping in southern Hispaniola and Jamaica along the plate boundary (Dürkefälden, Hoernle, Hauff, Wartho, et al., 2019;Escuder-Viruete et al., 2016;Hastie et al., 2008;Hoernle et al., 2002;Révillon et al., 2000;Sinton et al., 1998). Geochemically, the CLIP is known to be heterogeneous with both depleted and enriched compositions, consistent with other LIPs such as the Ontong Java and Manikil Plateaus. The islands across the North American-Caribbean plate boundary (Puerto Rico, Hispaniola, Jamaica, and Cuba) are dominantly formed from GAC tectonic terranes formed during the SW-dipping subduction of oceanic crust under the Caribbean plate from the late Cretaceous (∼100 Ma) through to early Eocene (∼50 Ma). P-to-S receiver functions show that crustal thicknesses across the islands typically vary from 25 to 30 km, increasing up to 40 km in central Hispaniola due to compression associated with the collision of the Caribbean plate with the Bahamas Platform (Corbeau et al., 2017;González et al., 2012).
The NW of the Caribbean plate is formed of the Chortìs block and the Lower and Upper Nicaraguan Rise. The Chortìs block is a Precambrian Craton that comprises much of Central America and extends into the marine setting in the form of the Nicaraguan Rise. The eastern extent is not well known, however, samples from the Lower Nicaraguan Rise have been inferred to be of CLIP origin (Case, 1991;Dürkefälden, Hoernle, Hauff, Werner, & Garbe-Schönberg, 2019;Mauffret & Leroy, 1997), suggesting the eastern limit may be along the Pedro Bank fault zone, which separates the Lower and Upper Nicaraguan Rise (James, 2007;Lewis et al., 2011;Sanchez et al., 2019). Seismic refraction surveys have shown the Upper Nicaraguan Rise has a crustal thickness of 20-25 km (Edgar et al., 1971), compared to 15-20 km for the Lower Nicaraguan Rise (Mauffret & Leroy, 1997). The eastern limit of the Lower Nicaraguan rise is defined by the Hess Fault zone.

Evolution of the Caribbean Plate
In most reconstructions, the Caribbean plate formed by the "Pacific model," in which the plate, GAC and CLIP formed in the Pacific Ocean above the Galápagos hot spot (Burke, 1988;Hastie & Kerr, 2010;Hoernle et al., 2002;J. L. Pindell, 1990). There is much observational evidence that supports this model, compared to the "inter-American" model (Ball et al., 1969;Meschede & Frisch, 1998), which cannot easily explain the relative age of arc magmatism, lack of recent sea-floor spreading in the Caribbean plate and the area of Caribbean lithosphere that has been subducted beneath South America.
In the "Pacific model," the GAC initially formed from NE-dipping subduction of the Farallon plate under a rifting seafloor between North and South America associated with the break-up of Pangea. An arc-polarity reversal then occurred with subduction switching to SW-dipping, The exact timing of this reversal is still debated with two end-member scenarios, the first being reversal during the Aptian-Albian (125-100 Ma) (J. Pindell et al., 2005; J. L. Pindell & Kennan, 2009), and the second being a Turonian-Campanian (93-70 Ma) reversal (Hastie et al., 2013). In this time frame, plume related activities associated with the Galápagos hot spot led to widespread volcanism and the formation of the CLIP on the Farallon/proto-Caribbean Plate. It is well accepted that the primary formation of the CLIP occurred at ∼89 Ma; likely secondary volcanic events have also been interpreted prior to the main event between 140-100 Ma (Whattam & Stern, 2015) and after, at ∼75 Ma (Hoernle et al., 2002(Hoernle et al., , 2004 and ∼62-65 Ma (Révillon et al., 2000).
If a later Turonian-Campanian reversal age is accepted, the collision of the buoyant CLIP lithosphere with the GAC is what may have led to the subduction reversal (Hastie et al., 2013). This explains the proximity of the CLIP to the GAC, with exposure of CLIP material seen on southern Hispaniola and Jamaica, through fragments of CLIP material being obducted onto the GAC during collision. In the earlier Aptian-Albian reversal scenario, the reversal is explained through inter-American transform faults. In this scenario, the observed present day proximity of the CLIP to the GAC is explained through either a "slab gap," associated with the subduction of the spreading ridge between North and South America (J. Pindell et al., 2012;J. L. Pindell & Kennan, 2009), which allowed the plume material to rise close to the GAC. Alternatively, it has been explained by the later collision and strike slip tectonics associated with the collision of the GAC with the Bahamas Platform during the Late Paleocene (∼55 Ma; J. Pindell et al., 2005). This collision also switched off SW-dipping subduction under the GAC (van Benthem et al., 2013). The subduction then migrated to a W-dipping system beneath the present-day Lesser Antilles arc and the Caribbean plate continued to move east to its present day position between North and South America.

Seismic Data
We used data from teleseismic earthquakes recorded by 78 seismic stations located on the North American-Caribbean plate boundary ( Figure 2). By using teleseismic events we are able to attain a good spatial coverage of seismic discontinuities across the plate boundary ( Figure 3). The stations are part of several seismic networks active along the plate boundary between the year 2000 and present day; a full list of network at stations used can be found in Table S1. Data were all accessed from IRIS-DMC using the Obspy-DMT software package (Hosseini & Sigloch, 2017). Events were located at an epicentral distance of 55-85° from the stations . In addition only events with a M W > 5.0 were considered for further analysis. In total, the data set consisted of 8,323 event-station pairs. POSSEE ET AL.

Station Orientation and Rotation
As a pre-processing step in the receiver function analysis, the waveforms were rotated into theoretical P and S components using a free surface transformation (Bostock, 1998;Kennett, 1991). To ensure this is done correctly, we first re-evaluated station orientations. To do this, we used 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). We applied a rotation correction to all stations where the misalignment was greater than the 4σ error determined using bootstrap analysis in DLOPy. A full list of results can be found in Table S1. We also determined the best-fitting near-surface parameters (Vp, Vp/Vs) for the free surface transformation. For the near-surface parameters, we applied a grid-search routine to find the best Vp, Vp/Vs that maximizes the ratio between the S-wave energy on the S component and the P-wave energy on the P component after the rotation. The grid-search space was stacked for all events at a particular station to find the best-fitting parameters for that station. Where a stable result was not achieved, we used the parameters from the nearest station with a stable result. We found that for all stations crustal parameters were between 3.5 and 4.5 km/s for Vp, 1.8-2.0 for Vp/Vs.

S-To-P Receiver Functions
Once rotated into P and S components, we removed all events where the signal-to-noise ratio (SNR) for the S waveform was <3. The SNR was calculated using a noise window prior to the S arrival and a signal window after the S arrival on the S component. This reduced the data set to 2,928 waveforms. S waveforms were then manually windowed before being deconvolved from the P component, using the extended-time multitaper deconvolution method (Helffrich, 2006;Rychert et al., 2012). We applied a zero-phase, fourth order Butterworth filter from 0.03 to 0.25 Hz. Although we tested other bands between 0.02 to 0.5 Hz, we found 0.03-0.25 Hz to be optimal for our data. The deconvolution itself was performed using three Slepian tapers, 50 s in duration with a 75% overlap as suggested by Shibutani et al. (2008). The deconvolutions were then manually assessed and removed from the data set if they showed persistent ringing, large, roughly periodic, similar amplitude signals, from 0 to 70 s in the time reversed records. This process eliminated a further 1,605 waveforms, leaving a total of 1,323 deconvolutions to be used in this study. Deconvolved waveforms were also multiplied by negative one to match the polarity of P-to-S (Ps) receiver functions.
Deconvolved waveforms were then migrated to depth in 3-D by backpropagating the waveform along the raypath of the Sp phase, before being stacked into a grid with 0.5° × 0.5° × 1 km spacing (Angus et al., , 2009Dueker & Sheehan, 1998;Hammond et al., 2011;Rychert et al., 2012). For the crust we used local 1-D velocity models determined in local earthquake studies for individual islands: Puerto Rico (Fischer & McCann, 1984), Hispaniola (Possee et al., 2019), Cuba (Moreno et al., 2002), Jamaica (Wiggins-Grandison, 2004). For the mantle, we assume velocities from IASP91 (Kennett, 1991). We correct for topography, and the depths of the migrated receiver functions are given with respect to sea level.
The grid was then smoothed using the Fresnel zone radius for each waveform, with a minimum width of 50 km. Grid points with fewer than 5 hits were discarded and not considered in the final receiver functions. The number of hits in each bin for conversion depths of 30 and 120 km are shown in Figure 3.
Representative 1-D receiver function profiles for different tectonic regions were extracted from the final 3-D grid by stacking all grid points within a specified region ( Figure 4).

Error Analysis
To estimate uncertainties in our discontinuity depths, we varied the shear wave velocities (Vs) in the models used for migration by ±0.2 km/s, leaving Vp fixed. This number was chosen to be consistent with the uncertainty on the velocity models themselves. The Moho, mid-lithospheric, and LAB phases were shown to vary by ±4, ±5, and ±7 km, respectively, when these changes were applied.
Uncertainties in the amplitude of the receiver functions were estimated using a bootstrap resampling procedure (Efron & Tibshirani, 1986). One hundred bootstrap samples were generated randomly from the data set, where each sample is the same size as the original data set (1,323 waveforms). Each bootstrap sample POSSEE ET AL.  receiver functions, gray shaded regions represent the 95% confidence limits. Depths are given relative to sea level (0 km). Black lines represent the data. Colored dashed lines show the best-fitting synthetic receiver functions, generated by the Vs velocity profile plotted to the right with a Vp/Vs of 1.8. Gray regions on the velocity profile represents the range of models that fit the stacked receiver functions profile within a 95% confidence interval. Red (positive) and blue (negative) shaded regions show where the data is significant from zero. Horizontal red dashed line shows the observed Moho, horizontal blue lines show where a lithosphere-asthenosphere boundary (LAB) phase is observed above error. Cu = Cuba, BaP = Bahamas Platform, NPR = North Puerto Rico, Hi = Hispaniola, SHi = Southern Hispaniola, SPR = Southern Puerto Rico. was migrated to depth using the procedure described above and then statistics were calculated on the bootstrap distribution. We used 2σ to define our 95% confidence interval for the amplitude of our receiver functions. We consider a phase robust where this confidence interval does not straddle zero.

Synthetic Waveform Modeling
We modeled six 1-D receiver function waveform stacks from representative locations using synthetic seismograms created with a reflectivity code (Shearer & Orcutt, 1987). The synthetic seismograms were then deconvolved using the same parameters and method as the data (Helffrich, 2006;Rychert et al., 2012). The ray parameter used for each region represents the mean for the data in each region. The locations chosen were Cuba, Bahamas Platform, Hispaniola, Southern Hispaniola, Northern Puerto Rico and Southern Puerto Rico ( Figure 4). The 1-D profiles also allow us to test which features of the model are robust at a 95% confidence interval.
To obtain the best-fitting velocity structure that could be used to forward model the receiver function, we applied a genetic algorithm (GA). The GA searches for a global minimum in the residual between the observed and modeled receiver function (Goldberg, 1989;Shibutani et al., 1996). We used a normalized Chi-Squared misfit function of the difference between the observed and predicted receiver function. To estimate uncertainty we save all models that fit the observed receiver function within a 95% confidence interval range. The upper and lower bound of these models is then taken to approximate the 95% confidence interval of the error space in our model. To reduce the search space, we only model the 1-D Vs profile and used a Vp/ Vs value of 1.8 to define the Vp model, which is consistent with Vp/Vs observations across the islands from Ps receiver functions (Corbeau et al., 2017). We pre-define the depths of the discontinuities based on the observed phases in the receiver functions. For simplicity, and to avoid overparameterization, we assume sharp discontinuities rather than considering velocity gradients. So, the result is effectively a minimum parameterization. Much work has been done to demonstrate the velocity gradients typically required by receiver functions using data from ocean bottom or ocean island data (Rychert, Harmon, & Armitage, 2018;Rychert, Harmon, & Ebinger, 2014;Rychert et al., 2013).

Significant Discontinuities
Regional Sp receiver function stacks demonstrate which phases are significant above the 95% error from bootstrap analysis (Figure 4). They show a shallow positive phase at 5-12 km that can be seen in all regions, although it is not above error beneath Cuba and Northern Puerto Rico. The largest amplitude positive phase is observed at a depth of 28-34 ± 4 km across all regions, which likely represents the Moho (Figure 4, Table 1). The amplitude of this feature varies considerably from 0.06 to 0.14 ± 0.05, with the lowest amplitude beneath Northern Puerto Rico and the largest beneath Cuba.
We also observe a deeper, positive phase that is significant at a depth of 65 ± 5 km, although it is only sporadically observed. It is only observed above error in the receiver function beneath Southern Hispaniola and Northern Puerto Rico (Figure 4b). The amplitude of this phase is less than the Moho, but is still significant at 0.04-0.06 ± 0.03. A significant negative phase can also be observed at much deeper depths. It is present beneath Cuba and Southern Puerto Rico at a depth of 105 ± 7 km and beneath Southern Hispaniola at a depth of 135 ± 7 km ( Figure 4b, Table 1). The amplitude of this phase ranges from 0.04 ± 0.03 to 0.06 ± 0.04. A negative phase can also be observed in these depth ranges across the other regions of the plate boundary ( Figure 4b). However, the amplitude is much lower at 0.02 ± 0.03 and is not above error. Note. Location of data used to model each region can be found in Figure 4.

Spatial Extent of Discontinuities
The spatial extent of the seismic discontinuities can be assessed using the 3-D depth migrated receiver functions. As was observed in the regionally stacked receiver functions, the Moho depth does not significantly vary spatially, which can be seen both in the N-S cross sections through our final depth migrated model ( Figure  The spatial variability of the strong, positive phase at 65 ± 5 km can also be seen ( Figure 5, profile BB' & CC′). It is present in the southern part of profile BB′, ending abruptly under southern Hispaniola. It also appears to not extend any further south than the Beata Ridge, however, with the caveat that this is approaching the limit of our data coverage. South of the plate boundary, the phase can also be observed under the POSSEE ET AL. Here the phase has an abrupt southern limit with coincides with the location of the subducting North American plate.
The negative phase observed at depths of 95-135 ± 7 km in the 1-D stacks is also seen to be variable across the plate boundary (Figures 7a and 7b). The deeper phase at a depth of 135 ± 7 km is only observed south of the plate boundary (Figure 7a), its western limit appears to be between the Lower and Upper Nicaraguan rise, while the eastern limit is under the Venezuelan basin. We cannot constrain the southern limit of the discontinuity as it extends south beyond the coverage of our land based stations. The shallowest negative phase is observed under the Upper Nicaraguan rise where it is at 95 ± 6 km depth. Elsewhere on the plate boundary, the negative phase is observed at a relatively consistent depth of 105-110 ± 7 km. No significant negative phase is observed north of Puerto Rico under the Atlantic oceanic crust.

Synthetic Waveforms
Modeling synthetic waveforms gives a possible velocity structure required to explain the observed receiver function stacks. Receiver functions are most sensitive to changes in velocity with depth rather than absolute velocity. Therefore, we focus on the magnitude of seismic discontinuities needed to explain the observed positive and negative phases within the receiver function stacks.
All waveforms in the upper crust (<15 km) show low amplitude positive phases. To account for this and for our modeled waveforms to remain within error, seismic layers with increasing velocity with depth need to be introduced within the crust (Figure 4b). Generally the velocity increases with depth fall within the range of 5%-8% and are consistent with local crustal velocity models derived for each region. To model the Moho beneath Cuba, a 28% ± 4% velocity increase is required at a depth of 32 km to produce the best fit. Beneath the Bahamas platform the velocity increase is 18% ± 2% at 31 km, for Northern Puerto Rico it is 9% ± 2% at 34 km, for Southern Hispaniola 18% ± 3% at 28 km, for central Hispaniola 21% ± 4% at 31 km and Southern Puerto Rico 24% ± 3% at 32 km ( Figure 4b).
The significant positive phase observed beneath Southern Hispaniola and Northern Puerto Rico at a depth of 64-66 km can be modeled using a velocity increase with depth with a strength of 11% ± 3% and 8% ± 2%, respectively (Figure 4b). The deeper negative phase, which was observed above error beneath Cuba, Southern Hispaniola and Southern Puerto Rico was modeled using a velocity decrease with depth. Beneath Cuba and Southern Puerto Rico this was a simple 9% ± 3% and 6% ± 2% velocity decrease at 105 and 104 km respectively. Beneath Southern Hispaniola to model the waveform a velocity decrease was modeled over POSSEE ET AL.

Crustal Thickness
The Moho depths observed in this study using Sp receiver functions are in good agreement with previous studies. In Cuba, our mean value of 33 ± 4 km agrees with previous estimates of 20-38 km from Ps receiver functions and is just deeper than the estimate of ∼20-25 km from gravity modeling (Arango-Arias et al., 2014; Bush & Shcherbakova, 1986;González et al., 2012;Toiran, 2003). A notable difference is that we do not observe a large difference in the crustal thickness of southern Cuba (∼20 km) and northern Cuba (>30 km) that has been observed by Ps receiver functions (Bush & Shcherbakova, 1986;Toiran, 2003). Our results only show a gradual increase from 31 km in the south, to 35 km in the north (Figures 5 and 6). While Ps and Sp receiver functions generally both see the sharpest part of the velocity gradient, Sp can be sensitive to broader velocity gradients, given their longer dominant periods (Rychert et al., 2007), which may explain the slight discrepancy between the two methods. The increase in thickness has been attributed to the collision of the GAC with the Bahamas platform, with a purely volcanic arc crust in the south, which has been emplaced over parts of the Bahamas platform in the north to produce the thicker crust during the collision Between Cuba and Jamaica, gravity modeling indicates that the crust should be highly extended continental crust ∼20-25 km thick (Arango-Arias et al., 2014; Ten Brink et al., 2002). While our Moho depths are ∼30 km (Figures 5 and 6), these are given with respect to sea level. Accounting for water depth, which is in excess of 5 km in this region, our estimate of crustal thickness is ∼25 km. This is in excellent agreement with the gravity modeling. The young rifted oceanic crust, 3-7 km thick, commonly associated with the Cayman Trough, only extends 400 km east of the Mid-Cayman Rise spreading center (79°W; Ten Brink et al., 2002). We would therefore not expect to image a discontinuity associated with such a shallow Moho given the location of our stations on Cuba and Jamaica, which is consistent with our results.
The crustal structure of the Bahamas is not well known, with estimates from gravity and refraction lines indicating it ranges from 20 to 30 km thick (Dale, 2013;Dietz et al., 1970;Sheridan, 1972). The origin of this crust is subject to debate, either being thinned continental crust (Dietz et al., 1970;Freeman-Lynde & Ryan, 1987) or thickened oceanic crust (Dale, 2013;Uchupi et al., 1971). Our crustal thickness estimate of 31 ± 4 km is in agreement with recent estimates of ∼30 km from gravity modeling for the southern Bahamas platform (Dale, 2013). We also observe a significant mid-crustal positive seismic discontinuity, an additional velocity increase with depth ( Figure 4). This likely represents the boundary between the carbonate platform, up to 10 km thick (Sheridan, 1972), and either the thinned continental or thickened oceanic crust. From our receiver functions, it is hard to determine which is more favorable, since both would produce a similar positive phase in our results, given that they reflect changes in velocity with depth rather than absolute velocity. However, it was recently suggested that thickened oceanic crust provides a better fit to the gravity data and that the thickened oceanic crust may relate to the Central Atlantic Magmatic Province (CAMP), which formed during the breakup of Pangea at ∼200 Ma (Dale, 2013).
On Jamaica, we find a mean Moho depth of 33 ± 4 km, which shows little variation across the island (Figure 5). Previous crustal thickness estimates on Jamaica are very limited and only local earthquake hypocentral depths have been used to estimate a minimum crustal thickness of 30 km (Wiggins-Grandison, 2004).
Our results are consistent with this, but provide a more robust estimate. However, as observed in Cuba we cannot rule out that there may be finer scale variations in crustal thickness below the sensitivity of this method.
Across Hispaniola crustal thickness has been shown to suggest three tectonic domains, the CLIP in the south, a thickened island arc crust in central Haiti and a forearc terrain in the North (Corbeau et al., 2017). Crustal thickness in the north and south is estimated at 20-30 km from Ps receiver functions, however, the central domain is thickened up to 40 km due to the compressional Trans-Haiti belt and possible under-thrusting of the CLIP under central Haiti (Corbeau et al., 2017). Refraction studies also suggested a thickened crust under central Hispaniola, seeing an increase from 24 to 30 km (Nuñez et al., 2015). The mean Moho depth we observe is 31 ± 4 km, this decreases to 28 ± 4 km across regions of northern Hispaniola and moving offshore to the south. This is very consistent with the refraction survey results (Nuñez et al., 2015). However, we do not observe the 40 km thick crust seen by Ps receiver functions (Corbeau et al., 2017), again potentially because of our broader lateral sensitivity and/or sensitivity to broader velocity gradients in depth.
In Puerto Rico, we observe a mean Moho depth of 31 ± 4 km, which agrees well with previous estimates of 30-40 km using Ps receiver functions (Vanacore et al., 2015). In northern Puerto Rico, the amplitude of the Moho, 0.06 ± 0.03, is significantly lower than in other regions of the plate boundary, with a required velocity increase of 9% ± 2% (Figure 5b). The crustal composition of Puerto Rico is of GAC origin, therefore given the similarity to the other islands, it is unlikely the weak Moho is a compositional effect. In addition, a strong Moho of 0.13 amplitude and velocity increase of 24% ± 3% is seen in Southern Puerto Rico. One reason for the weak Moho in Northern Puerto Rico may be the presence of the subducting North American slab beneath the crust in the north, which has an overall effect of reducing the velocity contrast across the Moho. Alternatively, a layered ocean island crust may also have an effect on the strength and depth of the Moho, which has been hypothesized for receiver function observations beneath the Lesser Antilles arc (Schlaphorst et al., 2018). Vanacore et al. (2015) observe possible evidence for a subducting slab at 60-80 km depth and we also observe a strong velocity increase with depth in this range at 65 km (Figure 4b). The velocity increase may represent the Moho of the subducting plate. In cross-section this discontinuity can be seen to terminate against the top of the subducting slab ( Figure 5, profile CC′). The apparent horizontal nature of this discontinuity is a result of the horizontal smoothing applied to the grid and not the structure of the slab itself, indicated by the observation that phase conversions are only present at grid points directly beneath the subducting slab prior to smoothing being applied.
Under the Caribbean Sea, refraction results image the Moho to be 10-20 km (Granja-Bruña et al., 2010;Mauffret & Leroy, 1997;Mauffret et al., 2001;Núñez et al., 2016). Given that we only use land based stations, we would not expect conversions from a 10-20 km discontinuity from offshore regions. This is consistent with our lack of resolved Moho in all but the closest regions to the island coastlines ( Figure 5).

The LAB
The regional 1-D stacked receiver functions show a negative seismic discontinuity above error at a depth of 105 ± 7 km beneath Cuba, Hispaniola and Southern Puerto Rico (Figure 4b). We interpret this feature to be the LAB, across which a reduction in seismic velocity is typically observed. The depth (105 ± 7 km) also agrees with the predicted LAB depth (100-110 km) for the corresponding age lithosphere (120-140 Myr) in models that relate surface wave velocities to the LAB via empirical relationships based on a thermal model (McKenzie et al., 2005;Nerlich et al., 2014). Thermal models would only predict a gradational reduction in temperature from the lithosphere to the asthenosphere, which experimental results suggest would translate to a similarly gradual reduction in seismic velocities, over 10s of kilometers (Jackson & Faul, 2010;Rychert et al., 2020;Tharimena et al., 2017). However, our LAB negative phases can only be modeled by a sharp velocity reduction of 6-9% ± 2% (Figure 4b). This may imply that while the LAB is predominantly thermally controlled across much of the plate boundary, there is also likely a change in properties which serves to enhance the strength of the observed LAB converted phase. One such mechanism could be the ponding of melt at the base of the lithosphere (Kawakatsu et al., 2009;Naif et al., 2013;Schmerr, 2012). This could have occurred through partial melting of the asthenosphere (Kawakatsu et al., 2009), or through the generation of melt during subduction beneath the GAC, which occurred from 140 Ma through to 55 Ma (J. L. Pindell, 1990).
Previous estimates of lithospheric thickness in Cuba range from 65 to 93 km across the island from the joint inversion of Raleigh waves and Ps receiver functions (González et al., 2012). Our estimate of 105 ± 7 km is deeper, however, González et al. (2012) observed gradual velocity reductions over depth ranges of 60-110 km, which are consistent with our results. We also observe a lower amplitude negative discontinuity at ∼60 km beneath Cuba ( Figure 5, profile AA′), which may relate to the depth at which the LAB has previously been interpreted. Ps receiver functions also have difficulty imaging deep discontinuities, since they can be masked by interference from crustal reverberations (González et al., 2012).
Beneath the Bahamas Platform we only observe a weak negative discontinuity at ∼115 km, which is not above error. This may reflect gradational or a subtle velocity change across the LAB that the receiver functions cannot resolve. This may be the same across the incoming North American plate, however, much of the plate is also outside of our resolvable region, even with a deep discontinuity such as the LAB (Figure 3).
To the south west of Jamaica, beneath the Upper Nicaraguan Rise, we were unable to attain a representative 1-D receiver function due to a lack of information at certain depth intervals. However, we do observe a negative seismic discontinuity at 90-95 km depth (Figure 7b). We interpret this to represent the LAB beneath the Chortìs block, which is part of a Precambrian Craton that forms much of Central America and extends offshore into the Caribbean Sea. The eastern limit of this discontinuity is along the Pedro Bank fault zone (PBFZ; Figure 7b), which confirms the hypothesis that this fault zone marks the boundary between the Chortìs block and the CLIP (Mauffret et al., 2001;Sanchez et al., 2019). The PBFZ, which has accumulated a significant lateral offset during the eastward migration of the Caribbean plate (Boschman et al., 2014;Sanchez et al., 2019), may also explain the relatively abrupt transition from thin lithosphere in the west to thicker lithosphere beneath the CLIP. Sp imaging can resolve sharp changes in depth, but not necessarily the slope of the change (Lekić & Fischer, 2017;Lekic et al., 2011).

A Dehydrated Mantle Root Beneath the CLIP
Beneath much of the Caribbean Sea and southern Hispaniola we observe a deep negative seismic discontinuity at 135 ± 7 km, which corresponds to a velocity decrease of 7% ± 3% at this depth. This is ∼30 km deeper than surrounding lithosphere of the same origin and age (120-140 Ma;Nerlich et al., 2014). We suggest this deep discontinuity is still the LAB, but propose that its deeper depth is related to its association with the formation of the CLIP.
The formation of the CLIP is thought to have been due to plume related activities associated with the Galápagos hot spot ranging in age from 140 Ma through to 60 Ma, with a peak event at ∼89 Ma (Hoernle et al., 2002(Hoernle et al., , 2004Révillon et al., 2000;Whattam & Stern, 2015). As a result, we interpret the seismic discontinuity to represent the base of the residual dehydrated root of the depleted mantle, whereby the discontinuity depth represents the base of melting during peak temperatures associated with the main event at ∼89 Ma. The extraction of melt from this layer leaves it more viscous than the hydrated mantle beneath, and thus means it is seismically fast and also relatively rigid (Yamamoto & Morgan, 2009). However, since compositional changes are not typically sufficient to explain strong sharp seismic discontinuities, some other mechanism is likely required. For instance, partial melt beneath the compositionally depleted root could further decrease velocities beneath the CLIP.
A thick layer of dehydrated and depleted mantle could be caused by melt extraction owing to a previous large melting event such as a past ridge or hot-spot volcanism (Gaherty et al., 1996;Morgan, 1997;Yamamoto & Morgan, 2009). An expectation of this layer is that it should thicken as potential mantle temperatures increase, since the elevated temperatures will increase the depth at which melt extraction and dehydration of the mantle occurs. This is consistent with what has been observed at modern day hot spots such as the Galápagos (Byrnes et al., 2015;Villagómez et al., 2014;Rychert, Harmon, & Ebinger, 2014) and Iceland (Rychert, Harmon, & Armitage, 2018), where the increased thickness of the dehydrated mantle layer has been linked to the elevated potential mantle temperatures. Dehydrated roots have also been interpreted beneath ancient large igneous provinces (LIP) such as the Ontong Java Plateau (OJP; Tharimena et al., 2017), suggesting the residual dehydrated mantle is stable over periods of time much longer than the thermal anomaly. Globally, sub-hotspot lithospheric thicknesses range from 50 to 90 km, without necessarily showing a trend with age or thinning owing to plume-related heating .
Assuming the base of this layer marks the depth of the anhydrous solidus of the mantle during the peak melting event for the CLIP, it can be used to estimate the peak potential mantle temperature during the formation of the CLIP. Using the interpretation that the anhydrous solidus was at 135 ± 7 km and solidus curves from Herzberg et al. (2000), this gives a potential mantle temperature of 1585°C ± 20°C. This is a remarkable agreement with the 1560°C-1620°C estimated by Herzberg and Gazel (2009) using petrological evidence from CLIP material. The modern day Galápagos hot spot is thought to be in the range 1465°C-1500°C from similar lines of evidence (Byrnes et al., 2015;Herzberg & Gazel, 2009;Rychert, Harmon, & Ebinger, 2014), therefore the plume head of the Galápagos hot spot responsible for the formation of the CLIP had an excess temperature of 60°C-165°C compared to the present day hot spot tail feeding the Galápagos. This is in general agreement with the notion that a massive, and therefore hotter, plume head arrived at Earth's surface, forming the CLIP, whereas today's hotspot reflects relatively cooler temperatures associated with a thinner plume stem.

A Mid-Lithospheric Discontinuity
Across parts of the plate boundary, we observe an intermittent seismic discontinuity at 65 ± 5 km, which is predominantly positive, but varies in strength and polarity (Figures 5 and 6). It is particularly evident as a strong positive phase offshore southern Hispaniola where the lithosphere is of CLIP origin and also offshore of northern Puerto Rico ( Figure 5). As discussed above, the velocity increase with depth observed beneath northern Puerto Rico can be explained by the presence of the subducting North American slab at this depth (Vanacore et al., 2015). However, under the CLIP to the south there is no such subducting slab to explain such a feature. The Muertos Trough, which has been suggested to be an early subduction zone (Byrne et al., 1985), has not experienced enough subduction for a slab to be present at 60 km depth (Granja-Bruña et al., 2010;Harris et al., 2018). In addition, more recently it has been suggested that the Muertos Trough is a retro-arc thrust fault associated with singular subduction of the North American plate (Ten Brink et al., 2009). In this section, we therefore postulate as to what might cause a significant seismic discontinuity at this depth beneath the CLIP.
There have been many mechanisms proposed to explain a mid-lithospheric discontinuity, which have included partial melt, compositional layering, elastically accommodated grain-boundary sliding (EABGS) and anisotropy (Karato et al., 2015). Of these we exclude EABGS from possible interpretations as this mechanism should only produce a 5% velocity decrease with depth (Karato, 2012), which is at odds with our observation of a predominant velocity increase with depth that is as large as 11% ± 3%. We also exclude partial melt, as volcanism beneath the Greater Antilles and the CLIP has been extinct for 50 myr. A velocity increase with depth could only indicate the base of a melt rich layer and there is no evidence for a shallower low-velocity layer, which indicates the melt rich-layer itself.
Our observations of a seismic discontinuity at a depth of 65 km beneath the CLIP are similar to those observed beneath other LIPs, such as the OJP, at depths of 60-80 km (Tharimena et al., 2017;Tonegawa et al., 2019). While Tonegawa et al. (2019) interpreted a velocity increase with depth using Ps receiver functions and Tharimena et al. (2017) interpreted a velocity decrease using SS waveform modeling, both attributed the discontinuity to compositional variations within the lithosphere. In particular, Tonegawa et al. (2019) suggested the phase transition from Spinel to Garnet, which represents a velocity increase with depth (Revenaugh & Jordan, 1991), may explain their observations. This mechanism may be able to reconcile our observations of a velocity increase at 65 km, however, a continuous compositional change across the CLIP does not explain a discontinuity that is variable in strength and sometimes polarity. To explain our observations any compositional change would need to be spatially variable and constraints from mantle xenoliths from CLIP material are lacking to test this hypothesis.
Anisotropy is another mechanism by which a seismic discontinuity can be introduced. The Greater Antilles region is known to show azimuthal anisotropy, which has been observed through shear wave splitting (Hodges & Miller, 2015;Possee et al., 2020). A ∼25 km thick anisotropic layer in the upper mantle with a fast direction parallel to relative plate motion is suggested from shear wave splitting results on Hispaniola (Possee et al., 2020). Azimuthal anisotropy adds complexity to Sp receiver functions (Farra & Vinnik, 2000). The polarity of scattered phases in the P-SV system varies depending on the back azimuth of the raypath and the polarization of the wave relative to the anisotropic medium, which has been demonstrated for underside reflections interacting with a hexagonal anisotropy approximation (Rychert, Harmon, & Schmerr, 2014). Globally, for thick continental lithosphere this has been ruled out since the discontinuity is uniformly a velocity decrease (Fischer et al., 2010;Karato et al., 2015). However, beneath the CLIP we observe a velocity increase across the mid-lithospheric discontinuity, with a strength that is not uniform and is negative in one region (68°W, Figure 7, profile AA′). This makes a change in azimuthal anisotropy with depth likely a good explanation for our result. However, modeling the exact anisotropic structure required by the data is complicated by both increased parameter space including olivine crystallographic orientations that could be variable in space and in two layers and also the source polarization and back azimuth of each earthquake. In this study, we do not have enough data to fully constrain such a parameter space, and therefore any model produced would be highly non-unique.

Plate Boundary Structure
In this section, we discuss the implication of this work for the structure and evolution of the North American-Caribbean plate boundary. Present day oblique subduction of the North American plate under Puerto Rico and through to eastern Hispaniola is well known and has been geophysically imaged (Harris et al., 2018;van Benthem et al., 2013;Vanacore et al., 2015). Our observation of a subducting slab at ∼65 km beneath northern Puerto Rico is therefore consistent with this. The fate of the slab westward of central Hispaniola is however, poorly understood. Tomographic models suggest the edge of the slab is under central Hispaniola (Harris et al., 2018;van Benthem et al., 2013), however, this longitude also corresponds to a dramatic reduction in station density and often marks the end of model resolution. Other studies have suggested the possible presence of a subducting slab or remnant slab beneath western Hispaniola (Corbeau et al., 2019) and Cuba (González et al., 2012). In our study, we see no evidence for a seismic discontinuity that can be attributed to a subducting or remnant slab westward of eastern Hispaniola.
The observation of a deeper LAB under CLIP material also allows us to place spatial constraints on the CLIP across the northern Caribbean plate and plate boundary. The western limit of the CLIP maps spatially to the Pedro Bank fault zone, which separates the Lower and upper Nicaraguan rise. This also makes this the likely boundary between the Chortìs block and the CLIP, which is consistent with previous hypotheses (Mauffret et al., 2001;Sanchez et al., 2019). Our results show that the northern limit of the deep LAB associated with the CLIP, extends north beyond the surface expression of the EPGF (Figure 7b). However, we still suggest that the EPGF marks the northern limit of the CLIP. A possible and simple explanation for our observation is that the EPGF has a northerly dip along parts of the plate boundary, allowing the CLIP LAB to be further north at depth. On southern Hispaniola, CLIP material has been observed outcropping on Haiti's southern peninsula and Massif de la Selle (Mann et al., 1995). However, Corbeau et al. (2017) also suggested that the CLIP may have been underthrust beneath central Hispaniola. Instead, our results of both the deep LAB and mid-lithospheric discontinuity, which we associate with the CLIP, suggest that the CLIP terminates against the Cul-de-Sac-Enriquillo (CSE) basin, as it not present at depth beneath central Hispaniola.
In the past, the compressional tectonics of Hispaniola have been solely attributed to the collision of the GAC with the Bahamas platform (Mann et al., 1995;Pubellier et al., 2000). Here we suggest that the thickened lithosphere of the CLIP may have contributed to the active compressional tectonics observed in southern Hispaniola and also to the current plate boundary segmentation. This is consistent with the active south dipping thrust faulting observed in southern Hispaniola, that sees CLIP material being actively thrust over central Hispaniola at the southern border of the CSE basin (Possee et al., 2019;Rodriguez et al., 2018;Symithe & Calais, 2016).
In our study, we observe a NW-SE eastern boundary to the CLIP starting at the western end of the Muertos Trough, which we infer from the deeper LAB beneath the CLIP. This is consistent with the location of a transition from thickened crust associated with the CLIP to normal oceanic crust within the Venezuelan basin (Mauffret et al., 2001). We are unable to resolve this boundary further to the south or east due to the limitations of the land based seismic station locations used in this study. However, this result confirms that the CLIP is not uniformly present under the Caribbean Sea and Venezuelan basin all the way to the Aves Ridge. To map the full spatial extent of this deep discontinuity will likely require future deployments of ocean bottom seismometers in the Caribbean Sea.

Conclusions
In this study, we have observed both crustal and mantle seismic discontinuities using Sp receiver functions across the North American-Caribbean plate boundary. To our knowledge, this is the first study to document the seismic discontinuity structure of the lithosphere at a regional scale across the plate boundary. We observe positive seismic discontinuities at consistent depths of 28-34 km beneath Puerto Rico, Hispaniola, Jamaica, Cuba and the Bahamas platform. These are interpreted to be the Moho and are consistent with previous estimates of crustal thickness across the GAC and the Bahamas platform.
In northern Puerto Rico, we also observe a positive seismic discontinuity at 65 ± 5 km, which is interpreted to be the top of the subducting North American slab. No similar discontinuity is seen to the west beneath central Hispaniola or Cuba to suggest a possible subducting or remnant slab in these locations. A seismic discontinuity at similar depths is also observed beneath parts of the CLIP in Southern Hispaniola and the Caribbean Sea. The intermittent nature of the discontinuity and its variable polarity may be mid-lithospheric, potentially caused by anisotropic effects, although further testing is required.
Across much of the plate boundary, the LAB is between 90 and 115 km, observed as a negative seismic discontinuity. The exception is beneath regions of the CLIP that show thickened crust, where we observe the negative discontinuity at a depth of 135 ± 7 km. We have interpreted this to be the base of a residual dehydrated mantle root produced during peak temperatures of the main CLIP forming event (∼89 Ma) as the Caribbean plate moved over the Galápagos hot spot. The depth of this discontinuity indicates potential mantle temperatures during the main event forming the CLIP were up to 1585°C ± 20°C, which is ∼100°C hotter than the present day hot spot. We also suggest this thicker lithosphere may contribute to the active compressional tectonics observed in southern Hispaniola and plate boundary segmentation.  (Wessel & Smith, 1998). We thank Alan Levander and the second anonymous reviewer for their constructive comments, which have helped improve this manuscript.