From intracrystalline distortion to plate motion: Unifying structural, kinematic, and textural analysis in heterogeneous shear zones through crystallographic orientation-dispersion methods

The localization of strain within the lithosphere is a fundamental process affecting the dynamic evolution of orogens. Within Earth’s crust and upper mantle, plastic heterogeneous deformation is commonly concentrated into shear zones—localized zones of deformation accommodated dominantly by ductile flow. Shear zones are significant in the tectonic evolution of many orogenic systems (e.g., Carreras, 2001; Dijkstra et al., 2002; Fusseis et al., 2006) because they are capable of accommodating large lithospheric displacements and are localized zones of high strain. In many crustal shear zones, strain is partly accommodated at the grain-scale by intracrystalline processes, such as dislocation creep; other mechanisms of strain accommodation include grain boundary migration, grain boundary sliding, and diffusion creep. Protracted deformation by dislocation creep progressively alters the texture and rheology of rocks (Brace and Kohlstedt, 1980), and these textural changes (e.g., grain size, grain shape, crystallographic preferred orientation [CPO]) can lead to grain size reduction by dynamic recrystallization (e.g., Handy et al., 2007; Platt and Behr, 2011) and the development of mylonites (e.g., Sibson, 1977; Lister and Snoke, 1984). During asymmetric shearing, mylonites develop microstructures that reflect the conditions of deformation and shear zone kinematics, including S-C structures, asymmetric porphyroclasts, and patterns of CPO (e.g., Poirier, 1980; Behrmann and Platt, 1982; Passchier and Simpson, 1986; Drury and Urai, 1990; Hirth and Tullis, 1992; Law et al., 2004). Microstructures and crystallographic textures of shear zone mylonites, coupled with field-based mapping and characterization, can therefore be used to recover information about the kinematics and conditions of lithospheric deformation in both active and exhumed orogenic systems.

Shearing of rocks at active tectonic boundaries induces material rotation at a range of scales, including the microscopic scale of the crystal lattice. For any scale of consideration, the degree of internal rotation can be represented by the mean kinematic vorticity number (W_m), which varies from 0 (perfectly coaxial deformation) to 1 (perfect simple shear). Traditional methods for estimating vorticity (e.g., Passchier, 1987; Simpson and De Paor, 1993, 1997; Wallis, 1992, 1995; Wallis et al., 1993) rely upon determining relationships between the shape (i.e., longest axis) and orientation of rigid particles, such as porphyroclasts, inferred to represent variable rotation in a flowing matrix (Jeffery, 1922). These methods provide a useful quantitative description of the relative contributions of pure and simple shear within a deforming zone, and they are best determined in the vorticity-normal plane. However, the vorticity-normal plane does not necessarily correspond to the plane parallel to lineation and perpendicular to foliation. In fact, during three-dimensional (3-D) deformation, the vorticity axis can be oriented perpendicular (e.g., simple shear), parallel (e.g., pure shear–dominated transpression; e.g., Tikoff and Greene, 1997), or oblique (e.g., triclinic shear zones; e.g., Jiang and Williams, 1998) to the lineation direction. This discrepancy highlights the fact that traditional methods of vorticity analysis depend on properly identifying the orientation of the vorticity axis, and that shape-based inferences about vorticity orientation (i.e., perpendicular to stretching lineation) are not always correct (Díaz-Azpiroz et al., 2018).

Recently developed methods for calculating shear-induced rotational axes from microscopic measurements in rocks collected from shear zones offer a possible solution to shortcomings inherent to traditional fabric-based methods for identifying vorticity axes. Crystallographic vorticity axis (CVA) analysis is a new method developed by Michels et al. (2015) that works by analyzing crystallographic distortion in the interiors of grains in sheared rocks. This method leverages rotational statistics to quantify an orientation-dispersion tensor of crystallographic orientations obtained using electron backscatter diffraction (EBSD). To date, CVA analysis has been applied primarily in monophase rocks, but it can also be used for polyphase rocks (e.g., Michels et al., 2015). CVA analysis is a unique method for studying deformational histories because it does not require assumptions about the significance of mineral shapes in determining the orientation of the vorticity tensor. Moreover, when combined with mineral shape information from field and microstructural studies, CVA analysis can be used to infer specifics about deformation histories, such as the nature of the deforming zone (e.g., monoclinic versus triclinic shear zones). Recent applications confirm that the results of CVA analyses are consistent with those of other modes of interpretation, even for shear zones that deviate from end-member simple shear (Michels et al., 2015; Giorgis et al., 2016; Schmidt et al., 2016).

The principal aim of this contribution is to combine field, microstructural, and textural analysis using EBSD to inform on the kinematics and conditions of deformation (i.e., temperature, stress, deformation mechanisms and kinematics) that facilitated strain localization along an enigmatic paleo–tectonic margin in eastern Massachusetts (northeastern USA)—the Burlington mylonite zone (BMZ). The BMZ is a NW-SW–trending, NW-dipping shear zone of greenschist- to amphibolite-facies metamorphic grade located entirely within the Boston Avalon terrane in eastern Massachusetts along the tectonic boundary with the Nashoba terrane (Fig. 1). The development of the BMZ is thought to represent the accretion of the Avalon terrane onto the Laurentian margin during the late Silurian–Early Devonian Acadian orogeny, but its timing and structural evolution remain largely unconstrained. This study integrates well-established field and microstructural analyses with CVA analysis in this polyphase, heterogenous shear zone as a means for recovering previously unknown aspects of the BMZ deformation history. Through this analysis we show that the BMZ can be dominantly characterized as a monoclinic shear zone that was active within the middle crust at or near the brittle-ductile transition. Additionally, this study further demonstrates the utility of CVA analysis in combination with traditional field and microstructural methods for deducing the sinistral transpressional deformation history within the BMZ.

Another innovative aspect of this study is that we extend our use of orientation statistics to the 3-D rock fabric itself. That is, for each measured field fabric, we define a single 3-D orientation that simultaneously represents a pair of lineation and foliation (pole) directions. This is a useful approach for our study of the BMZ—and likely for other fabric-based structural analyses of L-S tectonites—because it allows us to calculate a mean fabric orientation in a holistic manner, rather than by separately treating lineation directions and poles to foliation. Conversely, calculating directional means of lineation directions and foliation poles separately can result in two mean orientations that are not necessarily mutually perpendicular (i.e., the mean lineation does not lie within the mean foliation plane), despite the fact that the fabrics they are intended to represent are mutually perpendicular. Adopting a 3-D orientational representation of the fabric precludes such spurious results and is also useful for identifying outlier measurements in large data sets (e.g., Davis and Titus, 2017). Moreover, this approach permits us to calculate 3-D rotations between different analytical and data-acquisition reference frames so that we can compare results from microscopic analyses (acquired in a specimen’s fabric reference frame) to field measurements of fabric (acquired in a geographic reference frame).

