Widespread Hydration of the Back Arc and the Link to Variable Hydration of the Incoming Plate in the Lesser Antilles From Rayleigh Wave Imaging

Subduction zone dynamics are important for a better understanding of natural hazards, plate tectonics, and the evolution of the planet. Despite this, the factors dictating the location and style of volcanism are not well‐known. Here we present Rayleigh Wave imaging of the Lesser Antilles subduction zone using the ocean bottom and land seismic data collected as a part of the VoiLA experiment. This region is an important global endmember that represents a slow (<19 mm/yr) convergence rate of old (80–120 Ma), Atlantic lithosphere formed at a slow spreading ridge. We image the fast slab, the fast‐overriding plate and the slow mantle wedge across the entire arc. We find slow velocity anomalies (∼4.1 km/s) in the mantle wedge directly beneath the arc with local minima beneath Dominica/Martinique, Montserrat and the Grenadines. We observe that slow velocities in the wedge extend 200 km into the back arc west of Martinique. The slowest mantle wedge velocity anomaly is more muted than several global wedges, likely reflecting the lower temperatures and less partial melt predicted for the Antilles. Subducted fracture zones and plate boundaries are a potential source of hydration, since they are located near the anomalies, although not directly beneath them. To match our observations, geodynamic models with a broadly hydrated mantle wedge are required, which can be achieved via deep hydration of the slab, and fluid release further into the back arc. In addition, 3‐D flow and melt migration or ponding are required to explain the shape and location of our anomalies.

Frueh-Green, 2010), via fault systems such as transform faults/fracture zones, or at normal or detachment faulting (White et al., 1984). Further alteration may occur at the outer rise as the plate bends before descending into the subduction zone, reactivating older fault systems or creating new ones (Fujie et al., 2013;Ranero et al., 2003;Van Avendonk et al., 2011). At depth, once the fluid is released it may interact with the mantle wedge and flow through the system in a variety of ways. In the simplest models, fluids rise vertically from the slab due to their buoyancy, creating melt in the hotter mantle wedge above and melt flux at the surface is determined by the relative rates of melt generation and melt transport (e.g., Harmon & Blackman, 2010;Morgan, 1987). However, the solid flow in the mantle wedge can affect the fluid pathway and cause nonvertical transport, with the fluids being dragged down with the slab toward the back arc if the permeability is relatively low, which may occur, for example, in regions of small mantle grain sizes associated with relatively low temperatures or high stress or strain rates (Cagnioncle et al., 2007;Cerpa et al., 2017;Wada & Behn, 2015). This can have the effect of hydrating the mantle wedge well into the back arc as well as supplying melt to the arc itself. Compaction can create fluid and melt channels along the slab interface and beneath the upper plate, where the solidus acts as a permeability boundary, that can focus melt toward the arc (England & Katz, 2010;Ha et al., 2020;Wilson et al., 2014). In addition, weak sediments, fluids, and melt near the slab interface can create a low viscosity zone that can lead to Rayleigh-Taylor instabilities in the solid flow, causing rapid diapiric upwelling of material (Behn et al., 2011;Gerya & Yuen, 2003). These different processes would have different seismic structures associated with them, allowing us to examine which processes might be in play using seismic methods.
Over the past few decades our observational understanding has increased substantially. For instance, active source seismic results reveal alteration of downgoing plate (Davy et al., 2020;Ranero et al., 2003;Van Avendonk et al., 2011). Seismic waves image the mantle wedge, upper plate, and downgoing slab (Eberhart-Phillips & Reyners, 2012;Eberhart-Philips et al., 2020;Rychert et al., 2008;Syracuse et al., 2008;Wiens et al., 2006). Geochemistry reveals the mantle conditions that result in surface volcanism (Kelley & Cottrell, 2009;Plank et al., 2013). Geodynamic modeling allows us to test hypotheses regarding the controlling factors of the system (Cerpa et al., 2017;Davies et al., 2016;Gerya et al., 2006;Hacker, 2008;van Keken et al., 2011). In particular, the circum-Pacific subduction zones, such as New Zealand, Japan, Chile, Mexico, Central America, Marianas, and Tonga have been well studied. However, there remain a few key locations, such as the Lesser Antilles, that have yet to be studied at a scale sufficient to understand this system.
The Lesser Antilles is a global end member of subduction in terms of slow convergence of an old subducting plate (Figure 1). At present, the North and South American Plates are subducting beneath the Caribbean Plate at 19 mm/yr in a west-south-westerly direction (Kreemer et al., 2014), which is in stark contrast to convergence rates of old plates in the Pacific which can be > 70 mm/yr. Over the past 20 Myr or so this convergence rate has not varied much (Matthews et al., 2016) and the geometry of the system has been relatively stable (Boschman et al., 2014) as the Caribbean plate has continued to move relative to North and South America. Slow convergence is important in terms of our understanding of subduction zones because the slab is expected to warm up faster with depth due to conductive cooling, while the mantle wedge is expected to be cooler due to sluggish convection and cooling from the slab (Hacker, 2008;van Keken et al., 2011). Part of the water released from dehydration reactions is expected to occur at relatively shallow depths (Schmidt & Poli, 1998), possibly serpentinizing the cold nose (Abers et al., 2017;Pozgay et al., 2009), and melting in the mantle wedge should be strongly controlled by flux melting due to the volatiles as opposed to vigorous corner flow and decompression melting (Stern, 2002).
The Lesser Antilles is also a global end member for subduction of slow spread oceanic lithosphere. The downgoing lithosphere was generated at slow spreading centers, with the North American portion of the slab primarily being generated at the slow spreading Mid-Atlantic Ridge, while in the south, the plate is likely comprised of proto-Caribbean lithosphere which also was generated at a now extinct slow spreading center ( Figure 1) with a more N-S spreading direction (Cooper et al., 2020;Muller et al., 2019), which accommodated the separation of North and South America. There is also a difference in slab age with the Mid-Atlantic lithosphere being ∼80-85 Myr old at the trench and the proto-Caribbean lithosphere being ∼155 Myr old near the trench and likely >120 Myr in the subducted portion ( Figure 1) (Seton et al. 2020). Slow spread crust and oceanic lithosphere may be more altered/serpentinized. It has been suggested that up to 70% of the Mid Atlantic crust may be made up of detachment fault type ocean crust (Blackman et al., 1998;Escartin et al., 2008), where core complexes and irregular abyssal hill fabric expose lower crust and mantle rocks to sea water allowing them to be serpentinized, which may also enhance the potential for the water storage (Bach & Frueh-Green, 2010). In the lithosphere going into the Lesser Antilles there is some evidence that core complexes may be more altered at shallow levels, while in regions of abyssal hill fabric there may be a deeper alteration (Davy et al., 2020). There are also fracture zones on the incoming/ HARMON ET AL.   Becker et al. (2009). Seismic stations are shown as inverted red triangles, white, black, cyan and purple lines indicate tectonic features: thrust front of the Lesser and Greater Antilles subduction, white dashed lines indicate traces of major fracture zones (FZ) on the North American Plate generated at the Mid-Atlantic Ridge and black dashed lines indicate fracture zones on the proto-Caribbean lithosphere from Cooper et al., (2020). The North-South American plate boundary is a diffuse region roughly located between the Vema and 15-20 Fracture Zones and perhaps farther north (Bird, 2003). Purple dashed line indicates the inferred proto-Caribbean lithosphere boundary on the South American plate from Cooper et al. (2020) and cyan line extinct proto-Caribbean spreading ridge. Small black circles extending westward from the white dashed lines indicate the projection of fracture zones onto the subducted lithosphere after Harmon et al. (2019). Moment tensors for earthquakes located within the downgoing plate in the fore-arc are also shown. Island names are indicated: Saba (SB), Saint Kitts & Nevis (SK), Montserrat (MT), Guadeloupe (GD), Dominica (DM), Martinique (MQ), Saint Lucia (SL), Saint Vincent (SV), Grenada (GR), Barbados (BB), Tobago (TB) and Trinidad (TR) as well as the Aves Ridge (AR). downgoing plate that may play a role in hydrating the slab (Figure 1). At the outer rise, further alteration may occur as faults, which are generally more prevalent in lithosphere formed at slow spreading centers, are reactivated, in particular if the faults are favorably oriented, striking parallel to the trench. However, the sparsity of normal or reverse focal mechanisms suggests that the reactivation of ridge related normal faults might play less of a role in hydrating the slab. In the Lesser Antilles, the global centroid moment tensor catalog indicates a number of earthquakes with strike-slip mechanisms located between the Vema and Mercurious fracture zones, in other words different from normal faulting typically associated with the outer rise. In addition, the proto-Caribbean lithosphere and its transition to the South American continental margin may be more akin to amagmatic continent-ocean transitions such as the Galicia margin. In the Galicia margin, continental faulted crustal blocks gradually transition to fault dominated oceanic crust and, in some cases, exposed peridotites at the surface (Bayrakci et al., 2016) over a distance of tens to hundreds of km. The continental material may not be as tectonized with the alteration from hydrothermal circulation more focused on individual fault strands.
The Lesser Antilles has strong along strike variation in forearc sedimentation (Whittaker et al., 2013), arc chemistry (Macdonald et al., 2000) and seismicity (Bie et al., 2020;Wadge & Shepherd, 1984). Sedimentation rates in the fore and back arc are high in the southern part of the subduction zone due to outflow of the Orinoco river delta, which produces a > 10 km sediment package in the south obscuring the trench of the subduction system (Chichester, Rychert, Harmon, Allen, et al., 2020;Whittaker et al., 2013). The seismicity rates are higher in the north of the arc (Wadge & Shepherd, 1984). Higher levels of low magnitude earthquakes, or higher b-values, are located in the regions where long offset fracture zones have been subducted (Schlaphorst et al., 2016) (Figure 1), which may be caused by fluid release from the fracture zones. In addition, the δ 11 B record of the erupted lavas along arc has a peak centered beneath Dominica in the central part of the arc, suggesting an association with subducted fracture zones there (Cooper et al., 2020). Crustal production in the Lesser Antilles is low, 3 × 10 3 km 3 , in comparison to those in the circum-Pacific arcs, which are an order of magnitude greater (Wadge, 1984). The greatest crustal production in the Lesser Antilles over the past 100,000 years occurs in this central region of the arc, with Dominica having the greatest production of all the Islands (Wadge, 1984). Finally, the contributions to the primary mantle melts from fluids from the slab and sediments vary along strike, with sediment components, e.g. Ba/La,Nd isotope ratios decreasing and fluid components, e.g. Sr/Th ratios, increasing northward along the arc (Macdonald et al., 2000).
To understand the relationships between the along arc variations and arc volcanism, we investigate the mantle structure of the Lesser Antilles arc using Rayleigh wave tomography to produce a 3D shear velocity model across the entire subduction zone, from the incoming plate to the backarc, using both the land and ocean bottom seismic data from the Volatiles in the Lesser Antilles (VoiLA) project . We compare our seismic results to predictions of seismic velocity from a suite of geodynamic models to understand where and how fluids might be released into the mantle wedge to explain our results.

