Upper Plate Stress Controls the Distribution of Mariana Arc Volcanoes

We present a spatial analysis of volcano distribution and morphology in the young, intraoceanic Mariana Arc. Both the quality of fit to idealized models and the divergence from those ideals indicate that Mariana Arc volcanoes are arranged into five great circle segments, rather than a single small circle or multiple small circles. The alignment of magmatic centers suggests that magma transport is controlled by the stress regime in the deep crust and/or lithospheric mantle of the Philippine Sea Plate, into which the arc is emplaced, and that arc‐normal tension is the dominant process operating in the deep lithosphere along the whole arc. Volcano morphologies indicate that the stress regime in the shallow crust varies between arc‐normal tension and compression, which also implies that the stress field can vary with depth in the arc lithosphere. We show that this horizontal and vertical stress partitioning can be related to the changing dip of the subducting plate and the breadth of the zone where it is coupled with the overriding plate. The variation in stress regime is consistent with both the distribution of seismicity in the Philippine Sea Plate and with the structural fabrics of the nonvolcanic part of the plate margin to the south. Our analysis suggests that the upper plate exerts the principal control on the distribution of volcanoes in the Mariana Arc. Where tension in the deeper parts of arc lithosphere is sufficiently concentrated, then a distinct volcanic front is produced.


Introduction
Locations of volcanic edifices provide an opportunity to explore magma generation and transport beneath volcanic arcs (England & Katz, 2010). Processes occurring within the slab and mantle wedge, in particular through addition of fluid from the subducted slab to the wedge, make significant contributions to subduction zone magmas; therefore, some relationship between the locations of arc volcanoes and subduction dynamics can be anticipated England et al., 2004;Gill, 1981;Marsh, 1979;Syracuse & Abers, 2006;Tatsumi, 2005). England et al. (2004) approximated the distribution of arc volcanoes as small circles to define an average depth to slab (H) in each arc for comparison with other subduction dynamic parameters. This approach led them to propose that the locus of melting is related to the descent speed of slabs, from which it was inferred that the thermal structure of the mantle wedge is an additional key factor in localizing melting and, hence, volcano location (England et al., 2004;England & Katz, 2010). While the small circle approximation reveals correlated parameters for many arcs, the metamorphic reactions that release fluid from subducted slabs occur over a range of pressures and temperatures depending on the thermal and compositional profile of each slab (Grove et al., 2009;Schmidt & Poli, 1998). Thus, despite broadly similar H values within and between arcs (England & Katz, 2010;Jarrard, 1986;Wilson et al., 2014), with an average close to 105 km, it should be no surprise that there are wide variations, from about 60 to 200 km, in H within single subduction zones, often for volcanoes is close proximity (Pacey et al., 2013;Syracuse & Abers, 2006).
While mantle wedge hydration is clearly vital to the creation of arc volcano sources, an alternative view of volcano locations is that these are a function of structures or processes operating in the arc lithosphere, since the locations of volcanoes are the surface expression of magmatic pathways through the upper plate. The possibility that arc lithosphere might control arc volcano distribution was recognized from the early days of plate tectonics (Isacks et al., 1968) and substantial interplay of upper plate structural features with both spatial and temporal distributions of magmatism have been suggested in the Central American and Lesser Antilles margins (Bolge et al., 2009;Burkart & Self, 1985;Feuillet et al., 2002;Feuillet et al., 2010;Morgan et al., 2008;Weinberg, 1992). Pacey et al. (2013) demonstrated that volcanoes of the central Sunda Arc are aligned into a series of great circle segments, an arrangement that had previously been proposed to result from the locus of melt formation in many subduction systems (Marsh, 1979;Ranneft, 1979). In contrast to these previous studies, however, Pacey et al. (2013) attributed this arrangement to control of magma transport by the arc lithosphere because, with the exception of ocean island chains, most examples of great circle features on the Earth's surface are accepted as consequences of lithospheric control. The most notable examples of this are the products of tension, as seen to control alignment of magmatism at oceanic spreading centers and continental rifts, or association with lithospheric-scale fault systems, such as transverse and major normal faults. Therefore, recognition of great circle alignment in volcanic arcs may provide a means to determine the stress regime affecting arc lithosphere (Pacey et al., 2013). However, the margins mentioned above are predominantly continental systems, where structures inherited over protracted geological histories may also influence magmatic pathways.
To explore the distribution of arc volcanoes in an intraoceanic system, this paper examines the Mariana Arc (Figure 1), which initiated approximately 5 million years ago as the latest of several arcs to form in response to subduction of the Pacific Plate beneath the eastern margin of the Philippine Sea Plate (Fryer, 1996). The Mariana Arc covers geologically young, extensional basement and structures associated with rifting of the Mariana Ridge from the Mariana West Ridge (Bloomer et al., 1989;Hussong & Uyeda, 1981;Oakley et al., 2009;Yamazaki et al., 2003). Its youth, intraoceanic setting and, relatively, simple geological history mean that this margin is less susceptible to the structural and rheological complexities that influence continental arcs (Fryer, 1996). Thus, the Mariana Arc is well suited to understand the effect of current tectonic development upon arc volcano distributions. We demonstrate that the arrangement of Mariana Arc volcanoes is best described by a pattern of great circle segments. This segmentation is consistent with tensional forces dominating strain at the base of the arc lithosphere, thus focusing magma toward the volcanic arc. Comparison of these alignments with volcano ellipticity and seafloor fabrics indicates that stress is vertically partitioned in the arc lithosphere and that the nature of this partitioning varies along the arc. We develop a model for these stress variations, which is consistent with both earthquake focal mechanisms within the volcanic arc and the structural features of the nonvolcanic continuation of the Mariana margin to the southwest.