The Burlington mylonite zone (BMZ) is a NE-SW–trending, NW-dipping shear zone in eastern Massachusetts that extends ∼50 km along strike from near the town of Westborough to near the town of North Reading. Named for its outcrops in Burlington, Massachusetts (Castle et al., 1976), the BMZ is demarcated by a heterogeneous shear zone of quartzofeldspathic (± hornblende-bearing) mylonitic gneisses and quartzites that lie entirely within the Boston Avalon terrane, adjacent to the present-day tectonic boundary between the peri-Gondwanan Nashoba and Avalon terranes (Fig. 1). The Boston Avalon terrane is the name for the eastern Massachusetts portion of the more extensive Avalon terrane in the northern Appalachian Mountains that has its type locality in eastern Newfoundland (Canada). The BMZ is the most prominent ductile shear zone in this section of the Avalon-Nashoba terrane boundary, which forms a continuous ∼500–1300 m zone of anastomosed and distributed mylonitic deformation within Avalonian rocks adjacent to, and east of, the terrane boundary (Kohut and Hepburn, 2004). At its greatest extent, deformation associated with the BMZ extends ∼5 km across strike. Rocks within the BMZ are recognized by a well-developed mylonitic foliation that strikes ∼192° to 257° and dips ∼35° to 75° to the northwest; lineation plunges shallowly (∼20°–35°) down to the northeast or southwest (Figs. 1, 2; Supplemental Table S1¹).

The protoliths of the BMZ mylonites include rocks belonging to the Westboro Formation and the Lexington plutonic suite, as well as undifferentiated Proterozoic and Paleozoic plutonic rocks (Bell and Alvord, 1976; Zen et al., 1983; Rast and Skehan, 1995; Kohut, 1999; Kohut and Hepburn, 2004). Metasedimentary rocks belonging to the late Proterozoic Westboro Formation comprise primarily light gray quartzites that are the oldest rocks in the Boston Avalon terrane. U-Pb detrital zircon ages constrain the maximum age of the Westboro Formation to 912 Ma; its minimum age is constrained by the intrusion of the Dedham granite (ca. 609–606 Ma; Thompson et al., 2012). Rocks in the western portion of the Boston Avalon terrane comprise variably mylonitized Proterozoic metamorphosed mafic to felsic volcanic rocks (e.g., light gray to pink granites and granodiorites) and mylonitized Paleozoic mafic plutonic rocks consisting of diorite and gabbro (Fig. 1) (Kohut, 1999; Brady and Cheney, 2004; Kohut and Hepburn, 2004). The Lexington plutonic suite is similarly composed of diorite and gabbro (U-Pb age of 427 ± 2 Ma) intruded by pink granites (Hepburn et al., 1998; Kohut and Hepburn, 2004). Unmylonitized gabbros and diorites of the Cambridge Reservoir suite (378 ± 3 Ma) cut early formed mylonitic fabrics (Hepburn et al., 1998). While the exact timing of Paleozoic mylonitization within the BMZ is unknown, it postdates the intrusions of the Lexington plutonic suite and predates the intrusions of the Cambridge Reservoir suite, constraining the timing of ductile shearing to between mid-Silurian and mid-Devonian (Kohut and Hepburn, 2004).

The BMZ is related in structural position to two younger zones of localized deformation: the Kendal Green mylonite zone (KGMZ) and the Bloody Bluff fault zone (BBFZ) (Fig. 1). The KGMZ is a NE-striking, NW-dipping shear zone located near of the town of Weston, Massachusetts, ∼3–4 km southeast of the Nashoba-Avalon terrane boundary (Kohut and Hepburn, 2004). The KGMZ is separated from the BMZ by undeformed rocks and is characterized by an ∼100–250-m-wide zone of protomylonite to ultramylonite with injection veins of pseudotachylite, contrasting the mylonitic fabrics in the BMZ and suggesting that it likely formed at shallower crustal levels and/or higher strain rate conditions (Kohut and Hepburn, 2004). The northwestern edge of the BMZ is overprinted by the BBFZ, which today forms the Boston Avalon–Nashoba terrane boundary (Fig. 1). The BBFZ is a zone of NW-dipping brittle faulting and cataclastic deformation that has been interpreted to follow the angle of the paleo–subduction zone (i.e., the original tectonic boundary as represented by the BMZ) that placed the Nashoba terrane over the Avalon terrane by early oblique sinistral thrusting in the mid-Silurian to mid-Devonian (Grimes and Skehan, 1995; Skehan, 1997). Late, brittle fault motion of the BBFZ has been variously identified as right lateral, left lateral, and/or having normal displacement (Goldstein, 1989; Grimes and Skehan, 1995; Kohut and Hepburn, 2004; Walsh et al., 2011). This subsequent reactivation of the Nashoba-Avalon tectonic boundary likely occurred in the late Silurian to Early Devonian, though could be as young as Mesozoic in age (Goldstein, 1991; Roden-Tice and Wintsch, 2007).

Systematic mapping of the BMZ identified 16 field sites between the towns of Westborough and North Reading, Massachusetts (in the southwestern and northeastern parts of the study area, respectively), with the majority of shear zone exposures located near the town of Burlington, Massachusetts (Fig. 1; Table S1 [^{footnote 1}]). At each locality, we measured macroscopic structural fabrics (e.g., orientation of foliation and lineation) and, where possible, determined relative shear sense from deformed structural markers and kinematic indicators (e.g., sigma and delta porphyroclasts, S-C-C′ fabrics).

Mylonites of the BMZ predominantly developed within quartzofeldspathic gneisses, granite, and quartzite belonging to the Boston Avalon terrane (Figs. 1 and 2). In outcrop, the BMZ mylonites occur as coarse-grained mylonitic rocks (Figs. 2A–2C) and fine-grained ultramylonites (Fig. 2D) containing metamorphic assemblages indicative of deformation under amphibolite-facies conditions based on the presence of hornblende. Solid-state deformation within the BMZ is evidenced by stretched quartz reflecting crystal-plastic deformation (Figs. 2B–2D, 2I) and abundant porphyroclasts of feldspar with asymmetric recrystallization tails (Figs. 2A–2C). Compositional banding is well developed in many of the BMZ mylonites and is variable in thickness (Fig. 2). In some localities, quartzofeldspathic amphibolite-bearing gneisses (e.g., sample BMZ-MP-07B) are cut by discordant granitic layers (e.g., sample BMZ-MP-07A) that are also deformed in the solid state (Fig. 2I). Similarly, veins of quartz and granitic dikes are folded and variably transposed within mylonites, reflecting variations in the finite strain magnitude throughout the shear zone. Kinematic indicators observed in the field and in thin section (Fig. 3) suggest a dominantly oblique sinistral shear sense in the BMZ mylonites, with the northwestern side (Nashoba terrane) moving to the south over the southeastern side (Avalon terrane), consistent with interpretations from past studies in the region (Goldstein, 1989; Kohut and Hepburn, 2004).