Seismic Data Processing
For our Rayleigh wave tomography, we use the vertical component of the seismograms from the land and ocean bottom stations shown in Figure 1 for both teleseismic earthquakes and ambient noise cross correlation. All data have the mean, trend and instrument response removed and the records are decimated to 1 Hz. For the earthquake data, a time window of 4 h after the origin time is selected and dispersion is measured using frequency time analysis (Landisman et al., 1969;Levshin & Ritzwoller, 2001) of the fundamental mode Rayleigh wave. The data are visually inspected for the good signal-to-noise ratio and each event was checked for coherence across the array using a beamforming method. Incoherent events and stations with poor signal to noise ratios (<3) were removed from further analysis. Up to 2,486 dispersion measurements were used at a given period from 180 events from 20 to 111 s period in the teleseismic data set.
Ambient noise data for each station and for each day were cross correlated and stacked for a 1-year period following the method of Bensen et al. (Bensen et al., 2007) to generate the Noise Correlation Functions (NCF). Each NCF was visually inspected for good signal to noise and frequency time analysis (Landisman et al., 1969;Levshin & Ritzwoller, 2001) was used to isolate the fundamental mode of the Rayleigh wave and measure the dispersion. We use up to 445 NCF station-to-station pairs in our further analysis from 20 to 45 s period. Figure 2 shows example waveforms from the 2016 Christmas Day earthquake event in Chile as well as NCF as a function of interstation distance.