Mariana Arc
The intraoceanic Mariana Arc (Figure 1) is at the margin where subduction and volcanism have been ongoing since the Eocene (Hussong & Uyeda, 1981). Located at the eastern edge of the Philippine Sea Plate, where the Pacific Plate is subducted westward, the active Mariana Arc, which is dominated by basalt and basaltic andesite magmatism (Bloomer et al., 1989), has previously been described as comprising the Northern Seamount Province from 21°to 24°N, the Central Island Province from 16°(Anatahan) to 21°N (Uracas), and the Southern Seamount Province from 13°to 16°N (Dixon & Stern, 1983). Designation of boundaries between these provinces is arbitrary being based on volcano elevations with respect to sea level, which overlooks the presence of seamount volcanism within the Central Island Province. The whole arc, from Nikko in the north to Tracey in the south, comprises 60 groups of volcanic centers of which 26 are active. Twenty of the active centers are submarine (Baker et al., 2008). The active arc is located Figure 1. Distribution of Mariana Arc volcanoes on a gnomonic projection with center at 139.305°E, 17.333°N. Red triangles with black outline are subaerial volcanoes (Smithsonian's Institute Global Volcanism Program; our Data Set GVP04); those with no outline are submarine (abbreviations in the supporting information). Subaerial and submarine volcanoes between Farallon de Pajaros (FP) and Esmeralda (ES) are included in the second data set (GVP04+SM). All volcanoes from Nikko (NI) in the north to Tracey (TR) in the south comprise Data Set B08 (Baker et al., 2008). Digital elevation model in 30 arc second resolution (Becker et al., 2009) and additional 6 arc second resolution (Lim et al., 2013) from NOAA. Black dashed lines are backarc fracture zones. Thick red lines are sites of backarc spreading (Martinez & Taylor, 2003;Oakley et al., 2009;Yamazaki et al., 2003). Black lines with arrows are sinistral strike-slip faults (adapted from Stern & Smoot, 1998). Thin red to blue solid lines are slab contours in km beneath surface, from SLAB2.0 (Hayes et al., 2018). Green arrows illustrate backarc spreading directions inferred from motion of islands (Kato et al., 2003) with length proportional to rate. Blue arrows are the relative plate motion of the Pacific Plate (PA) to the stable Philippine Sea Plate (PS) with length proportional to rate that is annotated in mm/year (Argus et al., 2011). immediately to the west of Mariana Ridge, which is 100 km wide in the south and peaks at 300 m above sea level in Guam, narrowing to 20 km width near Sarigan where its crest is 1,100 m below sea level. This ridge forms a nonvolcanic island chain from Guam in the south to Asuncion in the north then becomes indistinct, further north (Bloomer et al., 1989; Figure 1). Hussong and Uyeda (1981) proposed that active edifices in the Mariana Arc are constructed on a basement of backarc crust or rifted arc crust that subsided along steeply dipping normal faults, to the west of the Mariana Ridge. Later, Bloomer et al. (1989) proposed that Mariana Arc edifices were aligned parallel to the West Mariana Ridge, Mariana Trough, and the forearc, from which they suggested that normal faults channel magma on its route to constructing the edifices. Accordingly, Bloomer et al. (1989) concluded that Mariana Arc volcanism is controlled by the structural development of the upper plate. Moreover, Oakley et al. (2009) interpreted scarps in the backarc, observed on seismic reflection profiles, as normal faults that formed by backarc spreading with most of the faults facing toward the spreading axes. Thus, it is also important to understand development of the backarc in order to understand the volcanism in the active arc.
The spreading-related, active faulting in the Mariana backarc is different to that in most ocean basins (Fryer, 1996), especially at fast or superfast spreading ridges where active faulting is concentrated within a narrow zone near the spreading axis (Edwards et al., 1991;Fornari et al., 2004). In contrast, diffuse extension of the Mariana Trough occurs as widely distributed normal faulting across the backarc even though a spreading axis is present and active (Fryer, 1995). Martínez et al. (1995) proposed that the diffuse deformation is caused by far-field strain as the upper plate deforms in response to subduction and described three types of structural development in the northern Mariana Trough: (1) asymmetric rifting between 22°15′N and 24°N, (2) localized rifting where spreading axis start to separate from the active arc between~21°N and 22°15′N, and (3) concentrated rifting from 20°N to 21°N to where the spreading axis is separated from active arc and forms deep grabens. This structural arrangement was later confirmed by identification of the rifting to spreading transition zone at about 22°N, through seafloor-spreading patterns in the bathymetry, magnetic field lineaments, and bulls-eye patterns in gravity data (Yamazaki et al., 2003). The same data also indicated diachronous initiation of spreading. Observations by Yamazaki et al. (2003) suggest that spreading between 19°N and 20°N started before 5 Ma then propagated to the north, which is compatible with the conclusion of Hussong and Uyeda (1981) who stated the spreading began after late Miocene.
Asymmetric spreading in the Mariana Trough has produced more backarc crust to the west of the spreading axis than to the east (Karig et al., 1978;Oakley et al., 2009). Yamazaki et al. (2003) suggested that this resulted from interaction of mantle upwelling beneath the active arc and the backarc spreading center. Deschamps and Fujiwara (2003) proposed it could be caused by preexisting magmatism in the east leading to asymmetry of crustal rheology, melting processes, and stress regime conditions, by resistance of the Pacific Plate slab in the eastern margin to the northwestward relative motion of the Philippine Sea Plate or by rollback of the Pacific Plate causing migration of the trench toward the east and southeast (Boutelier & Cruden, 2013;Faccenna et al., 2009). Slab rollback plays a significant role in backarc basin formation by causing hinge retreat and creating extension in the backarc (Macpherson & Hall, 1999;Macpherson & Hall, 2002); hence, rollback of the Pacific Plate caused trench retreat and began the opening of Mariana backarc basin (Faccenna et al., 2009).
The Mariana Trough spreading rate varies from north to the south. Based on modeling of magnetic anomalies and deep ocean drilling core analyses, the spreading half-rate to the west of the spreading center has varied from 2 to 3 cm/year since late Miocene (Hussong & Uyeda, 1981;Yamazaki et al., 2003). GPS observations conducted from 1991 to 1999 used the stable Eurasia reference frame for GPS stations in the island arc to determine the present backarc spreading rate relative to the Philippine Sea Plate. The backarc spreading rates are 15.9 ± 6.6 mm/year with an azimuth of 57.8°± 19.9°near Agrigan at the center of the arc, with a maximum rate of 44.6 ± 2.7 mm/year directed toward 97.1°± 4.1°near Guam in the south (Kato et al., 2003). This study also showed lateral, N-S motion in the residual data, which aligned with the model of spreading developing in the center of the trough propagating to the north and south (Martínez et al., 1995;Oakley et al., 2009).
Motion of the subducting Pacific Plate relative to the Philippine Sea Plate varies from north to south along the trench (Argus et al., 2011;Becker et al., 2015). At about 23°N the rate of convergence is 34 mm/year with an azimuth of N290°E (Argus et al., 2011). Convergence gradually decreases to the south, and the subduction direction rotates clockwise to 21 mm/year directed to N317°E at around 12°N ( Figure 1). Thus, as the trench azimuth rotates so does convergence obliquity from highly oblique in the north to orthogonal at about 13°N then oblique in the opposite sense further south. Stern and Smoot (1998) noted that this obliquity variation along the arc is manifest as the prevalence of left-lateral, strike-slip faults in the forearc north of 18°N contrasting with forearc grabens in the southern forearc ( Figure 1).
When averaged from 80 to 400 km depth, the Pacific Plate has a steep dip of approximately 75°beneath the Philippine Sea Plate (England et al., 2004), with some variation along the arc. Syracuse and Abers (2006) showed that between 50 and 250 km depth, the average dip decreases from 60°under Farallon de Pajaros (20.5°N) to 49°under the Esmeralda Bank (15°N). This is the same sense of variation as observed for depths greater than 125 km, where dips are 84°at 21°N compared to 73°at 12°N . However, at depths less than 125 km slab dips are steeper in the south, varying from 36°at about 21°N to 46°at 12°N (Fukao et al., 2001;Lallemand et al., 2005;Miller, Gorbatov, et al., 2006;. Variations in slab dip play no systematic role in variations in the depth from arc volcanoes to the slab, however. Figure 1 shows that in the southernmost arc there is an almost 60 km variation in depth to slab between Tracey and Northwest Rota seamounts. There is negligible variation in slab dip, or convergence rate and vector, between these volcanoes, which are separated by only 115 km. Gvirtzman and Stern (2004) used the term "plate coupling," which is different to "seismic coupling," to refer to the contact between upper and lower plates in a subduction zone. The steepening slab dip in the south, possibly related to tearing of the subducting slab and asthenospheric upwelling under the upper plate (Gvirtzman & Stern, 2004;Miller, Gorbatov, et al., 2006), also affects the zone in which such coupling occurs. The plate coupling zone is the horizontal breadth of this region of plate coupling at the surface, which they took to be represented by the distance from the trench to the front of the Mariana Ridge. This zone narrows southward from 170 km at 17.5°N, in the center of the arc, to 100 km at 11°N, in the south (Gvirtzman & Stern, 2004).

Data Sets
We employed three compilations of volcano locations to allow direct comparison of previous approaches with the method applied here (Figure 1). Following England et al. (2004), the first data set (GVP04), which largely equates to the subaerial Mariana Arc, takes the locations of 12 volcanoes from the Smithsonian Institute Global Volcanism Program catalog (Global Volcanism Program, 2013), from Esmeralda at the southern end to Farallon de Pajaros in the north. The second data set (GVP04+SM) includes all submarine volcanoes between Esmeralda and Farallon de Pajaros, which increases the number of edifices to 21. The third data set (B08) is taken from the study of Baker et al. (2008) on hydrothermal activity in the Mariana Arc. This includes 37 volcanoes from Nikko in the north to Tracey in the south, which geochemical studies suggest are all part of the volcanic arc (Baker et al., 2008;Pearce et al., 2005). The presence or absence of hydrothermal emissions (Baker et al., 2008) was used to constrain the locus of active volcanic craters. For centers lacking such emissions we used the least weathered crater morphologies to identify the most likely location of active or recent volcanism. The complete lists of volcanic center location data sets are listed in the supporting information.

Methods
We applied quantitative and objective methods to investigate the distribution of Mariana Arc volcanoes. Small circles were fitted to each data set (i.e., GVP04, GVP04+SM, and B08), and the Hough Transform method was applied to identify potential great circle alignments in B08 (Pacey et al., 2013). To determine whether a segmented small circle model might be more appropriate, best fit small circles were also obtained for the segments identified using the Hough Transform approach. Candidate great circle segments were compared to the structural lineament trends in the backarc basement and to the ellipticity of arc volcanoes (Nakamura, 1977). All results are integrated into a model that is evaluated against earthquake focal mechanisms and a structural understanding of the adjacent, amagmatic, southernmost part of the plate margin. A more detailed explanation of the methods is presented in the supporting information.
3.1. Geometric Fitting 3.1.1. Small Circles The geometry of a small circle is defined by the coordinates of its center (latitude and longitude) and its radius (r). The best fit small circles for the whole arc and for arc segments were obtained by varying those parameters to minimize the standard root-mean-square misfit (rms-misfit) of the volcanoes, as expressed by the equation where d n is the shortest (perpendicular) distance from each volcano to the small circle and n is the number of the volcanoes in the data set. This equation weights each volcano in the calculation equally, in contrast to alternative approaches that more heavily weight those volcanoes lying close to the small circle (e.g., England et al., 2004).

Great Circles: Hough Transform
To evaluate whether the Mariana Arc volcanoes are aligned as great circles, we used a Hough Transform approach (Ballard, 1981;Duda & Hart, 1972). This method has been used in Earth sciences to detect aligned structures, including volcanic vents and monogenetic volcanoes (Cebriá et al., 2011;Fernández-Álvarez et al., 2016;Martínez et al., 2000;von Veh & Németh, 2009;Wadge & Cross, 1988). Specific to arcs, Pacey et al. (2013) applied a Hough Transform approach to identify great circle segmentation in the central Sunda Arc. We have developed their method and applied this to determine potential alignment of Mariana Arc volcanoes. The quality of data fit to each potential great circle was, again, quantified by rmsmisfit. Initially, the endpoints of each segment were fixed as the locations of the two volcanoes at its ends (Pacey et al., 2013). Then, linear transformations were applied iteratively to the length, center point, and orientation of each potential alignment to minimize the misfit. Misfits for the overall arc were calculated for all possible segment combinations and weighted based on the number of volcanoes in each arc-segment. The best fit combination of great circle segments was then identified by minimizing the total number of segments, where possible associating each volcano with one segment only; maximizing the number of volcanoes on each segment; and minimizing the overall misfit.

Comparing Small and Great Circle Fits
We employed two approaches to compare the quality of fit of small and great circles. We used the Akaike Information Criterion (AIC) to compare the segmented small circle and segmented great circle cases because, although the number of segments are similar in each, the numbers of adjusted parameters (i.e., degrees of freedom) differ. A small circle has three degrees of freedom (radius, and central latitude and longitude) while a great circle has only two (central latitude and longitude). To allow direct comparison, we adapted the least squares fit case of the AIC (Akaike, 1974) where the model estimator is the rms-misfit that we obtained from fitting the geometric models (Banks & Joyner, 2017;Burnham & Anderson, 2004). The AIC is expressed as where n is the number of data points, K is the number of adjusted parameters, and b σ 2 is the estimator. Since we use the rms-misfit as model estimators, where d n is the residual or misfit from each volcano to its geometric model.
Given that the number of volcanoes in our data sets is small compared to the number of adjusted parameters in our geometric fitting, the AIC parameters should be corrected (corrected Akaike Information Criterion [AICc]) to prevent bias from the model with more adjusted parameters (Akpa & Unuabonah, 2011;Hurvich & Tsai, 1989). The AICc parameter formula is expressed as Since any AIC parameter represents the amount of information lost in fitting any model, models with the lowest AICc values should be preferred (Akaike, 1974;Burnham & Anderson, 2004).
A segmented model, whether of small or great circles, will involve considerably more degrees of freedom than for a single small circle; therefore, it is questionable how appropriate the AICc is for comparing a single small circle with a segmented great circle model. Instead, we analyzed the systematic changes in misfit (residuals) along the length of each segment with respect to ideal great and small circles. The principal behind this is outlined in Figure 2. Relative to a great circle datum, a chain of volcanoes that form a great circle on the surface of a sphere ( Figure 2a) would have residuals (ΔGC) that produce a y ≈ 0 regression line Figure 2. Use of residuals from fitting small and great circles to volcano distributions to determine true shape of arc segments. For volcanoes in great circle alignment (a) the residuals to a great circle (ΔGC) yield a linear regression with y ≈ 0 when plotted against distance along the segment (b), whereas the residuals for the best fit small circle (ΔSC) form a polynomial curve (c). The reverse is true for volcanoes distributed as a small circle (d) where ΔGC varies as a polynomial curve along the alignment (e), while ΔSC gives y ≈ 0 (f).
when plotted against distance along the segment (Figure 2b). Residuals for the same chain of volcanoes with respect to a small circle datum (ΔSC) would form a polynomial curve ( Figure 2c). Conversely, if the volcanoes are actually distributed as a small circle (Figure 2d), then the systematics of residual plots would be reversed; that is, ΔGC plotted against along-segment distance would generate a polynomial curve ( Figure 2e) and a line of y ≈ 0 would be produced for ΔSC (Figure 2f). The patterns predicted would be apparent regardless of whether a single small circle datum or segmented circle datums are used, although segmentation of either the data set or the datums would lead to inflections or reversals of the residual variation along the length of the arc.

Surface Strain Observation 3.2.1. Lineament Mapping
Extensional faults are distributed widely across the backarc and not concentrated within a narrow zone near the spreading axis (Fryer, 1996). We mapped individual normal faults, with the assumption that they are still in an early stage of development and can be considered as isolated structures (Cowie et al., 2000). Faults were identified on the bathymetric surface as lineaments across which the seafloor depth changes significantly. Backarc faults were assigned to a particular arc segment by projecting segment ends along the direction of backarc spreading (Kato et al., 2003). Nakamura (1977) proposed a method to approximate tectonic stress orientation from volcano morphologies by assuming that volcanoes experiencing a uniform, horizontal stress regime would produce a radially symmetrical central conduit and dike swarm in the shallow crust (Muller & Pollard, 1977). Using the Nakamura and Uyeda (1980) principle, that the vertical stress around the volcanic arc always forms the intermediate principal stress (σ 2 ), the two-dimensional near-surface stress regime perpendicular to the trench can be determined for any deviation from the ideal symmetrical arrangement toward an elliptical form. For a two-dimensional cross section perpendicular to the trench, trench-perpendicular compression would produce a maximum horizontal stress (σ Hmax ) parallel to the section. This, in turn, would produce volcanoes that are elongated perpendicular to the trend of the arc (Figures 3a and 3b). Conversely, trenchperpendicular tensional stress would lead to a minimum horizontal stress (σ Hmin ) perpendicular to the trench and cause volcanoes to be elongated parallel to the arc (Figures 3c and 3d).

Volcano Ellipticity
As proxies of elongation of the magmatic systems, Mariana Arc volcano ellipticity was determined for craters (Marliyani, 2016) or, in their absence, the footprint of the edifice (Bonini, 2012;Tibaldi, 1995). For some edifices the proximity of other volcanoes obscures the edifice shape; in which case, those volcanoes are not included in this analysis. Ellipticity was quantified by determining the azimuth and length of the longest and shortest axes, which are considered to approximate the orientation of horizontal maximum and minimum stresses (σ Hmax and σ Hmin ), respectively. Examples of volcano ellipticity observations are presented in the supporting information.

Model Evaluation 3.3.1. Focal Mechanism Distributions
Apperson (1991) employed earthquake focal mechanisms at shallow and intermediate depths in subduction zones to determine the seismic strain field of overriding plates. We applied the same method to focal mechanisms of earthquakes Mw ≥ 5.5 from the Global CMT catalog database (Dziewonski et al., 1981;Ekström et al., 2012) to evaluate the stress regime in the overriding plate. The analysis was carried out by projecting the individual focal mechanisms onto cross sections located at the center of each of the segments identified using the Hough Transform approach. Moment tensors were gathered in a swathe width equal to the corresponding linear segment. We focused on depths from 0 to 30 km where the relevant focal mechanisms appear.

Southern Mariana Structures
To further constrain the stress regime in the overriding plate, we studied the seafloor structures of the nonvolcanic, southwest continuation of the Mariana margin (Martinez et al., 2018). We extended the lineament map to the SSW of Tracey seamount into a part of the subduction system from 140°E to 146°E and from 11°N to 13°N that is referred to here as the Southern Mariana (Lim et al., 2013).

Geometric Fitting 4.1.1. Small Circles
The best fit small circle for GVP04 has an rms-misfit of 2.5 km, essentially identical to the 2 km value reported by England et al. (2004). Data set GVP04+SM, which shares geographical boundaries with GVP04 but includes submarine volcanoes, yielded a small circle with a higher misfit of 3.9 km. Data set B08, with geographical boundaries incorporating the whole Mariana Arc, gave the highest misfit of 8.4 km (Table 1). Thus, the effect of including all the edifices of data set B08 is a small circle with significantly larger misfit to the Mariana Arc than previously recognized.

Great Circles: Hough Transform
The Hough Transform method generated 12 potential great circle alignments within the Mariana Arc ( Figure 4a). These were optimized and combined, as discussed in section 3, to produce a refined pattern of five, best fit, great circle segments (Table 2 and Figure 4b): north segment (rms-misfit = 2.9 km), midnorth segment (rms-misfit = 3.2 km), central segment (rms-misfit = 2.7 km), midsouth segment (rms-misfit = 2.9 km), and south segment (rms-misfit = 0.7 km). The combined misfit values for all segments are 3.3 km before optimization and 2.7 km after optimization.

Comparing Small and Great Circle Fits
Once segments were identified from great circle fitting, their best fit small circles were also found. Comparison of the great and small circles for each segment with the AICc shows that segmented great circles are consistently better fitted to the arc than segmented small circles (Table 2). Treating the arc as a whole (supporting information), the AICc value for a series of great circle segments is 60.6, which is lower, and statistically preferable, to the AICc value for a series of small circle segments (89.9). Although we prefer a different method to compare the segmented great circle model with a single small circle (next paragraph), the AICc also returns a more favorable value for segmented great circles than for a single small circle (75.2). Therefore, despite the differences in the degrees of freedom, the Mariana Arc is better approximated as segments of great circles than segments of small circles.
Comparison of residual plots for fits to small circle (ΔSC) and great circle (ΔGC) datums also indicate that a segmented great circle model is preferable to either a single or segmented small circle model ( Figure 5).  (Marliyani, 2016;Nakamura & Uyeda, 1980). Vertical stress acts as the intermediate principal stress in the volcanic arc (σ v = σ 2 ), implying that the minimum and maximum principal stress are horizontal (σ Hmax = σ 1 and σ Hmin = σ 3 ). (a and b) Compression perpendicular to the trench produces maximum horizontal stress normal to arc trend, which is also the direction in which volcanoes are elongated. (c and d) Tension perpendicular to the trench leads to maximum horizontal stress that is parallel to the arc and volcanoes that are elongated along the trend of the arc.

Journal of Geophysical Research: Solid Earth
Residuals to the best fit, whole arc small circle for each data set do not show systematic variation for the whole arc ( Figure 5a). However, for shorter distance along the B08 data set, the small circle residuals show polynomial deviation. This is most evident in the south, from Tracey (TR) to West Saipan (WS), and the north, from Nikko (Ni) to South Daikoku (SD). Inflections or reversals in the residual variation occur    Figure 5. Residuals of small and great circle fitting plotted against distance along arc. (a) ΔSC of GVP04, GVP04+SM, and B08 data sets against distance from Tracey (TR). Full abbreviation list in Table S3. (b) ΔGC for Mariana Arc volcanoes from Data Set B08 using the best fit great circle alignments as the datum (Figure 2b). All segments plot as lines approaching y = 0; (c) ΔSC for Mariana Arc volcanoes from Data Set B08 using a single best fit small circle for the whole arc as the datum (Figure 2c; Table 1). (d) ΔSC for Mariana Arc volcanoes from Data Set B08 using the best fit small circle to each segment as the datum. All of the ΔGC and ΔSC residual plots suggest that Mariana Arc volcanoes follow great circle alignment.

10.1029/2019JB017391
Journal of Geophysical Research: Solid Earth close to the ends of segment identified using the Hough Transform. Figure 5c displays ΔSC residuals for the five segments identified using the Hough Transform method relative to the best fit, whole arc small circle. The residuals vary systematically along the length of each segment as anticipated for a distribution of volcanoes as a great circle. None of the residual plots in Figure 5c yield ΔSC distributions that can be better described as a straight line than a polynomial curve. This is corroborated by using a great circle datum for each segment, which produces straight lines that closely approach y = 0 (Figure 5b). This approach also replicates the result of the AICc approach by suggesting that Mariana Arc volcanoes are not arranged as multiple segments each constituting a different small circle (Figure 5d).
The analyses of distributions in this section lead us to reject the hypotheses that Mariana Arc volcanoes are distributed either as a single small circle arc or as segments made up of multiple small circles. In view of these observations, the simplest conclusion is that Mariana Arc volcanoes are distributed along five, great circle segments, which we shall refer to as arc segments. Not only is this most consistent with the Mariana Arc data, it is also most consistent with the Java arc system where a similar quantitative approach has been applied (Pacey et al., 2013), and to less quantitative analysis of multiple arcs worldwide (Marsh, 1979;Ranneft, 1979).

Surface Strain and Stress 4.2.1. Backarc Structure
Orientations of backarc faults associated with each arc segment ( Figure 6a) were compared to the backarc spreading direction, the segment azimuths, and the motion of the Pacific Plate relative to the Philippine Sea Plate. Rose diagrams show that fault orientations are perpendicular to the backarc spreading direction (Figure 7). This is consistent with the interpretation of Oakley et al. (2009) that the faults are associated with spreading. The arc-segment azimuths are subparallel to the fault orientations. In the northern segment the azimuths are slightly counter clockwise of the faults, and this offset gradually shifts to clockwise toward the south, but the deviation never exceeds 10°. In contrast, the subduction direction of the Pacific Plate does not show any systematic relationship to arc-segment azimuths.

Volcano Ellipticity
There are significant differences in the relationship between arc-segment azimuth and volcano ellipticity in the northern and southern segments (Figures 6b and 7). Volcano ellipticities in the north arc segment are elongated between N000°E and N080°E, the majority being strongly oblique to the segment azimuth (N316°E), the backarc faults, the convergence direction, and the trench. Similar, but stronger, relationships occur in the midnorth segment as the volcanoes are elongated between N040°E and N090°E while the segment trends toward N328°E. In contrast, the ellipticities of south arc-segment volcanoes are mostly elongated between N330°E and N020°E and, apart from West Saipan (WS), subparallel to the segment azimuth, which is oriented N031°E. Between these segments, the central and midsouth segments trend toward N350°E and N011°E, respectively. Volcano ellipticities in these two segments range between perpendicular to and parallel to the segment orientations, ranging from N340°E to N070°E in the central segment and from N045°E to N150°E in the midsouth segment.

Discussion
A relationship between averaged depths of subducted slabs beneath volcanoes (H) and descent speeds of slabs in some arcs has been used to infer a geodynamic significance for H (England et al., 2004;Syracuse & Abers, 2006). But several arcs, including Scotia, Java, Bali, Nicaragua, the Lesser Antilles, and Vanuatu, do not fit the global trend while other individual arcs show large fluctuations in H over distances where neither descent speed nor slab dip can vary (Syracuse & Abers, 2006). For example, in the central Sunda arc values of H change by more than 100 km, over horizontal separations of less than 150 km (Pacey et al., 2013). We have shown that Mariana is another example of this variation, since H for active arc volcanoes, determined from the SLAB2 model (Hayes et al., 2018), ranges from 110 to 180 km across the arc. This entire range is apparent over short distances in the south, from Tracey to Esmerelda, where there can be no change in slab descent speed (Figures 1 and S1). Furthermore, White et al. (2019) have suggested that earthquakes in the mantle wedge behind the Mariana Arc may reflect significant transport of slab derived fluid across a much greater breadth of the subduction margin than a simple, vertical conduit directly beneath the volcanic arc. This conclusion is consistent with modeling of fluid flow within the mantle wedge (Wilson et al., 2014). Both these approaches further illustrate that caution should be applied in relating H to the processes that localize volcanic and magmatic centers in arcs.
Several early studies of arcs that proposed that volcanoes form "linear"-in fact great circle-segments attributed segmentation to processes in the mantle wedge or to structures and processes at the slab surface Marsh, 1979; or were equivocal as to the causes of such distribution (e.g. Hughes et al., 1980;Ranneft, 1979). Pacey et al. (2013) documented great circle alignments of volcanoes in the central and eastern Sunda Arc but concluded that both the alignment and its segmentation into an en echelon arrangement resulted from magma transport being focused by the arc lithosphere. Close study of the Central American and Lesser Antilles systems has also demonstrated links between arc segmentation and structural development in the upper plate (Bolge et al., 2009;Burkart & Self, 1985;Feuillet et al., 2002;Morgan et al., 2008).
We propose that the great circle segmentation of the Mariana Arc reflects focusing of magma in the deep crust and/or lithospheric mantle of the upper plate, as has been suggested for the central Sunda Arc.

10.1029/2019JB017391
Journal of Geophysical Research: Solid Earth Pacey et al. (2013) inferred that great circle segmentation in central Sunda occurred due to magma exploiting stress-related, upper plate weaknesses and identified three principal mechanisms that might contribute to this: arc-normal tension, oblique tectonics, or upper plate flexure. In the remainder of this section, we evaluate these and potential alternative mechanisms in the Mariana Arc.
Alignment of contemporaneous volcanic, magmatic, and tectonic features in ocean ridge and rift settings is widely accepted to reflect magma channeling by the plates through which magma is transported (e.g., Crane & Ballard, 1981;Mazzarini, 2007;Rooney et al., 2011;Searle, 1992). Tension in the arc lithosphere might result from rollback of the subducting slab (Macpherson & Hall, 1999;Macpherson & Hall, 2002) and links between extension and volcanic productivity have been recognized in many arcs (Acocella & Funiciello, 2010;Smellie, 1995). The similarity of Mariana's arc-segment orientations to tension indicators in the immediate backarc (Figure 7) provides another indication that arc-normal tensional stress in the upper plate contributes to alignment of arc volcanoes. While tension, alone, would be consistent with the similar sense of deep and shallow lithospheric strain that we have determined from volcano alignment and volcano morphology, respectively, in the southern Mariana Arc, it cannot explain why the volcanoes of the north and midnorth segments have ellipticities consistent with arc-normal compression (Figure 7). Therefore, tension across the whole depth of the upper plate is unlikely to be the sole, or principal, cause of segmentation.

10.1029/2019JB017391
Journal of Geophysical Research: Solid Earth centers in the Lesser Antilles (Feuillet et al., 2002) and Central America (Bolge et al., 2009). Unlike the case in central Sunda, however, the Mariana arc-segments do not show the consistent en echelon step overs that led Pacey et al. (2013) to consider a role for oblique tectonics. Furthermore, the presence of arc segmentation is independent of changes in both the presence of strike-slip faulting in the Mariana forearc, and the overall convergence vector. Strike slip faults, oriented at about N314°E, cut the Mariana forearc north of 18°N where convergence becomes highly oblique (Stern & Smoot, 1998) but are absent further south (Figure 1), yet segmentation of the arc persists there. Indeed, the great circle segmentation persists along the entire length of the Mariana Arc as convergence varies from highly oblique in the north to entirely orthogonal at the latitude of the southernmost stratovolcanoes. Further south, convergence is again oblique but in a reversed sense, yet there is no strongly developed chain of arc volcanoes here (see below). Thus, we conclude that oblique convergence does not play a strong role in producing segmentation of volcanism in the Mariana Arc, although it may influence the extent of individual segments.
Downward flexure can produce compression in the shallow parts of a plate with simultaneous tension at greater depth (Hieronymus & Bercovici, 2000). Indeed, during the earliest days of plate tectonics arc-parallel tension at the base of arc lithosphere was considered as a dynamic response to flexure (Isacks et al., 1968). Husson (2006) modeled the development of negative dynamic topography in upper plates of subduction zones due to the presence of the dense, subjacent slab. In the Mariana margin this effect was predicted to produce maximum downward displacement of the Earth's surface in a broad, arc-parallel band between the Mariana Ridge and the backarc ridge axis (Husson, 2006), close to the location of the active arc. More generally, Hassani et al. (1997) demonstrated the development of lithosphere flexure at convergent margins due to the downward force from the slab as the hydrostatic suction maintained the coupling surface between the slab and upper plate. Crameri et al. (2017) showed that vertical deflection of the upper plate due to downgoing slab could occur up to a thousand kilometers from the trench; however, the wavelength of the flexure would depend on the upper plate's resistance to deformation based on its thickness and rheology (e.g., Meyer & Schellart, 2013;Sharples et al., 2014). In any case, horizontal tension at the base of upper plate would still occur regardless of the flexure wavelength as long as the coupling surface maintained the plate-to-plate contact and resistance (Hassani et al., 1997). In older arcs, flexure of the upper plate may be augmented by the load imposed by the arc, but this is unlikely to be a major effect in younger arcs such as Mariana (Waltham et al., 2008). In addition, Bonnardot et al. (2008) also predicted the possible contribution from mantle corner flow in enhancing the arc-normal tension at the base of flexed upper plate. Thus, the key feature to induce downward flexure of the upper plate is the operation of subduction, and, therefore, we consider this to be the most viable mechanism for producing arc-normal tension in the deeper parts of the upper plate along the length of the margin.
Other processes that may generate upper plate stress can also be evaluated at the Mariana margin. Lateral forces due to topographic or tectonic features may contribute to upper plate stress. For example, the topographic high of the Mariana Ridge (Bloomer et al., 1989) may exert tensional stress due to higher vertical loading or gravitational effects (Artyushkov, 1973;Bada et al., 2001;Bird, 1991). However, the form of the Mariana Ridge changes substantially along its length from around 100 km wide near Tracey seamount in the south, where it emerges as the island of Guam, to approximately 20 km width near Guguan, north of which it is not evident (Stern & Smoot, 1998; Figure 1). The presence of arc segmentation is not correlated with the presence or absence of the Mariana Ridge, or any other, topographic feature. Thus, we conclude that topographic effects have a negligible influence upon the development of segmentation in the Mariana Arc but are limited, perhaps, to the near surface stress regime.
To the west of the arc, backarc extension across a broad part of the Mariana Trough is suggested by the widespread distribution of normal faults (Fryer, 1995;Fryer, 1996;Martínez et al., 1995). A small ridge push effect from the spreading center, potentially close to zero due to the subdued topography in the backarc, may contribute to stress on the arc; however, a tensional stress regime induced as a passive response to the far-field rollback is likely to be dominant (Deschamps & Fujiwara, 2003;Macpherson & Hall, 2002;Nakakuki & Mura, 2013). Nonetheless, Mariana shares a great circle segmentation pattern with other arcs that lack backarc basins, including Sunda and Central America (Marsh, 1979;Pacey et al., 2013), suggesting that backarc spreading it not a primary control upon the development of arc segmentation.

Journal of Geophysical Research: Solid Earth
Consideration of all the factors above leads us to infer that arc-segmentation results from arc-normal tension in the deeper parts of the arc lithosphere, which is produced by tension and/or flexure of the plate upon which the Mariana Arc is constructed. This contrasts with a more complex variation of stress in the shallow upper plate, as indicated by volcano ellipticity (Figures 6 and 7) and mapped out in Figure 8 (Apperson, 1991;Nakamura & Uyeda, 1980;Oakley et al., 2009). In the north and midnorth arc segments, volcanoes are generally elongated subnormal to arc-segment azimuths (Figure 7), suggesting that horizontal maximum stress is perpendicular to the segments and that trench-perpendicular compression affects the shallow crust of the northern segment. While this is consistent with the general expectation of plate margin compression (Figure 3; Nakamura & Uyeda, 1980), it appears to contradict our interpretation of arc normal tension in the deep arc lithosphere. This contradiction would be resolved if the plate is flexing downward (Hassani et al., 1997;Hieronymus & Bercovici, 2001). In the south segment the volcano ellipticities suggest that horizontal maximum stress in the shallow arc crust is parallel to the arc segments (Figure 7), conforming to the conclusion reached for the deep crust but contrasting with the expectations from Nakamura and Uyeda (1980). In the central and midsouth segments, the orientation of the major axis of ellipticity is more variable suggesting a transition between trench-perpendicular compression and trench-perpendicular tension in the shallow crust.

Vertical Stress Partitioning
The apparent contradiction of inferred deep and shallow stress regimes in the lithosphere of the north and midnorth, and to a lesser extent in the central and midsouth, segments can be resolved by considering a model of vertical stress partitioning in the overriding plate. The gradual increase in vertical stress (σ v ), due to lithostatic pressure and rheological stratification, is recognized as an important control upon stress geometry and can potentially cause different stress regimes in the shallow and deep crust (Hasegawa et al., 1985;McGarr & Gay, 1978;Petrini & Podladchikov, 2000;Ranalli & Murphy, 1987). This increase in vertical stress means that the stress orientations approximated from volcano ellipticity should strictly be applied to the shallow lithosphere only and may not be the same as those in the deeper crust. In the shallow crust of an active arc, vertical stress acts as the intermediate principal stress (σ v = σ 2 ). However, increasing vertical stress, due to lithostatic stress, in the deeper crust might modify the principal stress geometry such that the vertical stress becomes the maximum principal stress (σ v = σ 1 ). Brace and Kohlstedt (1980) explained the mechanism of stress regime changes at depth using the limit of lithospheric strength, which is defined as the maximum difference between horizontal stress and vertical stress (σ H − σ v ). With a thermal gradient of 15°C/km in a crustal column with hydrostatic pore pressure σ H − σ v would reach a maximum value at about 15 km for a quartz rheology or 30 km for an olivine rheology, with shallower depths predicted for their dry equivalents. At deeper levels, σ H − σ v would decrease gradually toward zero at 25 km for the quartz rheology and 50 km for olivine. This change with depth provides a mechanism that could produce different stress regimes at different depths, even in the absence of forces external to the crustal volume. Therefore, it seems possible there is a mechanism that allows vertical partitioning of the stress regime in the upper plate, as most evident in the northern segments of the Mariana Arc. Furthermore, the relationships between shallow and deep stress appear to vary along the arc; they are aligned in the south but contradictory in the north. Thus, the stress regime must also be responding to changes along the length of the plate margin.
To further investigate vertical and horizontal stress partitioning in the Mariana margin, we investigated how the contrasting stress patterns vary with respect to along-arc changes in subduction dynamic parameters. The dip (δ) of the Pacific Plate beneath the Mariana Arc varies from south to north. Earthquake data show that the general dip of the subducting plate is less steep in the south (Syracuse & Abers, 2006) but at depths less than 125 km the slab dip (δ sh ) is actually steeper in the south (Fukao et al., 2001;Lallemand et al., 2005). In assessing plate coupling between the two plates, Gvirtzman and Stern (2004) focused on the shallow slab dip because this covers, and extends beyond, the depth range within which the Philippine Sea Plate and Pacific Plate are in contact. The average crustal thickness within the volcanic arc ranges from 14.5 to 20 km (Hughes & Mahood, 2011;Zellmer, 2008), and the average plate thickness is about 50 km (Gvirtzman & Stern, 2004). Therefore, δ sh is likely to be a more important influence than δ upon the stress regime of the overriding plate.

10.1029/2019JB017391
Journal of Geophysical Research: Solid Earth Figure 8. Trench-perpendicular stress regime in the shallow crust interpreted from volcano ellipticities, along with the backarc extension. Trench-perpendicular compression is marked by red shading, trench-perpendicular tension by yellow shading and transition zone by orange shading. Map is a gnomonic projection. Blue arrows show the PA-PS relative motion.

10.1029/2019JB017391
Journal of Geophysical Research: Solid Earth Gvirtzman and Stern (2004) estimated the extent of plate-to-plate coupling in the Mariana margin by measuring the horizontal, perpendicular distance from the trench to Mariana Ridge. The breadth of this plate-toplate coupling zone is about 170 km in the vicinity of our central arc segment and narrows southward to 100 km where we identify the south arc segment. Applying the Gvirtzman and Stern (2004) treatment to the margin at our north arc-segment suggests a coupling distance of up to 190 km. Assuming the thickness of the forearc is constant along the length of the arc and that coupling distances are defined in the horizontal plane, then the contact surface between the plates is approximated by the coupling distance divided by the cosine of the slab dip. Using these assumptions suggests that the contact surface between the plates would be wider in the north and narrower in the south.
Descent of the slab produces frictional resistance along the contact with the overriding plate. The downgoing slab movement itself is induced by a slab pull force (F sp ; Carlson et al., 1983). Gvirtzman and Stern (2004) proposed that the frictional resistance between the two plates and the downward movement of the subducting slab act to pull down the overriding plate. Since variations in this pull-down force (F pd ) are likely to be generated by changes in the surface resistance between the two plates, then F pd would be proportional to the breadth of the contact surface. Therefore, F pd would be stronger in the Mariana north arc segment and weaker in the south arc segment.
Variation of δ sh from north to south could also contribute further to modify the stress regime variation in the Mariana Arc. As the subducting slab descends below the overriding plate, a horizontal slab push (F up ) force is exerted upon the upper plate . Since velocity and force are proportional, F up will be proportional to the subducting velocity of the downgoing slab (v slab ) and the cosine of the shallow slab dip. Assuming that v slab , normal to the trench, does not vary significantly along the length of the arc, F up would be controlled by δ sh , and so the shallow slab dip in the north would generate stronger F up than the steeper dip in the south. Therefore, stronger compression would be exerted in the northern arc than toward the south.
Variation in the pull-down force may contribute to stronger upper plate flexure in the north compared to the south. In the north, vertical stress partitioning is most distinct where it may be enhanced by the strong compression from F up . Meanwhile, in the south, vertical stress partitioning is less obvious as the segment experiences stronger trench-perpendicular tension due to the seaward trench rollback and weaker F up . However, tension of deeper arc lithosphere occurs in all segments and is responsible for generating magma pathways in the deep crust of the Mariana Arc, which are manifest as arc segments.
Three-dimensional mechanical modeling of the strain mechanisms proposed above is beyond the scope of this study, but we can make a very simplified estimate of the possible magnitude of the operating stresses by assuming that the thickness of the upper plate (z) is a constant 50 km along the margin, comprising 20 km crust and 30 km lithospheric mantle (Gvirtzman & Stern, 2004). Then, the average shear stress (τ) operating along the plate interface is estimated as the product of the friction coefficient (μ), gravitational acceleration (g), density (ρ), and upper plate thickness (τ = μ.g.ρ.z). For this purpose, we assumed a typical friction coefficients (μ) for plate-to-plate coupling of 0.032 ± 0.006 for crust with a density of 2,800 kg/m 3 , and 0.019 ± 0.004 for lithospheric mantle with density of 3,300 kg/m 3 (Lamb, 2006). This results in an estimated shear stress along the plate interface of approximately 36 MPa. Although this is toward the upper end of estimates made for other subduction zones (e.g., Duarte et al., 2015;Holt et al., 2015;Lamb, 2006), it is not unreasonable in view of the limited constraints on real friction coefficients and on the thickness of lithosphere and crust. As the shear operates along the coupling surface, the horizontal and vertical stresses can be estimated as vector components of the mean value. The horizontal stress component, associated with F up , is the product of τ cos δ sh . Therefore, based on variation of δ sh along the margin, the horizontal stress is estimated to be about 29 MPa in the north and 25 MPa in the south. This is a small variation but one that is consistent with the changing forces we have proposed along the plate margin. The nonlithostatic component of vertical stress, arising from F pd , is τ sin δ sh leading to estimates of 21 MPa in the north and 25 MPa in the south. This decreases in vertical stress from north to south contrasts with the suggestion above but is, again, small and would be eliminated or even reversed if, for example, the arc lithosphere in the south were slightly thinner than and/or had proportionately more crust than lithosphere in the north. In general, basal traction exerted upon the lithosphere by the asthenosphere is thought to affect plates at a regional rather than at a local scale (Naliboff et al., 2009); therefore, we 10.1029/2019JB017391 Journal of Geophysical Research: Solid Earth do not expect plate scale forces of the sort described by Forte et al. (2010) and Ghosh et al. (2013) to be relevant to subduction zones. However, the more intensely focused flow of mantle generated by a subducting slab does have the potential to induce stress in the overriding plate (Bonnardot et al., 2008;Hassani et al., 1997). Quantification of such stress is highly dependent upon assumptions about the rheological properties of the two plates and the mantle wedge, but estimates are of comparable magnitude to those we have determined for the shear stress at the plate interface (Bonnardot et al., 2008). Figure 9 is a schematic illustration of how we interpret the stress regime of the overriding plate responds to changing dynamic parameters in the Mariana subduction system. The overall distribution of stress is similar to that envisaged due to plate flexure (Hieronymus & Bercovici, 2000). In the northern part of the arc the near-surface crust is dominated by trench-perpendicular compression that must give way to tension downward and to extension toward the backarc (Figure 9a). We envisage this being accommodated in a transition zone, which may have an intermediate stress state or may be composed of smaller domains of contrasting stress. In the southern arc ( Figure 9c) the trench-normal tensional regime is expressed as spreading in the backarc and tensional magma transport throughout most depths beneath the arc. The Mariana Ridge may represent compression of the upper plate with some component of bending at these latitudes, and there may also be some compression in the very shallow depths due to a volcano loading effect (Waltham et al., 2008). In the central parts of the arc the transition between tension and compression extends closer to the surface than in the north (Figure 9b). Like the north, the deeper crust beneath the central arc is under tension with the exact form of the volcano depending on the balance of tension and compression over the shallow depth range, where a loading effect may also come into play (Figure 9b). This geometry maintains a tensional regime at depth under all parts of the active arc allowing magma to ascend into the arc lithosphere where decreasing pressure and/or differentiation may aid further buoyant rise toward the surface regardless of the stress in the shallower crust.

Model Evaluation
Two independent sets of observations can be used to test our interpretation of stress distribution in Mariana Arc lithosphere. The first uses earthquake focal mechanisms as a direct measure of strain in the arc. The second uses the structures present in the southward, nonvolcanic continuation of the upper plate to explore the stresses that might be generated by subduction. Figure 10 shows trench-perpendicular cross sections for each of the five arc segments. Focal mechanisms of individual earthquakes are plotted on the sections, together with elevation profiles and our interpretation of stress regime boundaries from Figure 9. Figure 10c also shows the Moho profile of the central Mariana Arc from Takahashi et al. (2007). The location of these cross sections on the map view is in the supporting information, Figure S6.

Earthquake Focal Mechanisms
Most available focal mechanisms are concentrated at shallow and intermediate depths, from 10 to 22 km. In the north and midnorth segments, compressional focal mechanisms are located in the forearc, and extensional focal mechanisms occur in the backarc. Compressional mechanisms also occur near the trench, showing compressional stress fields both at the toe of the overriding plate and in the slab. In the south segment, extensional mechanisms are distributed widely from the backarc to the forearc, and compressional mechanisms are concentrated at the tip of the overriding plate. Between forearc and backarc compressional and extensional focal mechanisms are colocated in the south segment, and the arc lies within such a zone at the surface of the central and midsouth segments.
Seismic events are rare in the deep crust/shallow mantle zone beneath arcs. However, extensional events have been recorded at about 10 to 20 km depth beneath the north and midnorth arc segments (Figures 10a and 10b). In our model, this coincides with the deep expression of the great circle alignments, which define the arc segments that we infer to be the main supply routes for arc magma. Of the remaining events beneath the arc, most have an oblique sense but often with a strong normal component (Figures 10c  and 10d). The rarity of seismic events probably results from the rheology of the crust, its water content, and/or temperatures at that depth (Maggi et al., 2000;Meissner & Strehlau, 1982;Wortel, 1982) but could also be due to uncertainty in determining source locations of shallow to intermediate earthquakes (Ekström et al., 2012). In addition, dykes that intruding the deeper crust would involve tensile failure and so probably be aseismic ( , trench-perpendicular tension (yellow), and the transition zone between these (orange). The pull-down force (F pd ) on the overriding plate is controlled by plate coupling and shallow slab dip (δ sh ) where slab dip in the north is shallower than in the south (δ 1 < δ 2 < δ 3 ). Slab dip also controls the force to the upper plate (F up ) from the motion of the subducting slab. Plate flexure acts to cause tension in the deeper crust of the upper plate and focus magma flow into arc segments. There may also be a trench-perpendicular compression effect from volcano loading. The transition between the trench-perpendicular compression and tension may exist as subdomains of each type of stress or a gradual transition from compression to tension with depth. Dashed lines at the arc and inner trench wall show the upper and lower surfaces of the arc lithosphere before, or with less extreme, application of the subduction related forces. Figure 10. Individual earthquake focal mechanisms from the CMT catalog (Dziewonski et al., 1981;Ekström et al., 2012) projected on cross sections of each segment (adapted from Apperson, 1991) with bathymetric profiles. Moho depth details are only available in the central segment (Takahashi et al., 2007). Sections were plotted with GMT software (Wessel & Smith, 1998). Orange arrows describe the interpreted relative motion on the fault plane. Insets with gray background on sections (a) to (d) show the focal mechanisms in map view (map) and cross section (cs) to illustrate the strike-slip mechanisms. Abbreviations on insets describe the deformation type: ex = extension; cp = compression; ss = strike slip; sn = sinistral; dx = dextral; ob = oblique; nt = (strike) normal to trench.

10.1029/2019JB017391
Journal of Geophysical Research: Solid Earth mechanisms plots describe a strain pattern in the overriding plate that is consistent with our proposed model of stress regime fluctuation in the Mariana Arc lithosphere (Figure 9).
Combining our result with the findings of Apperson (1991) implies that our approach should reveal upper plate tension in other subduction zones. This is because Apperson (1991) concluded that subhorizontal, arc-normal tension was the principal stress regime in the overriding plates of most subduction zones, including Mariana. If we are correct to infer that the segmentation pattern of the Mariana Arc is controlled by variation in forces generated by the plate margin, then we predict that segmentation of other volcanic arcs will reflect changes in the forces generated at each margin.

Southern Mariana
Lineament analysis of the Southern Mariana area reveals short wavelength features in the bathymetry that are interpreted as normal faults (red lines on Figure 11) produced by extension perpendicular to the lineament orientations. South of the West Mariana Ridge the principal lineaments are oriented NE-SW at 140°E and ENE-WSW at 141°E, which is generally parallel to the trench and indicates NNW-SSE directed extension (σ 3 ). There are also a few lineaments near the trench that are oriented NNE-SSW, cross-cutting the major lineaments. The main lineament orientation remains ENE-WSW near 142°E and then deflects to NE-SW at about 143°E, which we infer means σ 3 in a NW-SE direction. The lineament orientation becomes parallel to the main NNE-SSW spreading axis near 144°E, consistent with σ 3 -oriented WNW-ESE.
Our study is consistent with the findings of Martínez et al. (2000) and Martinez et al. (2018), who showed that lineaments display contrasting orientation on either side of the spreading axis at its southern end. From the spreading axis to the West Mariana Ridge, as described above, most lineaments are parallel to spreading axis orientations while from the spreading axis to the Mariana Ridge lineaments are normal to the spreading axis and the trench (oriented in ENE-WSW and ESE-WNW directions). Ishihara et al. (2001) recognized these shorter lineaments at about 144°E as showing minor extension in a NNE-SSW direction, and GPS measurements also show a residual azimuth in N209°E direction near Guam (Kato et al., 2003). The trench-parallel lineaments formed during crustal accretion as a response to trench rollback in a southward direction such that these lineaments conform to the spreading direction and trench axis. Meanwhile, the trench-normal lineaments propagated to accommodate the increase of Mariana trench curvature (Martinez et al., 2018;Martínez et al., 2000) and thus appear as secondary structures on the surface. These variations of lineament orientation are part of diffuse extension as a result of weakening of the upper plate due to extensive hydration from the slab during early development of the subduction zone (Martinez et al., 2018). The diffuse extension in the Southern Mariana area developed rapidly and restricted the development of large, central arc-type volcanoes south of Tracey seamount (Stern et al., 2013).
Our lineament mapping in Southern Mariana (11-13°N) shows σ 3 mainly perpendicular to the trench, despite the relatively oblique convergence ( Figure 11). This is particularly evident between 140°N and 143.5°E consistent with our conclusion of trench-perpendicular tension in the adjacent (south) Mariana arc segment (Figures 7 and 8). A strong tensional regime is also consistent with the deep crustal tension implied from the volcano alignment of the south arc segment. Therefore, we infer that the trenchperpendicular tensional stress in the south arc-segment formed in a similar stress regime to structural lineaments of the Southern Mariana seafloor and was mainly caused by rollback of the subducting Pacific Plate. Martinez et al. (2018) have suggested that the absence of discrete central arc volcanoes here is due to the more diffuse focus of deep arc stress. This contrasts with the concentration of stress in the deep arc lithosphere between 13.5°and 23°N, which produces a more distinct arc front. Brounce et al. (2016) and Martinez et al. (2018) have postulated that the Fina Nagu volcanic complex, immediately west of the southern tip of the Mariana Ridge (Figure 11), may have been generated by the same mechanisms as arc volcanoes but that its elongate, multivent form may be due to focussing of magma by lithosphere deformation. Interestingly, the NE-SW arrangement of cones and vents of the Fina Nagu volcanic complex shows a similar degree of alignment to that of stratovolcanoes into arc segments that we have identified throughout the rest of the Mariana arc, albeit over a shorter distance, and is colinear with the principal lineaments immediately to its southwest ( Figure 11).
The tensional stress regime across Southern Mariana and the compressional stress regime at the plate margin can also be observed in the earthquake focal mechanisms. Figure 10f shows similar features to the south arc segment, where extensional focal mechanisms are widely distributed in the overriding plate while compressional mechanisms are concentrated at its trenchward tip. Strike-slip mechanisms are not observed from the spreading axis to the upper plate toe. This suggests that the stress regime depicted in our model is applicable to both the volcanic arc segments and the nonvolcanic region to the south west along the same plate margin.

Conclusions
Spatial analyses indicate that Mariana Arc volcanoes are distributed as a series of five great circle segments with lengths of 190 to 320 km. Two approaches demonstrate that the volcano locations are more consistent with great circle segments than any combination of small circles. First, the AICc values indicate that the great circle pattern (AICc = 60.6) is more consistent with volcano locations than either a series of small circles fit to multiple segments (AICc = 89.9) or a single small circle (AICc = 75.2). Second, deviation of volcano locations from best fit small and great circles show systematic variations consistent with the Mariana Arc volcanoes describing a segmented great circle distribution ( Figure 5).
The subduction direction and obliquity of convergence have no systematic relationship to the segmentation pattern. Arc-segment azimuths are parallel to the orientations of faults in their adjacent backarc and so to tension in the backarc. We infer that the arc segments are caused by tensional stress present in the deep arc lithosphere. This creates pathways that magma can exploit, leading to the alignment of magmatic pathways and the volcanic edifices that are their surface expression. Therefore, each arc segment represents a volume with relatively consistent tensional stress in the deep upper plate.
Elongation of volcanoes and/or their craters shows that the stress regime in the shallow crust of the north and midnorth arc segments is dominated by trench-perpendicular compression. We infer vertical and horizontal partitioning of stress in the arc lithosphere to explain both the apparent contradiction between indicators of deep (volcano alignment) and shallow (volcano elongation) stress in the north, and the change of shallow crustal stress orientation along the arc from north to south. Vertical and horizontal stress partitioning is a product of variations in dynamic subduction forces resulting from the coupling zone between the Pacific and Philippine Sea plates being wider in the north, where the slab dip is shallower, and narrower Figure 11. Lineament map of the Southern Mariana area to the southwest of Tracey seamount. Light blue arrows indicate the minimum horizontal stress (σ Hmin ) as the minimum principal stress (σ 3 ) direction identified from the mapped lineament while dark blue arrows indicate the maximum horizontal stress (σ Hmax ) as the maximum principal stress (σ 1 ) direction. Principal lineament orientations imply that extension occurs normal to the trench with secondary extension parallel to the trench to accommodate the strong curvature of the margin (Martinez et al., 2018). in the south. These differences reflect a stronger pull-down force on the overriding plate and stronger, trench-perpendicular compressional force in the shallow parts of the upper plate in the north and midnorth compared to the south segment. The pull-down force contributes to flexure of the overriding plate that produces tension in the lower crust beneath the arc (Figure 9). In the south, rollback of the Pacific Plate means that tension is prevalent throughout the full depth of the upper plate in the arc.