Thin section–scale observations of quartz show diagnostic microstructures indicative of dynamic recrystallization by subgrain rotation recrystallization and grain boundary migration (Figs. 3A–3F). Microstructures observed in thin section include (1) undulose extinction (Figs. 3B, 3F); (2) subgrain development (Fig. 3E); (3) oblique shape preferred orientations (SPOs); and (4) grain boundary migration (i.e., Fig. 3E) and dynamic recrystallization (Figs. 3C–3D) microstructures. Extensive recrystallization of quartz surrounding feldspar porphyroclasts with abundant intragrain cracks (Fig. 3C) and deformation tiling (Figs. 3D–3F) are indicative of deformation conditions at or near the brittle-ductile transition for feldspar. Subsequent deformation within the BMZ, but under lower-temperature conditions, is evidenced by the formation of overprinting localized chlorite- and epidote-bearing mylonite zones (Fig. 2F).

Twenty-four (24) samples were collected from the BMZ mylonites for microstructural and textural analyses, including 18 quartzofeldspathic gneisses, three quartzites, and three from transposed granitic layers. These samples reveal inter- and intracrystalline deformation features that allow for the qualitative and quantitative characterization of the conditions of deformation recorded in BMZ shear zone samples (e.g., temperature, stress history), determination of shear zone kinematics (e.g., from asymmetric porphyroclasts, CPO patterns), and deformation geometry using combined EBSD, CVA, and rigid grain net (RGN) methods. Polished thin sections were produced from oriented samples taken in the field. Thin sections were cut relative to the fabric reference frame such that they are oriented parallel to lineation and perpendicular to foliation (i.e., the X-Z plane of the fabric ellipsoid).

Microstructural analyses included both optical and electron microscopy methods. Petrographic analysis of thin sections using optical microscopy provide a critical first constraint on deformation based on the observed compositional and microstructural sample attributes. Full thin section backscattered electron (BSE) images were produced for all samples using a Tescan MIRA3 LMU Schottky field emission gun scanning electron microscope (FEG-SEM) housed within the Department of Earth and Environmental Sciences at Boston College (Chestnut Hill, Massachusetts). These BSE images provide detailed maps of all phase relationships and deformation microstructures in the samples that are in turn used to guide the selection of regions for subsequent EBSD textural analyses (Figs. 4, 5; Figs. S1–S3 [^{footnote 1}]) and the identification of porphyroclast phase boundaries critical to RGN analyses.

The RGN is composed of a set of semi-hyperbolas plotted in both positive and negative space at 0.025 intervals of W_m. These semi-hyperbolas transition into vertical lines to define the maximum B* as the critical threshold (R_c) between grains that rotate infinitely and those that reach a stable-sink position. In other words, grains with B* values below R_c rotate freely (i.e., they define a vertical boundary outside of the semi-hyperbola) and do not develop a preferred orientation, while grains with B* values above R_c have limited rotation (i.e., they lie within the semi-hyperbola). The mean kinematic vorticity number (W_m) is estimated directly from the RGN based on where the critical shape factor (B*) equals the aspect ratio of the porphyroclast (R; R = M_x/M_n) at the critical threshold (R_c). Therefore, each sample’s estimated R_c is equal to its estimated W_m.

The validity of RGN analysis relies on two assumptions: (1) the plane of reference for measurement is parallel to lineation and the pole to foliation; and (2) the plane of reference for measurement is normal to the vorticity axis. These assumptions define a monoclinic shear zone, and therefore require that the shear zone to which the methods are applied have a history dominated by monoclinic deformation. As mentioned earlier, this assumption is not always warranted for ancient shear zones due to the possibility of triclinic geological deformation that can result in oblique fabrics, or pure shear–dominated transpression that results in finite fabrics with lineation-parallel vorticity axes (e.g., Fossen and Tikoff, 1998). In the BMZ mylonites, a preponderance of field and microstructural evidence is compatible with a monoclinic deformation history, and therefore, the X-Z plane (i.e., the sample thin section surface) is initially assumed to be subparallel to the vorticity-normal section, satisfying the requirements of the RGN method. For the purposes of our study involving RGN analysis on two-dimensional sections in the X-Z plane, we therefore assume that W_m = W_k for samples dominated by monoclinic deformation. We revisit this assumption in subsequent sections with a follow-up quantitative and objective determination of the vorticity axis orientation using CVA analysis.

Electron backscatter diffraction (EBSD) is used to characterize microstructural relationships and quantify patterns of crystallographic orientations for all constituent mineral phases within the BMZ mylonites. EBSD analyses were performed in the Department of Earth and Environmental Sciences at Boston College using a Tescan MIRA3 LMU FEG-SEM equipped with an Oxford Instruments NordlysMax2 EBSD detector. EBSD data were acquired using an accelerating voltage of 25 kV and beam currents ranging from 50 to 100 nA. The analytical method is described in detail by Prior et al. (1999).

Large-area maps of crystallographic texture (Prior et al., 1999; Halfpenny et al., 2006; Maitland and Sitzman, 2007) were produced for all samples using Oxford Instruments AZtecHKL acquisition and analysis software (version 3.2). Areas were selected such that they would encompass the microstructural and grain size variations within each sample; as such, mapped regions vary in size from the scale of full thin sections (26 × 46 mm) to smaller regions of interest (Figs. 4, 5; Figs. S1–S3 [^{footnote 1}]). In order to ensure a high density of crystallographic orientation solutions within each grain, step sizes of 1.5 µm to 4 µm were used throughout the mapped regions. Indexing rates were high in all of the samples analyzed (typically >90%), and consequently complete EBSD data sets typically comprise >20 million individual solutions per sample map region.

EBSD data acquired using Oxford’s AZtecHKL software were subsequently processed using the freely available MTEX (version 4.5) toolbox for textural analysis in MATLAB (Bachmann et al., 2010). In order to ensure high data quality and account for potential indexing error, raw EBSD data sets were processed in MTEX to remove individual crystallographic orientations with median absolute deviation (MAD) values >1.0 (corresponding to a low confidence in the electron backscatter pattern solution). The rock microstructure of raw EBSD maps is revealed through a grain reconstruction algorithm that calculates the location of grain boundaries based on a misorientation angle of 10° relative to neighboring EBSD solutions. A grain set is calculated for all crystallographic orientations such that multiple individual crystallographic orientations are assigned to unique grains. Processed EBSD data sets from the BMZ samples typically comprise 5000–260,000 individual grains. Once the grain reconstruction is complete, a variety of microstructural information can be calculated from the grain set, including, but not limited to: (1) grain size data, (2) grain shape preferred orientation statistics, (3) mean grain-scale crystallographic orientations, (4) intragranular crystallographic misorientations, and (5) grain-scale crystallographic vorticity axes. For ease of comparison with experimental studies, and as required for subsequent paleostress calculations, grain size data are calculated as part of the MTEX grain reconstruction algorithm using the equivalent area diameter (2 × equivalent radius) method along with determination of the geometric mean of the full grain size distribution (Table S1; Figure S3 [^{footnote 1}]).