Rayleigh Wave Tomography Method
The Rayleigh wave tomography we use is a two-step process to produce a 3-D shear velocity model across the region. In the first step the amplitude and phase of the teleseismic dispersion and the phase of the ambient noise is inverted at each period to generate phase velocity maps. In the second step, the 3-D shear velocity volume is created by inverting for a 1-D shear velocity profile at each point in the phase velocity maps using all periods.
For the phase velocity map inversion, we use the two plane wave method (Forsyth & Li, 2005;Yang & Forsyth, 2006), which accounts for off-great circle propagation and finite frequency effects in teleseismic Rayleigh waves, modified to also invert NCF dispersion (Harmon & Rychert, 2016) in the period range where the two methods overlap. At the longer periods (>50 s) where there is no NCF data, we stabilize the inversion using the anomaly structure from the next lowest period as the starting model. Otherwise, we use the average phase velocity as the starting model. We assume an a priori standard error for the phase velocities of 0.2 km/s in the inversion. The inversions are parameterized on a set of nodes shown in Figure 3, which are spaced at 0.5° throughout the central part of the study region, with a more broadly spaced set of nodes 1.0° on the exterior where there is a greater a priori standard error. These outer nodes allow heterogeneity in earth structure outside the study region to be absorbed. The finite frequency kernels at 0.1° resolution are averaged over an 80 km-radius and interpolated onto the nodal structure for inversion. After the inversion, we then resample at 0.1° using a weighted distance interpolation with an 80 km-radius at each node to attain the phase velocity maps presented here. For greater detail on the two-plane wave method we refer the reader to the original work (Forsyth & Li, 2005;Yang & Forsyth, 2006). Raypaths for both the teleseismic and NCF data used at select periods are shown in Figure 3.
For the shear velocity inversion, we determine the best fitting regional average shear velocity model using an iterative damped least squares approach (Tarantola & Valette, 1982), assuming an a priori standard error of the model of 0.2 km/s. We parameterize the shear velocity model, as consisting of a variable thickness water layer (where applicable), sediment layer and crustal layer, followed by 5 km thick mantle layers down to a minimum of 400 km depth. For the 1-D average model, we use the average water depth, and a 2-layer crust, with a total initial thickness of 20 km. The model becomes the basis for our starting model. There are strong lateral variations in crustal and mantle structure across the region, which require the use of a HARMON ET AL.   priori knowledge to help seed the shear velocity inversion, as it is a very nonlinear problem. We use the water layer thickness from satellite altimetry derived bathymetry in the region (Becker et al., 2009), with the initial sediment thickness taken from the global sediment thickness map and local estimates of sediment from receiver functions (Chichester, Rychert, Harmon, Allen, et al., 2020;Whittaker et al., 2013), with shear velocities that scale linearly with thickness from 0.5 km/s to a maximum of 1.5 km/s. Crustal thickness in the starting model is determined by assuming Airy isostasy, which produces a maximum crustal thickness of 35 km, and we impose a minimum crustal thickness of 6 km. We then allow the inversion to adjust the crustal and sediment thickness and velocity to satisfy the data. The mantle structure is set to 4.5 km/s in the slab and upper plate. The upper plate is assumed to be 60-70 km thick based on a companion S-to-P receiver function study in the region , and we assume a low velocity zone of 4.3 km/s in the mantle wedge. We used two different models for the depth to the slab, Bie et al. (2020) and Slab2 (Hayes et al., 2018), generating two different inversion results. Figures S1-S4 show the starting model shear velocity structure and crustal structure for both slab geometries. We assume a Vp/Vs ratio of 1.78. We allow up to 18 iterations or terminate the inversion if the normalized chi-squared value drops below 1.

Checkerboard Tests
We estimate the resolution of our 2-D phase velocity model using checkerboard tests. We insert anomalies of ±0.05 km/s and ±0.10 km/s relative to the average velocity with 1 × 1° length scale and 2 × 2° length scale, respectively as our input velocity structure. We examine the recovery of this structure at 22, 45 and 100 s period, which spans the range of periods used in this study. We use the same ray paths shown in Figure 3 for the inversions including both ambient noise and teleseismic data when both are available and add a small percentage of noise to the synthetic data, which achieves a similar error in the checkerboard inversion to that in the observed data.
Our recovery test results are shown in Figure 4, where we show the recovery of the 1 × 1° length scale for 22 and 45 s period and 2 × 2° length scale for the 100 s period tests. Overall, at 22 and 45 s the recovery of the anomaly pattern is good, especially within the array. In the western portion of the model, there is an evidence for E-W smearing of anomalies. The magnitude of the anomaly at these periods is generally well recovered, within 2% in most cases, although the amplitude of the recovered anomalies in the forearc is lower than the input. At 100 s the recovery of the 1 × 1° length scale anomalies is poor so we do not show it. The 2 × 2° length scale recovery at 100 s is good, with the pattern well recovered within the array. Again, there is evidence for E-W smearing in the western most part of the model. The magnitude of the anomalies at 100 s is within 3% of the input model, but again the strength of the anomalies in the forearc in the recovery test is weaker than the input by up to 3%. We also tested whether or not high velocities beneath the arc would be resolved in a separate set of tests ( Figure S5), and found similar results. Overall, these results suggest that anomalies with lateral extents of ∼100 km or larger are well imaged within the array in the shallowest parts based on the sensitivity of Rayleigh waves at 22-45 s (<100 km depth sensitivity). The longer period data, which have peak sensitivity at >100 km depth, can resolve structures at 200 km length scales or longer.

