Mapping geologic features onto subducted slabs

Estimating the location of geologic and tectonic features on a subducting plate is important for interpreting their spatial relationships with other observables including seismicity, seismic velocity and attenuation anomalies, and the location of ore deposits and arc volcanism in the over-riding plate. Here we present two methods for estimating the location of predictable features such as seamounts, ridges and fracture zones on the slab. One uses kinematic reconstructions of plate motions, and the other uses multidimensional scaling to ﬂatten the slab onto the surface of the Earth. We demonstrate the methods using synthetic examples and also using the test case of fracture zones entering the Lesser Antilles subduction zone. The two methods produce results that are in good agreement with each other in both the synthetic and real examples. In the Lesser Antilles, the subducted fracture zones trend northwards of the surface projections. The two methods begin to diverge in regions where the multidimensional scaling method has its greatest likely error. Wider application of these methods may help to establish spatial correlations globally.


I N T RO D U C T I O N
The geologic features and tectonic structures on the downgoing plate of subduction zones potentially have an impact on seismogenesis, fluid release, ore formation and volcanism. For example, it has been suggested that fracture zones on the downgoing slab are linked to higher b-values, that is, greater amounts of small seismicity (Lange et al. 2010;Schlaphorst et al. 2016). Seamounts have been proposed to be asperities for large ruptures and/or weak zones during rupture, where their location has been identified with some certainty with seismic refraction (Kodaira et al. 2000). Seamounts have also been hypothesized to cause creep due to the evolution of the surrounding fault system (Wang & Bilek 2011) potentially explaining seismic gaps in some subduction systems such as the inferred location of the subduction of the Louisville Seamounts in Tonga (Bassett & Watts 2015). Furthermore, fracture zones have been proposed to be barriers for rupture, for example during the 2004 December 26 Sumatra megathrust earthquake (Robinson et al. 2006). The subduction of buoyant features such as oceanic plateau, ridges and seamounts has been linked to the formation of magmas capable of forming porphyry-type ore deposits (Cooke et al. 2005;Rohrlach & Loucks 2005;Rosenbaum et al. 2005) and the subduction of fracture zones has also been proposed as an important trigger for porphyry ore formation (Richards & Holm 2013). These connections are fundamental to our understanding of the driving forces and determining factors of observables like style of seismicity, conditions of melt production and migration to volcanoes. However, as the slab descends to greater depths, the simple projection of seamounts, fracture zones, and ridges becomes more tenuous, as do the links between features on the downgoing plate and subduction zone processes.
Here we examine the problem of estimating the location of geologic features on the downgoing plate in subduction zones. We focus on features such as fracture zones and linear seamount chains, given that their general linearity lends itself to a greater degree of predictability for projection onto the downgoing slab. We present two methods: the first is a kinematic approach and the second is a method based on multidimensional scaling (MDS). Both methods uphold conservation of surface area and mass of the plate assuming there is no internal deformation of the plate such as stretching or compaction. We present synthetic examples of the methods and an example from the Lesser Antilles subduction zone. The methodology can be generalized to any subduction setting.
We use the Lesser Antilles ( Fig. 1) as our real test case. It is an ideal location because it has several examples of predictable geologic features that have been subducted, a curved slab geometry, and a well constrained tectonic history for the plates interacting at the boundary. Specifically, there are several fracture zones entering the subduction zone, including the Fifteen-Twenty, Marathon, Mercurius, Vema and Doldrums fracture zones. These are longlived fracture zones and can be traced to either side of the Atlantic basin. The opening of the Atlantic is well characterized by magnetic anomalies and small circles of fracture zones (Klitgord & Schouten 1986;Cande et al. 1988;Muller & Roest 1992), meaning that the trends or flowlines of these fracture zones are well constrained. The strong curvature of the trench and subducted slab creates a scenario where tracing these features into the Earth is not necessarily intuitive, making it an ideal location to test our methods. Finally, despite some debate over the precise tectonic evolution of the Caribbean plate (Boschman et al. 2014), its motions relative to the North and South American plates are known (Matthews et al. 2016).