Pole figures of crystallographic orientations are plotted on lower-hemisphere stereonets where all constituent grains in the data set are represented as a unique single orientation on the stereonet by calculating the mean orientation of all indexed solutions within a grain (i.e., one point per grain is represented for all grains of all phases) (Fig. 7). Orientation distribution functions (ODFs) are a preferred method for quantitatively evaluating crystallographic texture and constructing hemispheric pole figures because they describe the 3-D frequency distribution of the measured crystal orientations in Euler space (Wenk and Wilde, 1972; Ismaïl and Mainprice, 1998). Calculations of ODFs are then interpreted with respect to the patterns of CPO that have been extensively documented in studies of experimentally and naturally deformed quartzofeldspathic rocks to assess deformation conditions and kinematics (e.g., Hirth and Tullis, 1992; Dunlap et al., 1997; Law et al., 2004; Toy et al., 2008). In all samples, CPO data are plotted in the fabric reference frame (i.e., relative to foliation and lineation). Finally, low-angle intragranular misorientation axes are calculated for quartz—the dominant phase in the BMZ mylonites—to aid in the determination of slip systems, CPO patterns, and deformation temperatures. Misorientation axes with angles between 2° and 10° are plotted on inverse pole figures as ODFs relative to the quartz crystal reference framework (Fig. 7). An example of the complete EBSD analytical method with associated calculated properties is provided in the Supplemental Material (Fig. S3 [^{footnote 1}]).

A recently developed method for analyzing EBSD data, crystallographic vorticity axis (CVA) analysis, which uses rotational statistics to analyze crystallographic orientations inside individual grains, has been shown to constrain the axis of kinematic vorticity within plastically deformed sheared rocks (Michels et al., 2015). CVA analysis differs from misorientation axis analysis in that CVA analysis calculates a single best-fit rotation axis from numerous crystallographic orientations within a grain, rather than comparing individual pairs of orientations (as in the calculation of a misorientation axis). The main advantage of CVA analysis is that it can be used to objectively determine a grain-scale lattice rotation tensor (an axis and a magnitude of rotation) that describes the rotational dispersion of mutually dependent crystallographic axes. Importantly, a single grain’s lattice rotation axis is insufficient to deduce larger-scale rotation axes because each grain is thought to respond to very localized stress conditions in a deforming polycrystalline aggregate. However, the dominant direction of numerous grain-scale dispersion axes within the sample provides a robust proxy for the rotational axis of the sample-scale deformation geometry (i.e., the axis of kinematic vorticity).

It is important to clarify what a “grain” represents in a sheared rock’s deformation history because most “grains” likely represent ephemeral features within a long-lived shear zone, such as the BMZ, although they define the entirety of the microstructural record of strain. In this context, a grain can be understood to represent an intact crystalline volume that has a crystal lattice orientation distinct from grains surrounding it. In other words, a grain is identified as different from surrounding grains based on a dominant crystal lattice orientation.

Grains in rocks that have been deformed typically exhibit evidence of internal lattice distortion that preserves a partial record of dislocation-accommodated strain at the intragranular scale (e.g., undulose extinction, subgrain development). Using an approach that is conceptually analogous to principal component analysis (principal geodesic analysis; Fletcher et al., 2004), CVA analysis is capable of calculating a rotational axis that best fits the 3-D distortion of the crystal lattice within an individual grain (Michels et al., 2015). Such a best-fit axis is presumed to record the orientation of a local, grain-scale vorticity axis. Because most grains are likely to be ephemeral features in a rock throughout a protracted deformation history, it is assumed that the lattice rotation axis (CVA) is reset relatively rapidly compared to timescales of geological deformations. Therefore, these axes are considered to represent a snapshot of effectively instantaneous rotation in each grain, and thought to likely record the most recent geological deformation of the crystal lattice in each grain. Furthermore, it is presumed that the sum of such instantaneous grain-scale deformations reflects a snapshot of the geometry of the larger-scale deformation, such that by calculating the lattice rotation axes within many grains in a sample, a dominant rotational axis can be statistically identified for the rock. This approach has been shown to be capable of tracking the deformation geometry in geological shear zones that are characterized by simple shear, subsimple shear, and transpression (Michels et al., 2015).