Geodynamic Modeling
We generate a set of 2-D kinematically driven thermal models of the Lesser Antilles Arc. We model temperature and mantle flow velocity by solving the coupled incompressible Stokes and energy equations in the Boussinesq approximation using the finite-element, control-volume Fluidity framework (Davies et al., 2011). Full details of the methods as well as the reproduction of kinematic subduction benchmarks can be found in Le Voci et al. (2014), Davies et al. (2016) and Perrin et al. (2018). Our model domain is between 500 and 800 km length (depending on the cross section) by 200 km deep and aligned in the direction of plate convergence. The slab geometry is prescribed, based on the seismicity relocated by Bie et al. (2020). The top 10 km of the slab is prescribed to move at 20 mm/yr, consistent with velocities for the North American-Caribbean plate convergence projected into the plane of our model. A zero velocity is prescribed in the upper plate down to 50 km, as well as in a 5-km thin layer on top of the slab down to the coupling transition depth, below which the slab couples to the convecting mantle wedge. Side boundaries and the bottom boundary below the subducting plate are stress free, while velocities equal to the subducting plate velocity are prescribed at the bottom boundary of the mantle wedge to mimic the effect of the deeper slab (Le Voci et al., 2014). The top surface of the model is held at a constant temperature of 0°C, the bottom and overriding plate side of the model are thermally insulating. The thermal structure of an 80 Myr plate from a half space cooling model with a potential mantle temperature of 1350°C is imposed on the slab side of the model, and a time-dependent half-space cooling model is imposed on the upper-plate side. In the mantle surrounding the prescribed slab, a composite dislocation-diffusion creep rheology is used, with parameters as given in Perrin et al. (2018). For the computation of viscosity and for the conversion to a seismic structure, an adiabatic gradient of 0.5°C/km is added to the Boussinesq temperature field. The models are run for 60 Myr. This is sufficiently long for the subducting plate to achieve a steady-state temperature, and by HARMON ET AL.  this time the upper plate reaches a thickness consistent with the receiver function results from  and heat flow from Manga et al. (2012). We also investigate the impact of a thinner upper plate, equivalent to an oceanic lithospheric age of 20-40 Myr by analyzing earlier stages of the model.
We convert the thermal models to seismic velocity after Goes et al. (2012). Crustal velocities for the upper plate are set according to the Antilles 1-D model from Bie et al. (2020). For the rest of the domain, we use the thermodynamic data base from Holland and Powell (1998) implemented in the Gibbs Free Energy minimization code PerPleX (Connolly, 2005) to compute the phase composition as a function of pressure and temperature and the elastic parameter compilation from Abers and Hacker (2016) to obtain the corresponding seismic velocities and density (e.g., Eeken et al., 2018). We consider a peridotitic composition for most of the domain, except for the top of the slab where we test either a hydrous or dry basalt crustal composition overlying a 14 km thick layer of hydrous or dry harzburgitic mantle lithosphere (Halpaap et al., 2019). This approach also provides estimates of the slab dehydration depths (as in Hacker, 2008). The effects of frequency-dependent anelastic attenuation are added either using the model from Jackson and Faul (2010) either in its original dry form (Case A), or with a correction for hydration following Behn et al. (2009) with 5000 H/10 6 Si water concentration in the mantle wedge (all other Cases B-D). It is debated whether volatiles have any direct effect on mantle wedge attenuation (Abers et al., 2014;Cline et al., 2018) or whether the low velocities and Q found in most subduction wedges are predominantly an effect of volatile-assisted melting (Abers et al., 2014). Temperatures throughout most of the wedge are such that the addition of volatiles would lead to melting. With our Q parameterization, we test whether the observed velocities could be explained without the presence of volatiles/melts (Case A) or require a mechanism that lowers both velocity and attenuation, which most plausibly would be the presence of either of these types of fluids. We assume that the amount of fluid in the mantle wedge reduces the velocity through attenuation effects only rather than directly affecting the shear and bulk moduli (Goes et al., 2012). Without additional seismic parameters (Vp or attenuation), we are not able to distinguish between volatiles or melts. We consider an end-member distribution where fluids are distributed uniformly throughout the wedge. In addition to the cases with a dry or fluid-rich wedge, we consider one case where we assume melts are present in the upper plate mantle lithosphere and lower crust in a vertical column below the arc. We impose 10 km wide Gaussian melt distribution with a maximum of 2% melt porosity, and use the melt derivatives from Hammond and Humphreys (2000) to estimate corresponding velocity decreases. Table 1 below summarizes the four cases we examine here and their pertinent parameterizations: To compare the high-resolution seismic predictions to our shear velocity models, we apply the procedures we use in the shear velocity inversion to the geodynamic prediction. We first forward calculate the phase velocity dispersion at the periods we use in this study at each lateral point in the geodynamic predicted velocity structure. We then invert each of these dispersion curves using the same parameters settings as described in Section 2.2 to generate the 2-D transects for comparison to the preferred shear velocity model here. We assume a 20 km thick crust in the inversion.