Fracture zone flowline estimation
In order to project geologic features onto a downgoing slab, an estimate of where those features existed on unsubducted plate must be determined first. In some cases, such as seamount chains, the trends of the unsubducted seamounts can be used to estimate the likely locations of the subducted seamounts. In the case of fracture zones, small circles or plate flowlines of their trends should be used, incorporating information from the spreading history and the conjugate margins. Here we detail how we estimate the likely trends of fracture zone locations prior to subduction for our Lesser Antilles example.
Beneath the Caribbean, material from both the North and South American plates is being subducted (Klitgord & Schouten 1986;Muller & Roest 1992;Muller & Smith 1993;Müller et al. 1999), but the present day boundary is not well constrained (DeMets et al. 2010). Muller & Roest (1992) previously showed that the Triple Junction between the North and South American and African plates moved northwards to its present position close to the Fifteen-Twenty Fracture Zone during the Cenozoic. Therefore, in order to accurately project the surface locations or flowlines of the fracture zones on the North and South American plates we must consider the time history of the plate boundary evolution and allow for changing poles of rotation and plate pairs.
We computed synthetic flowlines within the GPLATES software using combinations of northwest African and North and South American rotation poles as given in Seton et al. (2012). This was accomplished by choosing seed-points near the present-day fracture zones of interest (Fifteen-Twenty, Marathon, Mercurius, Vema, and Doldrums, Fig. 1) in GPLATES and then creating a 'Flowline' feature in time steps of 1 Myr from the present to 120 Ma, the approximate time of the opening of the southern Atlantic. The poles of rotation used can be adjusted manually within GPLATES. We assume symmetric spreading. Although other authors have noted asymmetric spreading at various time periods in the region (e.g. Bird et al. 2007;Schettino & Macchiavelli 2016), uncertainty in the asymmetric component is high. These synthetic flowlines were then compared to fracture zone traces observed from modern satellite altimetry data on both the American and African plates (Wessel et al. 2015). In addition, the synthetic flowlines were checked against magnetic anomaly identifications where they coincided. This was difficult in some areas because the locations were at low magnetic latitude or near closely spaced, large-offset transforms.
Our analysis confirms the interpretation of Muller & Roest (1992), with the older parts of both the Mercurius and Marathon fracture zones being best matched with North American-African poles, and the younger parts with South American-African poles. The African-North American-South American Triple Junction passed the Mercurius fracture zone at 60 Ma and Marathon fracture zone at 50 Ma (leading to the bends visible in Fig. 1b). Therefore, the Mercurius and Marathon fracture zones were modelled using hybrid poles accordingly. Fracture zones to the south (Vema and Doldrums) were modelled with South American-African poles, and the Fifteen-Twenty Fracture Zone with North American-African poles throughout. Given the age of Atlantic opening in this region (at Downloaded from https://academic.oup.com/gji/article-abstract/219/2/725/5531317 by Durham University user on 31 October 2019 most 120 Ma), the Doldrums Fracture Zone is not likely to underlie the present Lesser Antilles arc and Vema may only extend to the arc ( Fig. 1).

Kinematic projection of geologic features
In this approach we track the geologic feature through time as it enters the subduction zone and moves down along the slab based on the relative motion between the two plates, that is, using a kinematic reference frame. In the subduction zone, the local convergence direction is downdip of the slab. This method requires an estimation of the trench location and slab geometry through time, relative to the geologic feature of interest, and also the relative velocities between the plates through time. We assume the slab geometry and boundary have been stable through time, that is, while the slab may advance/retreat the geometry of the slab itself does not change. Although this could be debated, the assumption is not inconsistent with plate reconstructions (Matthews et al. 2016).
To illustrate the method, we present a synthetic example first ( Fig. 2a). We use a subducting slab surface described by a parabolic function (z = 0.001x 2 + 0.05x + 6 in km) defined for positive values of x and 6 km deep otherwise for the seafloor of the incoming plate. We introduce a linear geologic feature (e.g. a fracture zone) on the incoming plate outside of the subduction zone (negative x) at 60 Ma, with a trend defined by y = 0.25x. The downgoing plate has a convergence velocity (v c ) of 10 km Myr −1 along the positive x direction. Each part of the linear geologic feature then moves at each time step (1 Myr) at the convergence velocity. As the downgoing plate dips into the subduction zone the effective magnitude of the convergence velocity at each time step on the surface of the Earth is described trigonometrically by v x = |v c |cos(tan −1 ( v c • ∇z/|v c |)), where vertical bars indicate the scalar magnitude.x In our illustration, we find that the geologic feature has diverged >5 per cent from its original surface projection after 30 Myr of subduction, and after 60 Myr the surface projection of the subducted feature has been deflected significantly in the negative x direction, rotating counter-clockwise. Matlab scripts to reproduce this synthetic are available in the Supporting Information.
For the Lesser Antilles example, we estimate the relative motion of the present-day trench with respect to the incoming plate over the past 30 Myr. In our example, we use the relative motions between the South American Plate and the Caribbean Plate. We use the 6 km depth below sea level slab contour (approximately the trench) of the Lesser Antilles subduction zone as our overriding plate front (Syracuse & Abers 2006) and we use GPLATES software to track the motion paths of the points on our trench (Gurnis et al. 2012) through time using the reconstruction from Matthews et al. (2016) in 1 Myr time steps. To accomplish this in GPLATES, we first load the 'Data Bundle For Novices.proj' file which contains the Matthews et al. (2016) reconstruction. We then use our digitized trench locations in the present day as the seed points for 'Motion Path' features for the relative motion paths between the South American Plate and the Caribbean Plate from 0 to 30 Ma, assuming the South American Plate is fixed. The trench location through time is shown in Fig. 3(a). The outputs of the locations of the trench points are then used to determine the local linear velocity magnitude and azimuth at each time step.
In the second step, we determine where the fracture zones on the South American Plate enter the Lesser Antilles trench through time and then subduct the features down the slab. The fracture zone locations we use are those determined according to the procedure in Section 2.1. We map the Fifteen-Twenty, Marathon, Mercurius and Vema fracture zones to the trench locations. We use the motion paths calculated from GPLATES to determine the local convergence velocity magnitude and direction along the trench. Each point of the fracture zone is then propagated into the subduction zone with the local convergence velocity, resolved into the local surface of the slab at depth. The slab velocity into the Earth is equal the local convergence velocity at the free surface assuming that the slab is not being compressed or stretched. The streamlines of the points are tracked through time until the present day and the final location (longitude, latitude and depth) of each point on the fracture zone is determined. Matlab scripts and associated data for the Lesser Antilles example are available in the Supporting Information.

Brief review of mathematics of multidimensional scaling
The second approach we use is MDS (Kruskal 1964), which is a means of mapping an arbitrary surface defined in 3-D onto a surface confined to a plane or the Earth's ellipsoid-thus reducing the dimensionality. In our case it is taking a subducted slab defined in depth and transforming it to its unfolded position as a tectonic plate at the surface of the Earth. This is done via a transformation between the along-surface distance (referred to in the literature also as geodesic distance but not strictly confined to an ellipsoid) of the arbitrary slab surface to another distance (e.g. Euclidean, Geodesic) approximation confined to a 2-D coordinate system. We follow an approach adapted for texture mapping onto 3-D surfaces (Zigelman et al. 2002). This method effectively unfolds the slab and projects it back to its original position at the surface, nominally conserving surface area and/or mass without internal deformation of the plate. A similar approach has been applied for plate reconstructions when it is desirable to know where subducting slabs were in the past (Wu et al. 2016).
We briefly review the mathematics behind classical MDS and its relationship to other forms of MDS following Zigelman et al. (2002). Using classical MDS, we seek a 2-D Euclidian representation, X, of a surface defined in 3-D, given the squared distances along the surface between all possible pairs of points in 3-D in M. The size of X is n × 2, and that of M is n × n, where n is the number of points of interest on the surface. The squared Euclidian distance, E, between each point of X is given by where i and j are indices for the rows and columns of each point and a is the index for the two dimensions for the elements x for X. Or written compactly, where c is the vector of the magnitude squared of each element of X, and 1 is an n × 1 vector of ones and where T denotes transpose. E is n × n. We map M onto E for our transformation; however, we need to centre and normalize the distances in M so they are comparable to our presumed X, which has its centroid about zero. We define a matrix J, which we use to centre and normalize X. The matrix has the following properties, J1 and 1 T J gives a vector of zeros and for centred coordinates JX = X. J is given by where I is the identity matrix. Applying normalization by J to both M and E, multiplying by -0.5, and equating them yields we now have: Now we define a new matrix B: The classical MDS 'strain' minimization problem is given by where b ij is the i,j th element of B. We note that 'strain' is the nomenclature used to describe classical MDS, but it is not related to tectonic strains in the slab. A centred version of the new coordinates can be recovered by a singular value decomposition of B = U U T , where U is the matrix of eigenvectors and is the diagonal matrix of eigenvalues. X can then be recovered: where U is now the two eigenvectors corresponding to the two largest eigenvalues on the diagonals of . An alternative to classical MDS is given by 'stress' minimization, which is achieved through the following optimization of X: where d ij is the along-surface distance between points i and j, or (M ij ) .5 and f(x i, x j ) is a distance function for the new projected coordinates X, which can either be Euclidian distance in Cartesian coordinates or great circle distances on the Earth. We note that 'stress' minimization is the nomenclature used in the literature, but it is not related to the tectonic stress in the subducted plate. We examine both the Euclidian distance and great circle distance realization for f(x i, x j ) of the stress minimization. In the Euclidian distance stress minimization MDS, the classical MDS in Cartesian coordinates is used as the starting point for a gradient based iterative process to minimize the objective function (Kruskal 1964).
In the great circle distance stress minimization MDS, we optimize directly in longitude and latitude using our subducted slab longitude and latitude points as our starting point and constrain the points on the trench to remain stationary. We prefer the great circle distance stress minimization MDS in geographic coordinates and great circle distances and present the results below.

Multidimensional scaling projection of geologic features
We first illustrate the MDS method via a synthetic example, using the same geometry as in our kinematic synthetic example (Fig. 2, bottom) using only classical MDS and Euclidean distance stress minimization MDS, due to the Cartesian setup of the example. The grey surface at the top of the blue slab feature is the result from the MDS, but we have truncated the x limits at 500 km for presentation purposes. Both classical MDS and stress minimized MDS returned nearly identical results (Fig. 4), with little distortion of the slab in the y-direction and good agreement with the analytical expectation for the unfolding of the slab in the x-direction (i.e. the line integral of z). Errors were typically 1 per cent. Matlab scripts to reproduce this example are available in the Supporting Information. For our Lesser Antilles example we use the following steps for MDS to project geologic features. In the first step, the along-surface distances between all points of interest of the subducted slab are determined to construct M. We grid and triangulate the slab depths for our synthetic example and the Lesser Antilles example (contours from Syracuse & Abers 2006) at a 1 × 1 km horizontal resolution in a Cartesian coordinate system. The slab depths for the Lesser Antilles are transformed from geodetic coordinates to Cartesian using the WGS 84 reference ellipsoid and the median longitude and latitude as the origin. We then use a fast marching algorithm (Sethian 1999) on the 1 km resolution grid to estimate the alongsurface distances between every 20 × 20 km horizontally spaced point of the grid using an implementation in Matlab (Peyre & Cohen 2008). The 20 × 20 km spacing reduces the number of times the distance calculation needs to be carried out, but the 1 × 1 km grid ensures the distance calculations are sufficiently accurate. The squared distances are then compiled into a matrix M, where M ij is the squared along-surface distance between points i and j of the down sampled grid.
In the next step, the MDS is carried out to unfold the slab. We examined classical MDS and 'stress' minimization MDS using Euclidean distance and also great circle distances, as described in Section 2.3. In our case, X, the output from classical MDS and Euclidian distance stress minimization MDS in Cartesian coordinates is then rotated and translated so that the points that lie on the trench align with unscaled original points on the trench as best as possible. The points in the Cartesian coordinates are transformed into geographic coordinates (Fig. 3). For the MDS using great circle distances, the points are constrained at the trench, so no rotation or translation is required.
In the final step, the surface projection of the fracture zones is mapped down onto the subducted slab at depth via the interpolated proximity to the points of surface projection of the MDS unfolded slab. In other words, each point in the unfolded slab at the surface corresponds to a point on the subducted slab. This allows us to map the fracture zone flowlines at the surface to the subducted slab by their relative location to the unfolded slab points. Matlab scripts and associated data for the Lesser Antilles example are available in the Supporting Information.

R E S U LT S
The synthetic examples show good agreement with each other in their performance for projecting a linear feature onto the slab at the depth. We compare the two methods to the analytic prediction for the case of a linear subducted slab that only varies in x and z (Fig. 4). In this case, the surface projection our linear feature that is being subducted will only have its x-values modified as the slab descends into the Earth. The theoretical prediction is determined by rescaling the x-axis of the linear feature to the arc length integral of z (distance along the slab surface). In our case this has an analytic form for the distance along slab d as a function of x, d(x) = 250 * asinh(x/500 + 1/20) + 250 * (x/500 + 1/20) * ((x/500 + 1/20)∧2 + 1)∧(1/2), which effectively modifies the x-axis due to the bend in the slab, making the new linear feature projection y = 0.25 * d(x). There is no difference visible between the kinematic and theoretical predictions, while there is some small error visible with the MDS estimate. This error likely arises in the MDS due to the interpolation between the projected points. The MDS estimate is within <2 km of the theoretical prediction, suggesting for subduction systems with linear trenches both methods produce reasonably robust results. We compare the predictions for the surface projections of the unfolded slab in the Lesser Antilles from classical MDS, Euclidian distance stress minimized MDS and great circle distance stress minimized MDS (Fig. 5). Classical MDA and Euclidian distance stress minimized MDS tend to distort the shape of the trench, even after rotation and translation. This leads to differences in the unfolded slab trench location and likely errors in the final projections of the geologic features, specifically leading to slight shifts in the fracture zone flowline position at the trench when projected down onto the slab. The great circle distance stress minimized MDS preserves the location of the trench and the fracture zone flowline position at the trench when projected onto the subducted slab.
We assess the quality of the MDS stretching by determining how well the surface area of the downgoing slab is preserved after projection, both in each facet of the triangulation and the total surface area. The total slab surface area is well preserved: 94 per cent after classical MDS, 103 per cent after the Euclidian distance stress minimized MDS and 98 per cent for the great circle distance stress minimized MDS. Individual triangle areas are generally within 10  Syracuse & Abers (2006;SA06). The black polygon shows the unfolded slab from the great circle distance stress minimized MDS (GCMDS), the red polygon shows the prediction from the Euclidian distance stress minimized MDS (EMDS), the green polygon shows the prediction from the classical MDS (cMDS) and the cyan polygon shows the slab surface as determined by earthquakes (Slab SA06).
per cent of their original surface area for all methods. We show the error in the individual triangles for the great circle distance stress minimized MDS in Fig. 3(b). The regions with the highest fractional error are the deeper parts of the slab at the northern and southern extent of the subduction zone and the positive sign reflects the underprediction of the slab surface area. In our great circle distance stress minimization MDS the errors in the individual triangles are slightly higher by 0.05-0.10 fractional error relative to the classical and Euclidian distance stress minimization MDS, particularly at the edges of the unfolded slab. This difference likely due to the imposed constraint of the immobile trench points. We present results from great circle distance stress minimized MDS and discuss those results exclusively below. In our Lesser Antilles example, the projections of the fracture zones on the subducted slab do not follow the flowline projections of the fracture zones at the surface due to the changes in the relative convergence direction of the Caribbean and South American plates (Fig. 6). For the Vema, Marathon and Mercurius fracture zones, both methods predict trends (yellow circles/red crosses) that rotate 5-10 • clockwise from the fracture zone flowlines at the surface (black dashed lines). This is visible in predictions from both the kinematic and MDS methods. The Fifteen-Twenty Fracture Zone predictions are close to the trend of the fracture zone flowline at the surface for both methods. The MDS method predicts a counterclockwise deviation from the fracture zone flowline at the surface west of −63 • Longitude, whereas the kinematic estimate is nearly aligned with the flowline.
The predictions between the two methods are, in general, consistent with each other. The MDS estimates tend to be smoother than the kinematic estimates which display distinct kinks related to the changes in convergence direction and velocity, particularly in the last 4 Myr. The most significant differences in the predictions are for the Fifteen-Twenty Fracture Zone which is more oblique to the dip of the downgoing plate; this is also in the region where the fractional error in the MDS estimate is highest.

D I S C U S S I O N
We find generally good agreement between the projections from the kinematic and the MDS methods for both the synthetic and Lesser Antilles examples. Given that there are inherently different assumptions in the approaches, this corroborates the solutions. Overall, we are confident in the locations where the kinematic and MDS projections agree. In the locations where they differ, such as the westernmost part of the flowline for the Fifteen-Twenty Fracture Zone (blue circle versus red cross, Fig. 6), the difference constrains a range of possible projections. Lack of agreement with the flowline continuation (black dashed versus red crosses for MDS and blue circles for the kinematic method in Fig. 6) highlights the necessity of applying at least one of the methods.
Uncertainty in the feature mapping arises from the different approaches for conservation of the surface area/mass of the slab. The difference could either be from errors in the surface projection performed using the MDS approach or a poorly constrained plate tectonic convergence and changes in plate boundary and slab structure through time in the kinematic approach.
In the kinematic method, the geologic features are located by simply moving material forward down the slab at the convergence velocity. Our approach assumes a consistent shape of the plate boundary and slab through time. Although it would be possible to account for a variable shape in the boundary through time, the variability would also need to be well constrained. In addition, the kinematic approach is only as accurate as the plate tectonic reconstruction.
Uncertainty in the MDS method arises from imperfectly unfolding the slab and mismatch between the unfolded slab edge and the trench, which is ameliorated in our preferred great circle distance stress minimization MDS. In our example, 98 per cent of the surface area is retained, which suggests there may be some inaccuracy. Some of this inaccuracy may be physically based, especially if there is a discontinuity in the slab structure, such as a tear, that is smoothly mapped into the slab structure. In addition, the projection will only be as good as the slab model, which is dependent on being able to image the slab through seismicity. In the Lesser Antilles, a slab tear has been hypothesized near Dominica (van Benthem et al. 2013), where the North American and South American Plate Boundary may exist. In addition, the Barracuda Ridge and Tiburon Rise, outboard of the trench, are thought to be compressional features indicating convergence between the North and South American plates together with possible flexing as they approach the trench (Pichot et al. 2012). A sharp jump in the fractional error of the surface area in the near St. Vincent to Barbados may be in part due to a slight discontinuity of the slab due to a tear or compression. The distortion of the trench by the classical MDS and Euclidian stress minimization MDS (Fig. 5) might also be explained by internal deformation of the slab and/or discontinuities. Further analysis may be possible to determine if the reconstructed slab at the surface represents a continuous structure at depth or whether there has been deformation of the downgoing slab by analysis of the error of the individual facets of the triangulation as shown in Fig. 3(b). This is the subject of future work.
The Lesser Antilles is noted for the striking contrast in levels of seismicity between the northern and southern sectors of the arc. High b-values from Gutenberg-Richter earthquake magnitudefrequency relationships beneath the Lesser Antilles Arc have been attributed to deep release of water from fracture zones (Schlaphorst et al. 2016). We find good agreement between our projection of the fracture zones and the region of high b-values (>2) beneath Martinique for the Vema fracture zone. The broad region beneath Dominica and Martinique with b-values > 1.6 appears to have its northern limit demarcated with the Marathon and Mercurius fracture zones, suggesting the region in between the three large offset fracture zones is related to the high rates of seismicity. Therefore, the differences we observe in the fracture zone locations at depth relative to their surface flowlines would likely affect the interpretation of the high b-values in other locations.

C O N C L U S I O N S
We have presented two methods for mapping predictable geologic features on downgoing plates in subduction zones, one using a time varying kinematic approach and the other using MDS. We find good overall agreement between the predictions of the two methods, and a significant difference between these and simply using flowlines to project fracture zones in the Lesser Antilles. Small differences between the kinematic and MDS methods can be attributed to the fact that the two methods required either assumptions about the plate reconstruction and slab geometry through time or the geometry and continuity of the slab, respectively. Applying the methods together allowed us to cross check agreement and give greater confidence in the range of acceptable predictions. Future analysis using these methods will help to clarify the relationships between geologic features on the downgoing plate and arc volcanism and ore formation, seismicity, and seismic anomalies in the mantle wedge and crust.