CVA analysis is made possible by processing EBSD data using a combination of custom code and the MTEX (version 4.5) software toolbox for textural analysis (Bachmann et al., 2010). Open-source code for CVA analysis, developed by Michels et al. (2015), is available via download from the CVA repository on GitHub (https://github.com/zmichels/CVA). The code is written for the MATLAB programming platform and designed to integrate directly with the MTEX toolbox for microstructural analysis in MATLAB.

Individual grain-scale orientation-dispersion axes are calculated for every grain of all analyzed phases in the reconstructed EBSD data sets. Due to differences in grain sizes in rocks, each thin section yields a different number of grains for analysis. As such, the resultant populations of dispersion axes vary in size from sample to sample. The number of axes calculated in individual samples ranges from ∼150 to >15,000, with typical samples yielding >1000 dispersion axes.

A preferred dispersion axis is identified from each sample-scale population using kernel density estimation with a de la Vallée Poussin kernel and a 10° halfwidth. This procedure calculates the directional orientation density (i.e., multiples of uniform density) of a population of axes relative to an orientation grid space. A grid with a 1° resolution is used to calculate a single, antipodal axis that matches the greatest orientation density for each sample. This approach follows methods of sample-scale CVA analysis originally outlined by Michels et al. (2015) and subsequently applied in other recent studies (e.g., Giorgis et al., 2016; Schmidt et al., 2016).

Kinematic indicators are assessed both in the field and in thin sections cut parallel to lineation and perpendicular to foliation from oriented samples. A number of macroscopic and microscopic shear-sense indicators are present in the BMZ mylonites, including: (1) sigma and delta feldspar and hornblende porphyroclasts with asymmetric recrystallized tails (e.g., Figs. 3A–3C, 3E–3F); (2) low-angle and high-angle tiled feldspar (Figs. 3D–3F) and hornblende porphyroclasts (Figs. 3A and 4A, 4B, 4G); (3) shear band fabrics, such as S-C mylonitic fabrics (Fig. 5D); (4) en echelon veins (Fig. 3C); and (5) microfaulting (Fig. 4E). These structural fabrics exhibit dominant top-to-the-south, sinistral shear sense at both the macro- and microscopic scales of observation. However, mixed kinematic indicators (e.g., sinistral and dextral porphyroclasts) are present in all samples (Fig. 3A), indicating a substantial component of non-coaxial flow within the shear zone.

Ten samples were chosen for vorticity analysis using RGN methods. Full thin section polarized-light and BSE maps were spatially correlated such that porphyroclasts within each sample could be identified based on their shape and orientation relative to the plane of mylonitic foliation. As necessitated by rigid grain methods for vorticity analysis (e.g., Passchier, 1987; Simpson and De Paor, 1993, 1997; Wallis, 1992, 1995; Wallis et al., 1993; Jessup et al., 2007), porphyroclasts of feldspar and hornblende were treated as tailless and chosen based on lack of proximity to other porphyroclasts to ensure they were freely rotating within the quartz-rich matrix during deformation. Grain shapes for porphyroclasts suitable for RGN analysis were manually digitized by comparison of correlated full thin section optical and BSE images. Vectorized digital files of porphyroclasts shapes were, in turn, processed using a custom MATLAB script to calculate the long and short axes of each porphyroclast and the associated angle made between its long axis and the plane of mylonitic foliation. At least 200 clasts were characterized within each sample, except in one sample (BMZ-MP-15A-1) in which only 63 grains met the necessary criteria (Fig. 6). Sample data were plotted on the RGN (Jessup et al., 2007) in order to estimate the mean kinematic vorticity number (W_m). In all ten samples, the maximum B* value was identified as the critical threshold (R_c) between infinite rotation and the stable-sink position, providing an estimate for the kinematic vorticity. The estimated kinematic vorticity number within each of ten RGN-analyzed BMZ samples ranges from 0.25 to 0.55 (Fig. 6).

Deformation within shear zones accommodated by dislocation creep can result in the formation of a CPO pattern. The framework for interpreting patterns of CPO for quartz-rich mylonites is well established by experimental studies (e.g., Hirth and Tullis, 1992; Heilbronner and Tullis, 2002, 2006) and field studies of naturally deformed terrains (e.g., Behrmann and Platt, 1982; Dunlap et al., 1997; Law et al., 2004). Quartz—the dominant mineral by mode in all of the samples analyzed—has been shown to develop distinct CPO patterns based on the kinematics of deformation (e.g., sinistral versus dextral shear sense), temperature of deformation (activation of distinct slip systems: e.g., prism<a>, rhomb<a>, prism[c] slip), and strain geometry (coaxial versus non-coaxial deformation; cf. Passchier and Trouw, 2005). Therefore, we focus on interpreting CPO patterns in quartz to constrain the conditions of deformation within the BMZ.

Pole figure distributions of crystallographic orientations for quartz are plotted with respect to the macroscopic mylonitic foliation and lineation as determined in the field (Fig. 7). Quartz CPO patterns are highly variable within the BMZ mylonites, as recorded by c-axis distributions ranging from: (1) point maxima oriented within the plane of foliation, perpendicular to lineation (e.g., samples BMZ-MP-11A-1, BMZ-MP-20A-3; Fig. 7); (2) girdled distributions of various strength (e.g., samples BMZ-MP-04B-1, BMZ-MP-06A-2, BMZ-08A-1; Fig. 7); and (3) point maxima oriented perpendicular to the plane of foliation (e.g., sample BMZ-MP-PilgPath01; Fig. 7). Of these varieties, quartz c-axis distributions are primarily concentrated within the plane of the foliation perpendicular to lineation with corresponding a-axes forming weak point maxima roughly parallel to lineation or at high angle to the foliation—patterns broadly consistent with those expected for the activation of basal<a> and rhomb<a> slip or prism<a> slip in quartz (Fig. 7; cf. Schmid and Casey, 1986; Stipp et al., 2002; Passchier and Trouw, 2005; Toy et al., 2008; Barth et al., 2010).

To assess the strength of the CPOs, we evaluate two different texture indices, the J-index and M-index, that provide a relative measure of CPO intensity. The J-index ranges from a value of one (a completely random distribution) to infinity (a single crystal), whereas the M-index ranges from zero (a completely random distribution) to one (a single crystal) (Bunge, 1982; Barth et al., 2010; Mainprice et al., 2004; Skemer et al., 2005). The one-point-per-grain J-indices range from 1.2 to 4.1, and the one-point-per-grain M-indices range from 0.01 to 0.28 (Fig. 7).

In deformed grains, low-angle (2°–10°) boundaries develop during dislocation creep as a result of dislocation organization into planes of lower energy (Lloyd et al., 1997). The organization of edge dislocations produces tilt boundaries, while organization of screw dislocations forms twist boundaries. The identification of low-angle crystallographic misorientations allows for a more robust determination of active slip systems than being informed by patterns of CPO alone, and thus provides better constraints on the temperatures of deformation that accompany their differential activation.

Consequently, we identified all low-angle 2°–10° neighbor-to-neighbor misorientations in quartz grains. The misorientation axis plots for quartz are shown in Figure 7 along with their corresponding CPO patterns. All 24 BMZ samples preserve quartz misorientation axes with distributions that cluster near the c-axis, consistent with a predominance of prism<a> slip (Neumann, 2000). Several samples display secondary local maxima along the outer edge of misorientation inverse pole figures (e.g., samples BMZ-MP-04A-1, BMZ-MP-09A-1; Fig. 7). This dual-maxima pattern could result from: (1) a combination of prism<a>, basal<a>, and rhomb<a> slip-system activation; (2) activation of prism<a> and prism[c] slip systems (e.g., Lloyd et al., 1997; Neumann, 2000); or (3) a combination of basal<a>, rhomb<a>, prism<a>, and prism[c] slip-system activity. Although deformation conditions in some of the BMZ mylonites may have reached temperatures high enough for activation of the prism[c] slip system (e.g., 550–600 °C; cf. Toy et al., 2008, and references therein), quartz CPO patterns in the same samples do not exhibit c-axis maxima parallel to the shear direction, as is expected for a significant contribution of prism[c] slip during deformation. As such, we interpret that these patterns reflect a combination of activity of prism<a>, basal<a>, and rhomb<a> slip systems.

Using the CVA method, a single axis that represents a particular grain’s vorticity axis is determined from the 3-D dispersion of crystallographic orientations within the grain; this calculation, in turn, is repeated for all grains within each EBSD data set. The critical thesis of Michels et al. (2015) is that the sum of all of the dispersion axes within the sample is correlative with the vorticity axis of the aggregate-scale deformation geometry. While this method was originally tested in quartzites, conceptually there is no limitation for its application in other minerals and polyphase lithologies (Michels et al., 2015) such as those found in the BMZ mylonites. Consequently, the results of the CVA analysis are presented as pole figures of dispersion axes in all phases within a given sample (Fig. 8) and, for comparison, as phase-specific pole figures for individual constituent phases (Fig. 9; Fig. S4 [^{footnote 1}]).

CVA plots for each rock sample (i.e., bulk CVA; all grains, all phases analyzed together) are remarkably consistent throughout the BMZ. The bulk vorticity axes calculated for 23 of the BMZ samples plot within the plane of foliation and perpendicular to lineation, as indicated by consistent “bullseye” distributions in the center of the stereonets in the sample (i.e., fabric) reference frame (Fig. 8). This particular CVA pattern correlates to the orientation of the vorticity axis expected for simple shear and wrench-dominated transpression involving general shear (Michels et al., 2015, and references therein). Some variation of the bulk CVA centered point maxima pattern is evident in the form of slightly oblique or girdled maxima. The calculated bulk vorticity axis for sample BMZ-MP-PilgPath01—a quartzite—is the only one that plots as point maxima oriented normal to the foliation plane (Fig. 8). Notably, Figure 7 shows that this same sample is characterized by slip associated with rotation about the c-axis and a CPO with a c-axis point maxima oriented at high angle to the foliation plane, consistent with the CVA results. All bulk CVA pole figure distributions are also shown in their geographic reference frame in Figure 10, such that they are rotated with respect to the orientation of the macroscopic foliation and lineation observed in the field.

The major phases identified within the BMZ samples include quartz, anorthite, orthoclase, and (where present) hornblende. CVA plots for each of the major phases are remarkably consistent with the bulk CVA plots. Phase-specific CVA vorticity axes similarly plot within the foliation plane, displaying a bullseye point maxima perpendicular to lineation (i.e., parallel to the y-axis of the fabric ellipsoid) (Fig. 9; Fig. S4 [^{footnote 1}]). CVA analysis plots for quartz most strongly resemble the bulk CVA axis plots, with the exception of sample BMZ-MP-16A-1, where the anorthite CVA pattern more strongly resembles the bulk CVA pattern. Among the dominant phases, the orthoclase CVA patterns show the most variation from the bulk CVA patterns. These variations likely reflect the relative contribution of distinct phases to the bulk CVA based primarily on the mode, connectivity, and strength contrast of phases.

For the case of W_m = 0.25, we calculate an oblique flow apophysis inclined 75.5° from the foliation plane and instantaneous stretching axes inclined 7.2° (ISA₁) and 97.2° (ISA₃) from the lineation direction (Fig. 11B). For the case of W_m = 0.55, we calculate an oblique flow apophysis inclined 56.6° from the foliation plane and instantaneous stretching axes inclined 16.7° (ISA₁) and 106.7° (ISA₃) from the lineation direction (Fig. 11C). For each case, the other two flow apophyses of general shear are necessarily aligned with the lineation and the vorticity axis directions (Fig. 11D).

We combine the directions of the stretching lineation and the foliation pole from each field-based measurement to formulate a set of 3-D orientational representations of BMZ rock fabrics in the geographic reference frame (e.g., Davis and Titus, 2017). By doing so, we can define a rotation that maps results from the specimen reference frame to a unifying geographic reference frame. Using this approach, we rotate directional data from the results of our microanalyses (i.e., CVA axes, ISAs, flow apophyses)—which we calculate in separate specimen reference frames defined relative to the rock fabric in each sample—into a common geographic reference frame. Furthermore, this approach allows us to perform a variety of cursory orientation statistics within the geographic reference frame, such as calculate a dominant rock fabric orientation from the mode of an ODF calculated from 3-D fabric orientations—one that maintains the mutually perpendicular character of paired lineations and foliation poles. The dominant fabric orientation comprises a foliation plane that strikes 217.6° and dips 74.8° NW (Fig. 11E) and a lineation that trends 32.7° and plunges 17.4° (Fig. 11F). We use this dominant rock fabric orientation to present geographically contextualized versions of the calculated flow apophyses and ISAs, and thereby constrain a feasible tectonic-scale interpretation of the deformation geometry for the BMZ.

We also rotate vorticity axes from the CVA analyses based on each sample’s geographic orientation of rock fabric, as illustrated in Figure 10. In the geographic reference frame, we conduct a kernel density estimation of the axial orientations using a de la Vallée Poussin kernel with a halfwidth of 15°. From the density estimation, we calculate a single bulk vorticity axis for all samples to characterize the dominant deformation geometry of the BMZ at the tectonic scale (Fig. 11G). The geographic orientation of the derived vorticity axis trends 280° and plunges 69°. For this study, the main utility of the bulk geographic vorticity axis is for further constraining the geographic relevance of ISAs and flow apophyses derived from the results of RGN analysis (Fig. 11H), as we discuss in a following section.

Field, optical, and electron microscopy methods were used in this study for the purpose of determining the conditions of deformation preserved within ductilely deformed rocks comprising the Burlington mylonite zone (BMZ). Previous studies of the BMZ have focused primarily on the age and kinematics of deformation and have rarely addressed the conditions of deformation as revealed through quantitative microstructural methods (e.g., Castle et al., 1976: Rast and Skehan, 1995; Skehan, 1997; Kohut and Hepburn, 2004). In this study, we constrain temperature of deformation based on metamorphic assemblages, deformation microstructures, and crystallographic textures. We also constrain differential stresses within the BMZ mylonites by applying the recrystallized grain size piezometer for quartz.

All of the analyzed rocks of the BMZ are quartzofeldspathic gneisses, with the exception of volumetrically minor quartzites and granitic intrusions. Field observations of the BMZ, like in previous studies (Castle et al., 1976: Rast and Skehan, 1995; Kohut and Hepburn, 2004), suggest a complex history of mylonitic deformation followed by structural overprinting related to brittle deformation associated with the formation of the Bloody Bluff fault zone (BBFZ). Deformation under amphibolite-facies conditions is readily recognizable in the field by the presence of variably mylonitized (proto- to ultramylonite) hornblende-bearing quartzofeldspathic gneisses (Figs. 2–4). Subsequent shearing under lower-temperature conditions is evidenced by retrograde metamorphism involving the formation of chlorite- and epidote-bearing crosscutting shear zones (Fig. 2F). Given that the formation of the BBFZ represents a transition to brittle (cataclastic) deformation during subsequent structural overprinting at lower metamorphic grades, our discussion herein targets the significance of the high-temperature ductile history of the BMZ—the focus of this study.

Microstructures in the BMZ mylonites are indicative of subgrain rotation recrystallization and grain boundary migration dynamic recrystallization mechanisms (Fig. 3), consistent with the deformation conditions inferred from past studies and the observed metamorphic assemblages. It is important to note that these rock microstructures may not necessarily record the peak conditions of the BMZ mylonites (cf. Toy et al., 2008), but they provide a reasonable minimum estimate of the high-temperature conditions of ductile deformation.

Quartz CPO patterns are broadly consistent with experimental patterns expected for the activation of prism<a> slip based on c-axis maxima oriented within the plane of the foliation perpendicular to lineation, whereas girdled variations are also suggestive of the activation of basal<a> and rhomb<a> slip systems in quartz (Schmid and Casey, 1986) (Fig. 7). Further, the CPO pattern of the BMZ samples is as expected for wrench deformation geometries involving a component of flattening (as evidenced by the experimentally derived quartz CPO Flinn diagram; cf. Lister and Hobbs, 1980). Misorientation axis plots from all samples indicate dominant prism<a> slip in the BMZ mylonites, with some samples also indicating the activation of basal<a>, rhomb<a>, and—potentially—prism[c] slip. The activation of these slip systems agrees with the range of temperatures inferred from mineral parageneses and rock microstructures, suggesting that deformation occurred at moderate temperatures (∼450–600 °C) (Toy et al., 2008, and references therein). This inference is further supported by microstructural observations showing that feldspar in the BMZ mylonites predominantly deformed via fracturing, rather than by pervasive intracrystalline slip and thermally activated creep (Wilks and Carter, 1990; Stünitz et al., 2003; Rybacki and Dresen, 2000; Mehl and Hirth, 2008).

Finally, the piezometric relationship calculated for quartz using EBSD-generated grain size data indicates paleostress conditions in the BMZ ranging from ∼44 to 92 MPa. This observation is consistent with stress estimates for other crustal shear zones formed at or near the brittle-ductile transition (e.g., Behr and Platt, 2011, 2014, and references therein).

Based on the orientation of structural fabrics mapped in the field (i.e., foliation and lineation), the BMZ represents a dominantly NE-SW–striking, NW-dipping shear zone (Figs. 1, 11E–11F). Optical and electron microscopic analyses of asymmetric porphyroclasts and other kinematic indicators of 24 samples within the BMZ are consistent and preserve dominantly sinistral, reverse (top-to-the-south) kinematics associated with amphibolite-facies ductile deformation in the BMZ. This substantiates previous studies of the BMZ, which have identified NE-SW–striking, NW-dipping mylonitic foliation and microstructural evidence of oblique sinistral shearing and reverse northwest-over-southeast motion (Kohut and Hepburn, 2004). Microstructural evidence for overprinting brittle deformation is abundant both in the field and in thin section, indicating later top-to-the-northwest, normal motion associated with the formation of the BBFZ (e.g., Fig. 4E).

CVA analyses allow for a major advantage in this study—the results confirm that the orientation of the X-Z plane is approximately parallel to the vorticity-normal surface, supporting the validity of W_m estimates from RGN analyses. In the specimen reference frame (thin section orientation, parallel to the X-Z section of the finite strain ellipsoid), almost all samples have bulk CVA dispersion axes that plot within the foliation plane, perpendicular to the lineation (i.e., parallel to the y-axis of the finite strain ellipsoid). This relationship is consistent with a monoclinic deformation geometry in the shear zone, and places a unique constraint on the interpretation of W_m estimates.

Modeling of structural fabrics in shear zones (Tikoff and Fossen, 1993, 1995; Fossen and Tikoff, 1993, 1998) have established the numerical and geometric relationships between kinematic vorticity and structural fabrics that develop due to changes in the orientation of the finite strain ellipsoid, and changes in the convergence angle of blocks on either side of the deforming zone. The “standard” model for fabric development in shear zones assumes that the lineation is parallel to the x-axis of the finite strain ellipsoid and the pole to foliation is parallel to the z-axis. Assuming perfect simple shear, the vorticity axis lies within the plane of the foliation, perpendicular to the lineation and parallel the y-axis of the finite strain ellipsoid. However, as previously mentioned, results of strain models highlight that foliation and lineation orientation are not always diagnostic for determining vorticity as they do not necessarily correspond to the plane perpendicular to the vorticity axis. In fact, during 3-D deformation, the vorticity axis can be oriented perpendicular (e.g., simple shear), parallel (e.g., pure shear–dominated transpression; e.g., Tikoff and Greene, 1997), or oblique (e.g., triclinic shear zones; e.g., Jiang and Williams, 1998) to the lineation direction.

The results of the CVA analysis, when interpreted in the context of the relatively consistent field observations of the BMZ orientation (NW-dipping, NE-SW–striking foliation plane), suggest that the BMZ corresponds to a shear zone system with an overall monoclinic geometry (Fig. 11). Indeed, both the phase-specific (Fig. 9; Fig. S4 [^{footnote 1}]) and bulk CVA patterns (Fig. 8) correlate to the orientation of the vorticity axis expected for simple shear and wrench-dominated transpression involving general shear (Michels et al., 2015).

This observation is further supported by constraints on the mean kinematic vorticity number (W_m) estimated using RGN analysis (Fig. 6), which ranges from 0.25 to 0.55, corresponding to general shear. The relationship between foliation and lineation and the vorticity axis combined with the mean kinematic vorticity number therefore provide a critical new constraint on the deformation geometry of the BMZ (Fig. 11). In aggregate, these data show that crustal strain localization along the Avalon-Nashoba terrane boundary involved the combined effects of pure and simple shear in a monoclinic, simple shear and/or wrench-dominated sinistral, transpressional shear zone formed at or near the brittle-ductile transition (∼450–600 °C, ∼44–92 MPa).

In order to exploit the determination of ISAs and flow apophyses as calculated from the results of RGN analysis, we use the bulk geographic vorticity axis that we calculated from CVA analyses in the geographic reference frame to present a feasible range of deformation geometries for the BMZ. Specifically, we use the angles α and θ that we calculated for minimum and maximum W_k values to define geographic vectors for ISAs and flow apophyses relative to the dominant foliation, lineation, and vorticity axis. Of particular interest is the geographic orientation of the flow apophysis nearest the orientation of ISA₃, as the trend of this apophysis corresponds to an estimate of the geographic tectonic convergence direction associated with deformation in the BMZ. In this context, we constrain the relative angle of tectonic convergence between the Nashoba and Avalon terranes to ∼56°–75° (Fig. 11A) with the geographic orientation of convergence trending 141.6°–159.2° and plunging 3.2°–10.4°, respectively (Fig. 11H).

The way structural geologists typically conceptualize mineral-shape fabrics is at odds with the manner in which the fabrics are often analyzed. The main discrepancy is that such a fabric can (but not always) exhibit both planar and linear alignments of grain shapes, yet these features are typically recorded and treated separately during analyses. As such, any presumed genetic link between these structures is effectively severed by the convention of measuring and analyzing them separately.

Instead of separating planar and linear features, we can use the concept of a 3-D orientation to simultaneously represent the direction of lines within planes and the facing direction of the planes in which the lines exist. The analytical tools necessary to work with these “combined” representations of rock fabrics are described by the field of orientation statistics. Importantly, the term “orientation” in this context describes a complete 3-D orientation of a body, and can be applied to numerous geological structures such as a foliation-lineation pair, a crystal lattice, a fold, an ellipsoid, etc. For a useful review of methods for analyzing orientation data in structural geological applications, see Davis and Titus (2017).

This study benefits from orientational representations of field-scale mineral-shape fabrics and crystal lattices. Specifically, we define orientations that comprise foliation-lineation pairs (lines within planes) relative to a geographic reference coordinate system. The main advantages of this approach in our study are that: (1) results from the specimen reference frame (i.e., CVA results, ISAs, CPO patterns) can be easily rotated to a geographic reference frame and analyzed or compared within the tectonic framework; and (2) it allows us to calculate an average fabric orientation in a holistic manner—that is, we can calculate a mean of the foliation-lineation pairs, rather than calculate the mean directions of lineations and foliation poles separately (e.g., Figs. 11E–11F). The second point is crucial for constraining the dominant fabric orientation, which allows us to estimate the geographic convergence direction associated with the deformation in the BMZ based on results from RGN analyses. To conduct this type of analysis, we use custom code that leverages functions in the MTEX toolbox for MATLAB. We provide open access to our code and example scripts in an online public repository at https://github.com/zmichels/Fabrica.

This research also highlights the use of CVA analysis as an objective check on the validity of surfaces selected for RGN vorticity analysis. The methods of RGN analysis require that the analyzed section surface be oriented normal to the vorticity vector of the associated non-coaxial deformation geometry. Furthermore, RGN methods are also only valid for monoclinic deformation geometries (i.e., simple shear and transpression-transtension). The CVA analyses in this study were conducted on EBSD maps from sections that were cut parallel to lineation and perpendicular to foliation (i.e., the “standard” fabric reference frame as determined from samples collected in the field). As mentioned above, our results indicate that the sections are well aligned with the vorticity-normal section. Although this result is not surprising (i.e., it matches the common presumption about the significance of foliation and lineation in many mylonitic tectonites), it allows us to objectively rule out the possibility that the fabrics may have resulted from pure shear–dominated transpression (for which the lineation direction can be parallel to the vorticity vector) or triclinic (for which the relationship between lineation and vorticity may be oblique) deformation geometries.

For the sake of shear zone analyses, CVA provides a unique advantage because it is scalable and provides a direct link between the microkinematic behavior of crystals during deformation and the larger-scale deformation geometry of tectonic shear zones. The main premise of CVA methods is that the bulk rotational axis of a deformation history can be inferred from the aggregate of many smaller-scale rotational axes (e.g., Michels et al., 2015). This premise suggests that at the grain-scale, any individual grain may not yield a representative lattice-rotation axis; however, in bulk, the polycrystalline aggregate must conform to the larger-scale kinematics. This premise may be carried to field-scale mineral fabrics and the deformation geometries of tectonic-scale shear zones. In this manner, we consider that results of CVA analyses can likely be extended beyond the framework of a single sample and combined from numerous samples in a geographic reference frame to better constrain the rotational axes and deformation geometries associated with plate tectonic–scale deformations, as we show in our analysis (Fig. 11). Furthermore, this premise also suggests that the inverse scaling relationship is ripe for exploration. Specifically, if the preferred crystal lattice rotation axis in any given sample is controlled by the bulk deformation geometry, then perhaps the differences in the orientations of those axes can be used to quantitatively investigate the kinematics of strain partitioning within samples, or at larger spatial scales.

Finally, this study provides the first detailed demonstration of the utility of CVA analyses in polymineralic aggregates. The debut of CVA analysis presented by Michels et al. (2015) focused on quartzites in order to demonstrate the potential for the method in a relatively simplistic rock type. However, one of the main advantages of CVA analysis purported by the authors is that it can be applied to any crystalline mineral phase in a deformed material, permitting a unified approach for analyzing the rotation axes of grain interiors in many rocks. Indeed, the results of our CVA analyses indicate a strong agreement between the orientations of lattice rotation axes from major mineral phases in the BMZ. This observation supports two important inferences: (1) CVA analyses likely can be applied to numerous mineral phases—not only quartz; and (2) CVA analysis may represent a useful new tool for exploring the roles that different mineral phases may have in sample-scale partitioning of strain during deformation of polymineralic rocks. Moreover, to our knowledge, this is the first study in which rigid grain methods (for the determination of mean kinematic vorticity) have been used in conjunction with orientation-dispersion methods (specifically CVA analysis) to uniquely constrain not only the deformation geometry of an enigmatic shear zone, but also the plate-scale convergence of adjacent tectonic blocks.

Field and microstructural analyses of the BMZ provide new constraints on crustal strain localization along the Boston Avalon–Nashoba terrane boundary in eastern Massachusetts. The quartz-rich mylonites of the BMZ preserve microstructural, kinematic, and vorticity patterns that are interpreted to reveal previously unknown aspects of the deformation history and shear zone kinematics. The BMZ rocks analyzed comprise variably mylonitized quartzofeldspathic gneisses, quartzites, and granites. The presence of hornblende-bearing samples implies metamorphism up to amphibolite-facies conditions. Quartz CPO patterns and misorientation axes suggest the dominant activation of prism<a> slip, with minor evidence for contributions of basal<a>, rhomb<a>, and prism[c] slip, suggesting that the BMZ mylonites record intermediate temperatures of deformation ranging from ∼450 to 600 °C consistent with observed metamorphic assemblages. Thin section–scale observations of quartz further support this interpretation, as they reveal diagnostic microstructures indicative of subgrain rotation recrystallization and grain boundary migration dynamic recrystallization, consistent with temperatures inferred for deformation dominated by dislocation creep. EBSD-generated grain size data applied to the recrystallized grain size piezometer for quartz yield differential stress estimates for the BMZ mylonites ranging from ∼44 to 92 MPa, further suggesting that deformation within the BMZ occurred under relatively high stress conditions and moderate temperatures at or near the brittle-ductile transition.

Field observations and microstructural relationships indicate dominant sinistral top-to-the-south shear-sense kinematics within the BMZ. The RGN-estimated mean kinematic vorticity number ranges from 0.25 to 0.55, indicating general shear. Using the angles α and θ calculated for the minimum and maximum W_k values, we constrain the angle of tectonic convergence between the Nashoba and Avalon terranes in eastern Massachusetts to ∼56°–75° (Fig. 11A) with the geographic orientation of convergence trending ∼142°–160° and plunging ∼3°–10° (Fig. 11H). Patterns of crystallographic dispersion calculated using CVA analyses indicate that the axis of kinematic vorticity lies within the plane of foliation, perpendicular to lineation—a pattern consistent with monoclinic deformation geometries involving simple shear and/or wrench-dominated transpression. Combined field, microstructural, and vorticity analyses are therefore interpreted to suggest that crustal strain localization along the Avalon-Nashoba boundary, as recorded in the BMZ mylonites, involved the combined effects of pure and simple shear in a predominantly sinistral reverse top-to-the-south simple shear and/or wrench-dominated transpressional shear zone. If previous studies indicating that the BMZ represents the accretionary boundary of the Avalon terrane with the Nashoba terrane located along the Laurentian margin are correct, this indicates that the Avalon terrane moved into and below the Nashoba terrane at a shallow angle from the southeast. Further, because the Acadian orogeny in the northern Appalachians is thought to have resulted from collision between the Avalon terrane and the Laurentian margin (e.g., van Staal et al., 2009; Fyffe et al., 2011; Pollock et al., 2012), this study defines, for first time, the motion of the driving force for the Acadian orogeny, at least in southeastern New England.

This study demonstrates that the new CVA methods can be used in conjunction with more traditional methods of structural and vorticity analysis to reveal valuable insights into deformation geometries, without assumptions about the relationship between kinematic vorticity and fabric-forming elements. Further, the successful application of the CVA method to determine both bulk and major-phase dispersion axes in a polyphase shear zone is promising for future studies of large-scale crustal shear zones, particularly given that many aspects of microstructural interpretation are limited to inferences of a specific phase, such as quartz.

We would like to thank Elena A. Miranda and Jeffrey Rahl for detailed and constructive reviews that greatly improved this manuscript. We thank Shanaka de Silva and Allen McGrew for their editorial handling of this paper as part of the special volume honoring Arthur Snoke. Seth C. Kruckenberg would like to personally thank Art Snoke for introducing him to structural geology, metamorphic core complexes, and integrated field-based studies, thereby shaping the future career of yet another student who was lucky enough to pass through his doors at the University of Wyoming and fall under his mentorship. The authors would like to further acknowledge his enumerable contributions to the field of structural geology and to the tectonics community through his education and research efforts, which continue to serve as a great source of inspiration.