Rayleigh Wave Tomography Results
In the phase velocity maps there are distinct anomalies associated with tectonic and geologic features (Figure 5, with 2x standard error maps shown in Figure 6). The phase velocity maps are masked using the HARMON ET AL. 0.09 km/s 2x standard error contour, which is the region we consider to be well resolved. At 22 s period, the lowest velocities (<3.5 km/s) are visible east of the arc islands near the deepest water in the regions near the trench, with relatively high velocities in the back arc (>3.8 km/s). At 22 s period and shorter, the Rayleigh waves are sensitive to water depth, with deeper water yielding lower velocities, which is likely a contributing factor to the low velocities, in the east. The higher velocities beneath the arc, likely reflect the combination of variations in water depth, sediment thickness and crustal thickness and will be discussed in more detail using the shear velocity model. At 33 s period, the highest velocities (∼4.00 ± 0.04 km/s) are found in the fore arc/outboard of the trench and in the back arc. The lowest velocities are visible in the back arc in the north, (3.75 ± 0.03 km/s), west of Guadeloupe. This pattern likely represents the change from oceanic mantle lithosphere to the east of the arc and the transition to a thickened crust beneath the arc and back arc regions. In this period range, the phase velocity anomaly patterns from teleseismic earthquakes are similar to those from noise cross correlation, and the magnitudes of the anomalies are within error of each other ( Figure S6). From 45 to 100 s period, a high velocity region (ranging from 4.10-4.20 km/s ± 0.03-0.05 km/s) is visible beneath the arc and is associated with the downgoing plate. In the same period range the backarc region has the lowest velocities (3.80-4.10, ± 0.03-0.05 km/s) with minima located west of Guadeloupe, Dominica and Martinique, and another minimum behind the Grenadines, although the Grenadines anomaly only extends from 45 to 81 s).
We also present regional averages for phase and shear velocity for completeness, bearing in mind the expected lateral heterogeneity in the region. The average phase velocity is shown in Figure 7a, along with the best-fit 1D shear velocity in Figure 7b.  in the upper crust (2.6-7.6 km depth) and 2.97 ± 0.05 km/s in the lower crust (7.6-16.1 km depth). The shear velocity model has a high velocity lid down to 40 km depth beneath the sea surface, with a maximum shear velocity of 4.58 ± 0.06 km/s and a low velocity zone beneath between 40 and 100 km depth with a minimum velocity of 4.33 ± 0.05 km/s. The model also has a high velocity region centered at 100 km depth. Overall, this is an unusual structure but could reflect averaging across a heterogeneous region that includes the presence of a slab and mantle wedge. The depth sensitivity kernels (Figure 7c) suggest for the period range used here, our resolution likely extends to depths greater than 150 km depth given the sensitivity of 100 s period extends deeper. However, our best resolution is in the upper 150 km.
We present the average parameters of our crustal structure from our 3-D shear velocity inversion, specifically the thickness of the combined sediment and crystalline crust layers and the weighted average velocity of these layers (Figure 8). There is thickened crust beneath the arc up to 32 km, and thinner crust in the fore arc and back arc. Average crustal velocities are faster beneath the arc varying from 3.60 ± 0.05 to 4.00 ± 0.05 km/s, while in the back and fore arc velocities are slower, down to 2.30-2.60 ± 0.05 km/s likely reflecting thick sediments in this region. The greatest crustal thickness, up to 32 km thick, is found beneath the platform beneath Antigua and St Kitts while the thinnest crust is found in the center of the forearc in the north and the southern back arc at 15°N ranging from 15 to 18 km. Crustal thickness is greater near the Aves ridge, >25 km, with velocities on the lower end, ∼3.3 km/s, compared to those observed within the arc.
Transects and depth slices through the 3D shear velocity model for the mantle are shown in Figures 9-11 for the two different starting models, one with slab depth from Slab2 (Hayes et al., 2018) and one with slab depth from Bie et al. (2020). The comparison of the two models allows us to assess the effect of the choice of the initial model on the final model and the range of acceptable possible structures. The fit for both models across the region is nearly identical, with normalized chi-squared residuals <1 typically, fitting the phase dispersion within the error. We compare and contrast the main features between the two models below.
HARMON ET AL.
The transects A-A′, B-B′, C-C′ and D-D' (Figures 9 and 11) show the classic structures in subduction zones, fast velocities in the downgoing slab and overriding plate and slow velocities in the mantle wedge. The overriding plate is also visible as a high velocity region >4.4 km/s in the upper 60-70 km of the model above, west of the slab. The high velocity overriding plate generally gets thinner and/or slower near the arc region, with the 4.4 km/s contour shallowing to 30-40 km within 50-100 km of the arc. In the southernmost profile D-D′ the overriding plate is less visible in comparison to the deeper part of the mantle wedge, but there is a thinning beneath the arc. The downgoing slab is well imaged and has a velocity typically >4.5 km/s in all transects with some variability of the maximum, on the order of +0.1 km/s. The model that used the slab from Bie et al. (2020) in the start model is shown in Figure 9, with the slab top (red dashed) generally following the 4.5 km/s contour at the top of the slab, although the high velocities of the slab extend some 10s of km above the slab model. The slab decreases in velocity by ∼0.1 km/s near 80 km depth in all of the transects shown using this model. The Slab2 model (solid black line) generally brackets the upper limit of the high velocity region which may be expected as this model is designed to encompass the thermal anomaly associated with the slabs rather than just follow the top of the seismicity (Hayes et al., 2018). Our inversion using Slab2 in the starting model ( Figure 11) produces a more continuous looking slab with the 4.5 km/s contour tending to follow the slab contour (black solid line). At 80 km depth, there is also a decrease in the slab velocity of ∼0.1 km/s visible in all of the transects except B-B'. The model using the Slab2 produces a sharper wedge/slab boundary and a continuous slab structure more akin to geodynamic models. Therefore, it may be the best explanation for the surface wave data. However, the model using the Bie et al. (2020) (Figure 11a). These differences highlight some of the trade-offs between the velocity structure in depth with due to the changes in crustal structure and slab structure.
The depth slices allow us to examine the lateral distribution of mantle anomalies compared with the arc, back arc and their associations with projected features on the slab. At 40-50 km depth (Figures 9a, 10a and 10b) we observe low velocity regions more or less centered beneath the arc, with the slowest velocities beneath Martinique and Dominica with small minima beneath the Grenadines and Montserrat. The slowest velocities are found between the projection of the Marathon/Mercurius and the end of the Vema Fracture Zones (∼4.1 km/s). In addition, the projection of the inferred fracture zone on the proto-Caribbean plate is located directly beneath the anomaly in the Grenadines. There is no feature directly aligned with the Montserrat anomaly. At 60-75 km depth (Figures 9b, 10c and 10d) the slowest velocities are found in the back arc behind Guadeloupe, Dominica, and Martinique, in a region that extends ∼200 km west of the arc, ∼4.15-4.2 km/s. There is also a small local minimum beneath the backarc of the Grenadines. High velocities to the east of the arc are visible at these depths and are likely associated with a thick incoming plate and possibly the cold nose. At greater depth, 100-125 km, the low velocity anomalies are ∼4.3-4.4 km/s and broadly distributed across the region. A strong high velocity region (>4.6 km/s) is visible beneath the arc, which is the subducted slab (Figure 10e and 10f).

Geodynamic Modeling
The different Cases A-D predict measurable differences in velocity structure that can be observed both before and after the inversion for the shear velocity model. The Case A provides a baseline model with a "dry" mantle wedge, where prior to the inversion process ( Figure 12a) the crust of the slab is visible and the low velocity mantle wedge has a minimum shear velocity of ∼4.3 km/s. After the inversion process (Figure 12e), the low velocity associated with the slab crust is no longer visible, and the mantle wedge velocity is faster, ∼4.4 km/s. The slab geometry and cold nose are also preserved after the inversion, as are the thicknesses of the incoming and overriding plates (Figure 12e). For Cases B and C (Figure 12b and 12c vs. Figure 12f and 12g), with a hydrated mantle wedge, the velocities in the wedge are lower, with a minimum wedge velocity of 4.2 km/s in the original model and a velocity of ∼4.3 km/s in the inversion "filtered" model. For Case C, with a sub-arc low velocity zone, the sub arc feature is preserved through the inversion process, although again with a more muted velocity anomaly after inversion. Finally, Case D (Figure 12d vs. Figure 12h) shows the effect of a thin overriding lithosphere and a hydrated mantle wedge, which results in a further velocity reduction in the mantle wedge, with velocities as low as ∼4.1 km/s prior to the inversion process and 4.25 km/s afterward.
The comparison between the seismic velocity predictions from geodynamic modeling of Cases A-D, filtered through the inversion process of our tomography shows good general agreement in the major features visible in our observed tomographic profiles. Specifically, the downgoing slab velocity in the geodynamic seismic velocity prediction of 4.65 km/s ( Figure 12) agrees well with our observation of 4.60 ± 0.05 km/s. In addition, the velocities of the overriding plate (4.6 km/s predicted vs. 4.45 ± 0.05 km/s observed) are also in good general agreement. In the mantle wedge Case A with a dry wedge is faster than the observation velocities we observe (down to ∼4.0 km/s) as can a thinner plate in Case D (down to ∼4.22 km/s). In Case D however, the lowest velocities are distributed across the entire wedge and are not focused beneath the arc, unlike our observation. The sub-arc melt Case C has a low velocity over a narrower region in comparison to our observed anomaly beneath the arc. This suggests that we require melt directly beneath the arc as in Case C, but over a broader region than imposed in the geodynamic model predictions. Since the anomaly in Case C was artificially imposed, it suggests that additional physics such as melt migration and ponding or the effects of melt on elastic moduli need to be incorporated in future modeling efforts of seismic anomalies.
The recovery test also suggests that the resolution of fine scale features such as the low velocity region on top of the slab due to the oceanic crust will not be resolved with the present methodology. In other words, we have sensitivity to relatively large-scale structure, but resolving sharp and fine features requires higher frequency methods or substantial structural information to seed the inversion. The tests also suggest that we are not likely to map thin slow features at depth like the slab top into broader shallower anomalies in the mantle wedge as they are thin compared to the depth sensitivity of the Rayleigh waves (e.g., Figure 7c). Therefore, we will limit our interpretation to the larger scale structures and trends in the model.

Crust
The crustal structure from our inversion is in good general agreement with previous work in the area. For example, our range of crustal thicknesses (27-32 km) falls within the range reported from P-to-S receiver function imaging using H-k stacking of the crust along the arc (24-37 km) (Schlaphorst et al., 2018). The locations of the thickest crust are similarly in agreement with those from receiver functions, occurring beneath the Saba/St Kitts region (30 km here vs. 37 km from H-k) and Guadeloupe (31 km here vs. 37 km from H-k) (Schlaphorst et al., 2018). More muted variability in crustal thicknesses from surface waves likely reflects their lateral averaging. Our crustal thickness pattern of thicker crust beneath the Aves Ridge (>25) and thinner crust in the Grenada basin and Tobago Trough (<20 km) also agrees well with a seismic refraction experiments across the region which found ∼26 km thick crust beneath the Aves Ridge and 17-20 km thick crust including sediments beneath the Grenada basin and Tobago Basins (Allen et al., 2019;Christeson et al., 2008). Thicker crust in the south is likely at least in part owing to a larger sediment package there, in agreement with the global sediment model and also our constraints from P-to-S conversions (Chichester, Rychert, Harmon, Allen, et al., 2020;Whittaker et al., 2013). In addition, our observed low average crustal velocities in the south may be explained by thicker sediment, as also supported by P-to-S converted phases that found >2 km of sediment south of Martinique in both the Grenada basin and the Tobago Trough, and generally thinner sediments to the north (Chichester, Rychert, Harmon, Allen, et al., 2020). A companion study, using the shorter period Rayleigh wave data from ambient noise also has a good general agreement with our findings, but provides the greater detail of the regional crustal structure (Schlaphorst et al., 2021). The general agreement with a range of independent observations gives us confidence in our crustal model, which in turn increases our confidence that our mantle structures are not biased by the shallow structure.

Overriding and Incoming Plate Structure
The convergence of the Caribbean plate and the North and South American plates provides an opportunity to investigate these different oceanic plates and make comparisons. The depth extent of the overriding Caribbean plate we observe here (60-70 km below sealevel) based on the depth to the 4.4 km/s contour is generally consistent with the predicted depth extent from thermal models of the Caribbean lithosphere for a > 50 Myr old plate near the arc and a > 120 Myr old plate behind the Aves Ridge ( Figure 1). Specifically, the depth extent is less than the 90-100 km that would be expected from the plate or Chablis cooling models (Hasterok, 2013), but in the good agreement with the 60-80 km from global compilations of lithospheric thickness from scattered wave imaging . The depth extent of the overriding plate is similar to the thickness of the >85 Myr incoming plate (∼70-80 km based on the depth to the 4.4 km/s contour) and is similar to that observed using Rayleigh waves (∼55-85 km) near the equatorial Mid-Atlantic for 40-80 Myr seafloor Rychert et al., 2021). The velocities observed at the equatorial Mid-Atlantic lithosphere and the incoming plate are higher than the overriding plate by up to 0.2 km/s.
The general difference between the slab and overriding plate may be caused by the previous period of arc volcanism at the Aves ridge, which could have thermally altered the lithosphere by the injection of hotter magma. However, it seems unlikely, given the ∼75 Myr age of the Aves ridge and the more recent extension at 50-60 Myr in the fore-arc and back arc basins of the Aves, that a thermal anomaly persists (Allen et al., 2019;Neill et al., 2011). Magma injection could also simultaneously alter the chemistry of the lithosphere via metasomatism from fluids and melt-rock reaction. Therefore, we could be seeing a predominantly metasomatic signature, which has also been proposed for volcanism associated with ocean islands (Park & Rye, 2019). The velocities at ∼80 km depth in the slab are similar to the velocities observed in the overriding plate; however, here the velocity decrease may be related to the presence of fluids, which are expected to be released from 60-100 km depth from dewatering of the oceanic crust at the top of the slab (Figure 13).

Mantle Wedge Structure Beneath the Arc
In the mantle wedge beneath the arc, the punctuated regions of low velocities (<4.3 km/s at 60-80 km depth, Figures 9-11) are likely regions of greater fluids and/or melt. As we have shown earlier, the observed shear velocity anomalies are slower than the geodynamic model shear velocity predictions with fairly large volumes of fluid present (Case B) but is similar to the models where melt is imposed in the upper plate beneath the arc (Case C, Figure 12). The fact that the seismic anomaly is larger than that predicted by the geodynamic model suggests the amount of melt is even greater and/or more widespread than that imposed in the model. The amount of fluid imposed to match the background seismic velocity throughout the wedge would also necessarily result in melt for wedge temperatures >1000°C (Katz et al., 2003). Therefore, melt likely explains the wedge anomalies beneath Montserrat, Dominica, Martinique, and the Grenadines.
Melt in these regions is in general agreement with other seismic, geologic, and geochemical observations. The low velocity region beneath Dominica and Martinique has been observed before in local body wave tomography (Paulatto et al., 2017) and also interpreted in terms of greater amounts of hydration and melt beneath the arc. The locations of our anomalies are also in good agreement with regions that have elevated rates of low magnitude seismicity (b-values > 1.0) at intermediate depths beneath the arc (Schlaphorst et al., 2016). Elevated b-values are typically associated with smaller earthquakes caused by dehydration processes in the slab (Schlaphorst et al., 2016). The arc beneath Martinique and Dominica is also the region of the greatest volume of arc volcanism over the past 100 kyr (Wadge, 1984), while Montserrat, St. Vincent and Kick-em Jenny in Grenada are actively erupting today. Furthermore, there are geochemical proxies for fluids such as B/Nb that indicate relatively high fluid input into the central arc (Cooper et al., 2020), with another notable anomaly near Grenada. The fact that the sub-arc anomaly in the center of the arc is stronger than that beneath the Grenadines and Montserrat could reflect different amounts of melt and fluids fluxing into the system.

Relationship of Wedge Low Velocities to Projected/Inferred Slab Features
The punctuated mantle anomalies could be caused by variability in fluid released into the arc from the slab, and there are several possible causes for variability. First, highly tectonized fracture zones may deliver greater hydration and may have more alteration to greater depth. The location of the sub-arc anomaly beneath the Grenadines aligns with the predicted current location of an ancient fracture zone on the proto-Caribbean plate (Cooper et al., 2020). The anomaly beneath Martinique and Dominica is centered between two fracture zones and potentially the proto-Caribbean spreading ridge intersection, but not directly above either. There is no fracture zone or plate boundary associated with the anomaly beneath Montserrat.
A second possibility is fluid may be supplied from the boundary of the proto-Caribbean lithosphere with the Mid-Atlantic lithosphere and/or the extinct and subducted proto-Caribbean spreading ridge. The region may be highly tectonized as it is the result of the interaction of a spreading system and a long offset fracture zone. In addition, this part of the plate may be weaker in general and may be reactivated at the outer rise, allowing further hydration at the outer rise. Specifically, we note that the seismicity on the incoming plate today has a strong strike slip component in the region south of the Marathon/Mercurius Fracture Zones and has a strike direction that is in good agreement with the inferred proto-Caribbean boundary strike and the strike of the fracture zones of the proto-Caribbean ridge ( Figure 1). Our largest anomaly is centered just south of this boundary, and so again not directly aligned, but possibly related.
A third possibility is that hydration of the plate occurs at the ridge and varies according to ridge processes that are thought to vary temporally at slow-spreading centers like the MAR creating abyssal hill and detachment fault morphology (Escartin et al., 2008;Smith et al., 2006). The variation in morphology does appear to be linked to variable alteration of the upper 20 km of the incoming plate from the VoiLA active source study (Davy et al., 2020). The deepest alteration is found in the ridge segments adjacent to the Marathon Fracture zone rather than directly in the fracture zone itself (Davy et al., 2020). There is no way to know the ridge processes and resulting morphology of the subducted plate, and so direct connection with our result is speculative.
A fourth possibility is that aside from discrete tectonic features on the downgoing plate, variations could also be caused by differences in the origins, histories, and therefore likely compositional and/or tectonic differences in the downgoing plate along the arc. In the north the plate originates from the mid-Atlantic ridge, whereas less is known about the subducted plate in the south, which is thought to be proto-Caribbean and material from the edge of the continental margin. South of the Grenadines the lack of seismicity and reduced fluid flux may be related to there being a significant amount of South American continental margin material on the proto-Caribbean plate from the initial separation of North and South America. This material may be less altered as it may not be as tectonized or more discretely tectonized along fault systems. Similar observations have been recorded on the Galicia margin (Bayrakci et al., 2016). In contrast, in the north, the downgoing plate may be more hydrated, given that it was formed at a slow spread oceanic HARMON ET AL.
10.1029/2021GC009707 21 of 27  (2008), Phase stability was calculated using the data base from Holland and Powell (1998) using PerPleX (Connolly, 2005) and the solid-solution model selection from Hacker (2008). Slab geotherms at the slab top from thermal model shown in Figure 12j, the slab Moho (6 km depth) and 20 km deep into the slab show that the steeper slab below Dominica (Do) can retain water to larger depths than the more shallowly dipping slab below Grenada (Gr).
lithosphere. However, our strongest wedge anomaly occurs above the proto-Caribbean slab. Therefore, this explanation by itself is not well supported.
A fifth possibility is that the different origin of the northern versus the southern slab may affect hydration entering the plate at the outer rise, since the expected strike direction of the southerly proto-Caribbean ridge generated fault systems may be more oblique to the subduction direction and may not be reactivated at the outer rise as easily in the south. However, we note that the outboard seismicity is sparse here, and so this may not be an important factor (e.g., Figure 1).
The sixth and final possibility is that variation in sediment input along arc could influence the amount of water in the wedge. However, the fact that most sediment is in the far south and our strongest anomaly is in the central part of the wedge again rules out a direct relationship. Finally, sediment could influence fluid input by acting as a caprock, limiting the amount of water that can enter the plate at normal faults at the outer rise in the south. This could explain why our strongest anomaly is in the central part of the arc, but not the fact that there is a lack of a strong anomaly in the north.

Back Arc Anomalies
The low velocities extending >200 km into the back arc west of Dominica/Martinique, likely explained by melt beneath the lithosphere-asthenosphere boundary, are not predicted from geodynamic models that assume vertical ascent of fluids and melt ( Figure 10d). Thermal models for the Lesser Antilles suggest that there should be two primary pulses of volatiles from the slab, one at ∼80 km depth from the dehydration of oceanic crust directly below the model coupling transition depth and one at ∼150 km depth from the dehydration of peridotite directly below the Moho (solid and long-dashed lines, Figure 13) (van Keken et al., 2011). Volatiles that are released and travel vertically from these depths roughly bracket the fore-arc and extend ∼50 km into the back arc, in other words ∼150 km shy of the extent of our anomaly, depending on the slab top model. The mantle wedge is thought to be relatively cold due to the low convergence rates of the system which result in strong conductive cooling from the slab and overriding plates, likely with mantle wedge temperatures <1350°C (Syracuse et al., 2010, this study). In order for melt to be generated in the backarc in the ambient mantle, the wedge would need to be more hydrated than typical ambient mantle ∼100 ppm water (Bell & Rossman, 1992) and/or melt and fluid pathways are not vertical.
Hydration of the back arc could be achieved in a variety of ways. First, the mild association of the anomalies with fracture zones and the proto-Caribbean lithosphere boundary, suggests that these regions might be regions of deep hydration (>20 km into the slab typically assumed) which would allow them to carry water to greater depths before crossing the 600°C isotherm and dehydrating (Figure 13a and 13b short-dashed lines) (Schmidt & Poli, 1998). This could extend the release of the volatiles 150 km into the backarc wedge at depth to slab of 180 km, which is in excellent agreement with the western extent of our back arc anomaly, depending on the slab top model. There is some evidence for this from the VoiLA active source study out board of the trench of the Lesser Antilles, which found lower 5%-10% P-wave reduction to ∼20 km depth in the ridge segment adjacent to Marathon fracture zone indicating serpentinization deeper than the Moho (Davy et al., 2020).
Another geometric effect of changing slab dip might also allow for variation in the depth of the release of volatiles. Specifically, slab geotherms for this slow subduction zone almost parallel dehydration boundaries (Figure 13), so water may be carried up to 50 km deeper along the steepest slab in the central arc than in the significantly more shallowly dipping slab at the southern and northern end of the arc.
Along with the geometric effect, the 3-D corner flow may be particularly important here given the strong curvature of the arc along with diapiric flows. In addition, given that the slab has rolled back in its history, excess fluids or melt in the backarc could have been delivered in the past, at a time when the slab extended further into the backarc (Allen et al., 2019;Neill et al., 2011), although it is not clear if fluids and melt would linger over long geologic timescales. Similarly, diapiric flows of fluids, sediments and melts of the slab might be released from the deeper parts of the slab as cold plumes of material from 180 to 200 km depth to slab (Gerya & Yuen, 2003). These different types of solid flows might also aid or enhance fluid/melt migration.
Finally, fluid/melt porous flow dynamics of the system might effectively mix mildly hydrated mantle peridotite into the back arc. Several geodynamic models of porous flow in subduction suggest complex behaviors as the solid flow field interacts with the fluid/melt phase flow, rather than simple vertical buoyancy driven flow of the fluid phases. For example, if permeability is relatively low due to smaller grain sizes, the downgoing slab flow can entrain the fluids and drag them deeper into the system, which could be brought up further into the backarc (Cagnioncle et al., 2007;Wada & Behn, 2015). Small grain sizes are encouraged by low temperatures, as expected in the Antilles subduction zone, and are likely given the fast velocities observed above the Benioff zone. Compaction effects can also allow fluids released from the slab to migrate upwards along the slab interface but produce a widely hydrated mantle wedge when permeability is moderate to low . Fluids also can ascend through the mantle wedge and form high porosity channels at the base of the overriding lithosphere , similar to what we infer from our observed structures in the back arc and also from S-to-P Imaging . Finally, there could be some degree of ponding of melt beneath the upper plate. This could be relatively unique to the Lesser Antilles. Typically, the upper plate is not well-imaged in subduction environments, either because it is very thin or hot owing to thermal erosion, although it is more distinct in our model of the Lesser Antilles.
Overall, our velocity model supports the influence on fluids of fracture zones and/or plate boundaries but with additional influence from 3-dimensional flow and/or the effects of diapiric flow, melt compaction or ponding. There is no apparent direct relationship between the amount of sediment on the plate or the different tectonic origins of the northern versus the southern slabs. The lack of normal faulting at the outer rise does not necessarily support the importance of fluid input at the outer rise. Therefore, variability in anomaly structure along arc could be influenced in the fabric of the tectonic plate as it is created near the ridge, although these conditions are unknown for the subducted slab. Alternatively, the correlation between the location of the fracture zone in the south with the anomaly beneath Grenadines and the proximity of the plate boundary and fracture zones to the anomaly beneath Dominica and Martinique suggests that fracture zones play an important role. This notion is also supported by geochemical evidence for dewatering of serpentinized oceanic lithosphere in the vicinity of the seismic anomalies (Cooper et al., 2020). However, the fact that there is not a 1-to-1 relationship and also that there is a low velocity anomaly in the back arc that is not easily explained suggests that factors such as 3-D flow, particularly in the vicinity of the strongest arc curvature, and/or entrained sedimentary plumes, melt compaction, and/or ponding also play a role.

Global Comparison
Our results add to the global catalog of subduction zone images, and the comparisons between our velocity structures particularly in the mantle wedge yield some insights into the effects of slow subduction of old slow spread lithosphere on subduction processes. Specifically, while we observe low velocities in the mantle wedge, down to 4.1 km/s, the minimum velocity is relatively high compared to faster converging systems where the same method has been used such as Central America (Harmon et al., 2013), the Aleutians (Wang & Tape, 2014), Marianas (Pyle et al., 2010) and Tonga Lau (Wei et al., 2015), which have minimum velocities of 4.05, 4.00, 4.00, and 3.60 km/s respectively. These four studies used the same methodology, so the absolute velocities should be generally comparable. All of these studies invoked melt in the mantle wedge to explain their low velocity anomalies, which is likely given the tectonic environment, and it also suggests that there may be different amounts of melt present in these systems. In addition, the along-arc length scales of our sub-arc anomalies are generally more focused here than in other mantle wedges where they may extend >200 km along the arc. In most of these studies low velocities were located directly beneath the arc and did not extend more than 50 km into the back arc in comparison to the 200 km extent in the Lesser Antilles. The exception is Tonga Lau, which had strong anomalies in the back arc associated with a back arc spreading center, in contrast to the lack of an active spreading center in the Lesser Antilles. In the fast-converging systems, temperatures may be higher, melt production may be greater, and permeability higher allowing vertical ascent of melt to dominate. Whereas in the Lesser Antilles lower melt production may lead to lower permeabilities, allowing for greater hydration into the back arc by the solid flow. Therefore, even if melt is ephemeral (Sim et al., 2020) and seismic tomography is imaging melt at different stages Rychert et al., 2020;Wang et al., 2020), a general relationship of the increasing convergence rate with increasing temperature, melt, and volcanism holds.

Conclusions
We have presented shear wave velocity model of the Lesser Antilles arc, imaging both the incoming plate, subducted slab, mantle wedge, and upper plate. We observe the lowest velocity region (4.1 km/s) in the mantle wedge beneath the arc near Dominca and Martinique, that is likely related to the release of volatiles from fracture zones on the downgoing slab and the boundary between proto-Caribbean and Mid-Atlantic generated lithosphere, and also associated melt. Lack of a direct correlation with several variations in observables related to fluid along arc suggests the influence of 3-dimensional mantle flow at least near the central region of the arc, i.e., that with strongest curvature. The anomaly also extends 200 km into the back-arc suggesting either deep release of fluids from the downgoing slab in the region, possibly due to deep hydration of the downgoing slab at the fracture zones, ponding/compaction of ascending melt, or sediment plumes. Both active source imaging of deep hydration of the incoming plate, and geodynamic modeling that predicts full dewatering of the slab near the western extent of our anomaly lend support to deep release of fluid. There is also a more focused and muted low velocity anomaly in the southern arc beneath the Grenadines, which is directly above an inferred fracture zone on the proto-Caribbean plate, which again suggests greater fluid flux from a deeply hydrated fracture zone on the downgoing plate. Compared to other subduction zones, our lowest velocities in the mantle wedge are more focused beneath the arc in the alongarc direction and faster by 1%-12% than the lowest velocities observed in other subduction zones around the Pacific, suggesting less melt in the Lesser Antilles. This notion is generally consistent with low arc volumes in the Lesser Antilles and the inference of a cooler mantle wedge due to the slow convergence.

Data Availability Statement
Data used in this project is available from the IRIS DMC (www.iris.edu/dmc) and we used of data from the