Fabric transitions in quartz via viscoplastic self-consistent modeling part I: Axial compression and simple shear under constant strain

Quartz is a common crustal mineral that deforms plastically in a wide range of temperatures and pressures, leading to the development of different types of crystallographic preferred orientation (CPO) patterns. In this contribution we present the results of an extensive modeling of quartz fabric transitions via a viscoplastic self-consistent (VPSC) approach. For that, we have performed systematic simulations using different sets of relative critical resolved shear stress of the main quartz slip systems. We have performed these simulations in axial compression and simple shear regimes under constant Von Mises equivalent strain of 100% (γ = 1.73), assuming that the aggregates deformed exclusively by dislocation glide. Some of the predicted CPOs patterns are similar to those observed in naturally and experimentally deformed quartz. Nevertheless, some classical CPO patterns usually interpreted as result from dislocation glide (e.g. Y-maxima due to prism slip) are clearly not developed in the simulated conditions. In addition we reported new potential preferred orientation patterns that might happen in high temperature conditions, both in axial compression and simple shear. We have demonstrated that CPOs generated under axial compression are usually stronger that those predicted under simple shear, due to the continuous rotation observed in the later simulations. The fabric strength depends essentially on the dominant active slip system, and normally the stronger CPOs result from dominant basal slip in , followed by rhomb and prism [c] slip, whereas prism slip does not produce strong fabrics. The opening angle of quartz [0001] fabric used as a proxy of temperature seems to be reliable for deformation temperatures of ~ 400 °C, when the main slip systems have similar behaviors.


Introduction
Quartz is a common crustal mineral that deforms plastically in a wide range of temperatures and pressures, from very low metamorphic grades just above diagenesis, to sub-solidus temperatures. Such a broad range of metamorphic conditions, together with different imposed deformation geometries that aggregates may experience naturally and experimentally, leads to the development of different microstructures and crystallographic preferred orientation (CPO) distributions (e.g. Tullis, 1992, 1994;Law and Johnson, 2010;Schmid and Casey, 1986;Tullis et al., 1973). The development of CPO in quartz might have multiple origins, from mechanical rotation of crystals in very low-grade metamorphic rocks (e.g. Stallar and Shelley, 1995) to dissolution-precipitation (e.g. Hippertt, 1994;Hippertt and Egydio-Silva, 1996). In most cases however, the CPO of quartz is interpreted in terms of dislocation activity and the activation of different crystal slip systems.
Numerical modeling of CPO development has been conducted in the last four decades by different approaches. These studies were initially conducted via kinematic (Etchecopar, 1977) or relaxed-constrain Taylor (Lister and Paterson, 1979;Lister et al., 1978) approaches. However, the 1990s saw an evolution in the development of different modeling methods that simulated the development of preferred orientations in different types of materials. Such evolution was due mainly to improved understanding of the physical processes underpinning the development of CPO in crystalline materials, combined with a significant improvement in computational processes that allowed the development of complex codes capable of incorporating a number of different parameters that influence CPO development. These approaches include equilibrium-based models (e.g. Chastel et al., 1993), isotropic models (Takeshita et al., 1990), n-site models (Wenk et al., 1991), and viscoplastic self-consistent models (Lebensohn and Tomé, 1993). All of these methods have in common the assumption that deformation is accommodated by dislocation glide on a limited number of crystallographic planes in specific crystallographic directions. Furthermore, they all consider also that polycrystalline behavior can be averaged by considering the response of individual microscopic (grain) responses. Thus, the lower bound state assumes that the stress in the aggregate is homogeneous, whereas in the upper bound state it is the strain in an aggregate that is homogeneous (Sachs, 1928;Taylor, 1938).
An important aspect considered in only a few of the numerical models of CPO development is dynamic recrystallization (e.g. Etchecopar, 1977;Lebensohn et al., 1998;Wenk et al., 1997). Although the occurrence of dynamic recrystallization in parallel with crystal-plastic mechanisms is well-know from experiments and nature, the mechanisms by which the recovery occurs (i.e. subgrain rotation recrystallization, grain boundary migration, etc.) and their impact on the CPO development remain poorly constrained. More recently, integration of viscoplastic self-consistent methods with bi-dimensional microstructural network configurations defined by computer codes such as ELLE has provided new opportunities to better study the effects of dynamic recovery mechanisms concurrent with crystal plasticity (e.g. Griera et al., 2011).
Other processes that are normally not considered in these approaches but are of great importance for the development of CPOs is the effect of dislocation climb and twinning (in quartz, more specifically, Dauphiné twinning). Dislocation climb modeling has to consider the effect of diffusion of atoms in the lattice, as well as the shear components of applied stress in addition to the resolved shear stress described by the diagonal components of the applied stress matrix (e.g. Lebensohn et al., 2012), leading to complex computations. Dauphiné twinning is an important mechanism because it is potentially a grain weakening mechanism due to its effect on stiffness/compliance of quartz single crystals, as it switches − ve & + ve orientations, such as the b aN-axes and the rhombs (e.g. Lloyd, 2004;Menegon et al., 2010;Pehl and Wenk, 2005;Wenk et al., 2006Wenk et al., , 2007Wenk et al., , 2009a. In this contribution, the viscoplastic self-consistent approach (Lebensohn and Tomé, 1993) is used to investigate crystallographic preferred orientation transitions in quartz via a systematic modeling of the influence of multiple slip systems in quartz preferred orientations. For that, we assumed systematic changes in the relative critical resolved shear stress (CRSS) for the main slip systems in quartz during axial compression and simple shear. We assume that the main control on the activation of different slip systems in quartz is temperature, and the quartz slip systems are activated in the following order with increasing temperature: baN(c) → baN{r/z} → baN{m} → [c]{m}. This assumption is justified by the classical works of Hobbs (1968), Baeta and Ashbee (1970) and Tullis et al. (1973) that indicate that prismatic slip in quartz is dominant over basal slip in higher temperatures (and vice-versa) and that the increasing temperature leads to a switch between basal baN to rhomb b a N slip. According to Lister and Dornsiepen (1982) and Mainprice et al. (1986), at 'normal' geological strain rates, the change from basal to prismatic glide occurs at 600-700°C. The increase in temperature nevertheless does not imply that "low" temperature slip systems become inactive, but that "low" and "high" temperature slip systems may occur together (e.g. Egydio-Silva et al., 2002;Kruhl, 1996;Menegon et al., 2011;Tommasi et al., 1994). It is important to note that this classical view has been questioned in the last few years (see review by Toy et al. (2008)) and some authors clearly demonstrate that CPO in quartz may also be strain-dependent (e.g. Heilbronner and Tullis, 2006), and this will be discussed in the second part of this paper where we will present the detailed fabric transitions in quartz using the same systematic approach of slip system combinations under high strain simple shear. Finally, in the third part of the paper we will discuss the characterization of orientation distribution functions, orientation matrices and different parameters of CPO quantification resulting from the simulations presented in parts 1 and 2.

Viscoplastic self-consistent modeling (VPSC)
In this paper, we explore the possibilities of the viscoplastic selfconsistent approach (Lebensohn and Tomé, 1993) using the software VPSC7 compiled for Mac OS X. Through this approach, the overall response of the aggregate is averaged from known properties of the constituent grains, and the effective medium assumption is considered to calculate the interactions between each grain with the surrounding ones. In the VPSC approach, strain compatibility and stress equilibrium are averaged at the aggregate scale, and the microscopic stress and strain rate may differ significantly from the macroscopic values. The plastic anisotropy is also considered in this approach. Although this factor has a minor influence in quartz because of the relatively large number of slip systems, it is significant in materials with few slip systems.
The VPSC approach considers that deformation takes place when one or more slip systems become active. The shear rate induced in a slip system (s) by a given local deviatoric stress (s ij ) is described by a viscoplastic law of form (Eq. (1)) where γ 0 is the reference strain rate, n s is the stress exponent (i.e. for diffusion creep, n = 1, for dislocation creep in silicates, 3 ≤ n ≤ 5), τ r s and τ 0 s are the CRSS and the activation stress respectively for a given crystal slip system and r s is the Schmid tensor, which defines the orientation of the slip system in relation to the stress axes. The sum of the shear rates over all slip systems gives the microscopic strain rate (ε ) per grain. To determine the microscopic state (s, ε ) of each Fig. 1. Reference frames used in this paper, for axial compression and simple shear simulations. In the axial compression pole figures, the crystallographic orientation data is plotted in relation to the axial compression direction (Z) and the flattening plane (E-W vertical plane). In simple shear, the data is plotted in relation to the strain ellipsoid main features (foliation/lineation-E-W vertical plane/E-W horizontal pole) and in relation to the imposed shear plane/shear direction (NE-SW vertical plane/NE-SW horizontal pole). In the inverse pole figures the axial compressional direction is plotted in relation to the main crystallographic forms of quartz.
grain in an aggregate, together with the volume average (Σ; E -Eq. (2)) that determines the response of the aggregate, the 'one-site' simplification (Lebensohn and Tomé, 1993;Molinari et al., 1987) is used (Eq. (2)) This approach does not take into account the interaction between individual grains. Rather, the interaction between each grain and its surrounding grains is substituted progressively by the interaction between an inclusion with an equivalent crystallographic orientation and an equivalent infinite homogeneous medium (HEM). The behavior of this artificial medium is the weighted average of the behavior of the aggregate. To solve the problem of a deformed ellipsoidal domain in an infinite homogeneous body, the approach of Eshelby (1957) is used, which assumes ellipsoidal grains and a tangential HEM, where e M is the interaction tensor that depends on the rheological properties of the aggregate and the shapes of grains and α is a constant (range zero to infinity) used to parameterize the interaction between the homogeneous medium and the grains, which imposes essentially kinematically rigid conditions on to an aggregate. A zero value for α corresponds to the homogeneous strain in the aggregate and is equivalent to the Taylor-Bishop-Hill modeling approach (e.g. Lister et al., 1978;Taylor, 1938) in which all grains deform at the same rate and have the same shape at each imposed strain step. Such behavior leads to stress incompatibilities at grain boundaries, which have to be corrected by elastic strains (Wenk et al., 1989). A value of infinity for α means that all crystals experience the same state of stress (the so-called Sachs model; Sachs, 1928) but this does lead to microstructural incompatibilities on the grain scale as compatibility is obtained only on the aggregate scale.
The rotation (ω ij ) of the crystals in relation to the external reference frame is described by In this equation, Ω ij is the anti-symmetric component of the imposed macroscopic velocity gradient, which defines the type of imposed deformation we want to simulate. The summation term defines the plastic spin, where b is the Burgers vector, n is the normal direction of the slip plane s and e ω is the reorientation rate of the inclusion. This last term depends on the difference between the strain rates of the individual grain and the whole aggregate and increases with deformation. The velocity gradients (L ij ) that control the type of imposed deformation in the simulations are constant and are given by the following matrices, for axial compression and simple shear, respectively: For the axial compression simulations, the pole figure external reference frame is defined as a horizontal N-S axial compression direction and the vertical E-W flattening plane (Fig. 1). For the simple shear experiments, the finite strain is given by the vertical E-W foliation and the horizontal E-W lineation, whereas the shear plane is given by the NE-SW vertical plane and the shear direction given by the horizontal NE-SW pole (Fig. 1).

Modeling conditions
VPSC modeling was performed assuming 1000 initially spherical grains with random crystallographic orientation. After the initial setup, the parameters that can be changed and investigated in different models are: (i) the (combinations of) crystal slip systems that operate; (ii) the relative CRSS of each slip system; (iii) the stress exponent (n s ); and (iv) the interaction between the grains and the HEM (i.e. α). The models reported here considered constant values of n s = 3, and α = 1. In contrast, the relative values of CRSS were varied systematically from 1 to 10, where the former defines a slip system that is easily activated and the later a slip system that is hard to activate. All the models were carried out up to a Von Mises equivalent strain (ε eq ) of 100%, meaning a shear strain of γ = 1.73. The presentation of 'static' CPOs related to this specific strain value is justified by the fact that the crystallographic preferred orientations in VPSC models evolve very fast, essentially because only dislocation glide is considered. Normally at equivalent strains of 100% the CPO patterns are already very clear, and essentially the fabric just gets stronger. In simple shear nevertheless, the quartz fabric still evolves for γ N 1.73 (e.g. Keller and Stipp, 2011) and the static CPOs at this shear strain may not properly represent high strain CPOs. For that reason, we performed high-strain simple shear simulations using the same approach of systematic variation of the relative CRSS of different slip systems in quartz and present the results in the form of CPO evolution with strain in part 2 of this paper.  The ε eq is defined by where the Von Mises criterion D eq is described by and the transformation from equivalent strain to shear strain is given by The crystal slip systems selected for use in the simulations have all been defined in previous single crystal experiments and via transmission electron microscopy analyses (e.g. Baeta and Ashbee, 1969;Griggs and Blacic, 1965;Kronenberg and Tullis, 1984;Linker et al., 1984;Mainprice and Jaoul, 2009;Paterson, 1984, 2005;Rutter and Brodie, 2004). These slip systems are listed in Table 1. Although systems such as dipyramids with slip directions oblique to b cN and baN directions are relatively hard to activate and therefore are virtually inactive in the VPSC simulations, they have been listed in Table 1. However, acute rhomb {±π} slip in the baN direction, which has been recognized in some quartz ultra-mylonites deformed in simple shear (e.g. Law et al., 1990;Lloyd, 2004;Lloyd et al., 1992;Morales et al., 2011), has not been considered in the models due to their angular proximity with the standard rhomb slip systems {r/z} and {z} in b aN direction. Orientation distribution functions and pole figures were calculated using the software PFvpsc (Mainprice, 1990), whereas inverse pole figures were plotted using the MTEX toolbox for Matlab (Hielscher and Schäben, 2008).
Prior to the simulations with combined slip systems, we performed a series of models to illustrate the development of preferred orientations when the slip systems act alone (Fig. 2). The relative CRSS values for the simulations with multiple slip systems are represented as simple ratios. For example (see the x-axis in Fig. 3), if baN(c)/baN{r/z} = 1/10, the CRSS for b aN{r/z} slip is 10 times stronger than the CRSS for baN(c) slip, such that the latter system dominates the CPO. In contrast, if baN(c)/b aN{r/z} = 1/1, the values of CRSS are identical for both systems. Obviously, if b aN(c)/baN{r/z} = 10/1, then the CRSS for basal ba N slip is ten times higher than that for rhomb b aN slip. Similarly, if only a single slip system is plotted (e.g. {m}[c] slip on the y-axis in Fig. 3), then a value of 10 means that the system is virtually inactive, whereas a value of 1 means that it has effectively the same relative CRSS as the other systems under consideration.

Transitions between baN(c) and bbaN{r/z} at constant [c]{m}
When baN(c) is the dominant slip system (i.e. b aN(c)/baN{r/z} = 1/10 and [c]{m} = 10 in Fig. 3), the strongest concentration of [0001] axes is subparallel to the axial compression (Z) axis in the pole figure.
Keeping the relative CRSS of the [c]{m} slip system constant while decreasing that of the baN{r/z} system induces the opening and widening of the conical girdle around Z (direction of maximum compression) (Fig. 3). For example, the tight conical c-axis girdle about Z opens progressively from~20°when the relative CRSS for b aN(c) slip and b aN{r/z} slip are equal, to between 45°and 90°depending on the relative CRSS of [c] {m}. The stronger effect of [c]{m} is better observed when basal b aN is equally weak. The extreme condition is the development of a c-axes single girdle parallel to the flattening plane (XY).
In the equivalent inverse pole figure (Fig. 4), the axial compression direction is not exactly parallel to (0001) when b aN(c) is the dominant slip system. The decrease of relative CRSS for b aN{r/z} slip induces first the weakening of the CPO around (0001), where then the direction of maximum compression becomes parallel to (10-12) and (10-11) and progressively outwards across most orientations except for adjacent to the {11-20} a-axis. A progressive decrease of CRSS on [c]{m} induces a similar effect of spreading the orientations from of the axial compression direction from (10-12) and (10-11) to the first order (10-10) and (01-10) prismatic planes. Two other trends are observed in the inverse pole figures (Fig. 4). Firstly, the direction of maximum compression is never parallel to either the pole of the (11-20) plane or the c-axis, maintaining a low angle with the latter at least when the equivalent strain is relatively low (i.e.~1). Secondly, the concentration of Z-axes parallel to the rhomb form is weaker compared with the concentrations close to the basal and prismatic planes.
The activity of baN(c), baN{r/z}, baN{m} and [c]{m} present nonlinear variations following different setups of relative CRSS used in our modeling. In the transition basal → rhomb ba N slip under axial compression, about 80% of the strain is accommodated by basal slip in ba N, whereas rhomb slip in b a N accommodates the other 20%. (Fig. 5a), assuming that prismatic slip in [c] is virtually inactive. When the CRSS of basal and rhomb slips in baN is equally weaker the activity of basal baN increases from 50 to 62% of strain accommodation, whereas rhomb b aN decreases from 50 to 37% (Fig. 5c). When rhomb b aN slip is dominant, its activity decreases from 80% to about 70% with increasing strain, followed by 10% of activity of both basal ba N and prism [c] slip (Fig. 5e). The activity of these two slip systems is the same when they are equally weak (relative CRSS for b aN(c) and [c]{m} = 1) and decrease slightly from~40 to~35%, followed by an increase of rhomb ba N activity from~17 to~30% of accommodated strain with progressive deformation (Fig. 5b). When b aN(c), b aN{r/z} and [c]{m} are equally weak, the activity of rhomb slip is dominant and constant (~45%), and the basal b a N and prism [c] have constant activities of~25% with increasing strain (Fig. 5d). When the CRSS of basal slip is high, the activity of rhomb slip is dominant (~50%) and increases smoothly with increasing strain, followed by a decrease on the [c]{m} slip activity up to 40% (Fig. 5f).

Transitions between baN{r/z} and b aN{m} with constant [c]{m}
A dominance of rhomb b aN over prismatic b aN slip results in conical c-axis girdle with relatively large opening angles in relation to the axis of compression (Z) (Fig. 6). While a progressive decrease in the relative CRSS of b aN{m} slip does not modify significantly the opening angle of the conical girdles up to a ratio of baN{r/z}/baN{m} of 1/1, the CPO changes considerably when baN{m} is the dominant slip system (Fig. 6), leading to development of c-axis girdle parallel to the flattening plane. The decrease in the relative CRSS for [c]{m} from 10 to 1 always results in the same type of girdle parallel to parallel to the flattening (XY) plane.
The activities related to the b aN{r/z} → baN{m} transition on the other hand are very different. When prismatic slip in [c] and b aN directions are virtually inactive, the deformation is essentially accommodated by rhomb slip in baN (80%), with minor activity of the prismatic slip systems with increasing strain (Fig. 8a). If baN{r/z} and b aN{m} are equally weak, the activity of the former increases from 50 to 70%, while the later decreases from 50% to~30% with increasing strain (Fig. 8c). The deformation is essentially accommodated by prism b aN slip when this is the weak slip system, although prism [c] and rhomb ba N accommodated 20% of the imposed deformation (Fig. 8e). If baN{r/z} and [c]{m} are equally weak, the activity of the former increases from 50 to 60%, while the later decreases from 50 to 40% with increasing strain (Fig. 8b). When the CRSS of baN{r/z}, baN{m} and [c]{m} is equally weak, the activities of rhomb b aN and prismatic [c] slip are about the same and together accommodate~80% of the deformation, while prismatic slip accommodated 20% (Fig. 8d). The deformation is partitioned between prism baN and prism [c] with increasing strain when both have lower CRSS than rhomb baN slip, but the activity of the former increases with imposed strain in relation to the later (Fig. 8f).

Transitions between baN(c) and baN{m}, with constant b aN{r/z}
The transition from basal b aN to prismatic b aN ( Fig. 9) is similar to that for baN(c) → baN{r/z} (Section 3.1.1) but the increasing of opening angle of the conical girdles around Z is less pronounced (see Fig. 3). The  progressive decrease on the relative CRSS ratio between b aN(c)/b aN{m} from 10/1 to 1/1 tends to weaken the conical girdles only slightly, but changes considerably when the ratio is 1/5 or 1/10 (Fig. 9). In this case a c-axis girdle parallel to the flattening (XY) plane develops, with secondary concentration of c-axes remaining about Z.
In the inverse pole figures, the transition from b aN(c) to b aN{m} slip is marked by a strong concentration of compression axis parallel to the c-axis, which then moves progressively toward the rhomb {10-12} and {10-11} (Fig. 10) and finally to prismatic planes (01-10) and (10-10). The decrease of the relative CRSS of baN{r/z} displaces slightly the axial compression axis toward the (10-12) rhombs. The concomitant decrease on the relative CRSS ratio of baN(c)/baN{m} from 1/10 to 10/1 weakens the CPO but removes completely the concentrations about (01-10) and (10-10).

Transitions between baN(c) and baN{r/z} with constant [c]{m}
Irrespective of the relative CRSS ratios between baN(c) and b aN{r/z} slip systems, all c-axis distributions are characterized by girdles normal  10 (a, b), the CRSS for baN{r/z} is 10 times stronger than the CRSS of baN(c) slip. When the relation is 1/1 (c, d), both slip systems are easily activated, whereas 10/1 (e, f) means that basal slip in baN is 10 times stronger than rhomb slip in baN.
to the foliation and oblique to the shear plane when [c]{m} slip is not favored (Fig. 11a). In detail, when baN(c)/baN{r/z} ratio is 1/10, the maximum concentration of [c]-axes is parallel to Z but with the decrease of the relative CRSS for baN{r/z} the maxima moves to about~30°from Z in the YZ plane. The decrease in the relative CRSS for [c]{m} induces the widening of the c-axis girdles, then the development four c-axis maxima located at~90°to each other within the XZ plane and positioned symmetrically at~45°with respect to the foliation (XY) plane (and the lineation, X). Finally, two maxima migrating to a position that bisects the foliation and shear planes are developed when both baN{r/z} and [c]{m} have relative CRSS of 1.
When the relative CRSS for [c]{m} is 10, poles to {10-10} planes define an (incomplete) girdle parallel to the foliation and an obvious maximum parallel to X (Fig. 11b). The orientation of the shear plane appears to exert no influence on these pole figures. The changing ratios between the CRSS of b aN(c)/b aN{r/z} from 10 to 1 to 1/1 lead to the CPO weakening and the development of diffuse conical girdles about X. When the relative CRSS of [c]{m} is 1 and the CRSS ratio for baN(c)/baN{r/z} 1/ 10, four 'orthogonal' maxima develop in the XZ plane, but as the ratio changes to 10/1, two of these maxima become dominant and eventually define 'cross-girdles' distributions (Fig. 11b). The poles to rhomb planes {10-11} are distributed along broad, weak girdles  symmetrical about X with two maxima sub-parallel to the shear direction when baN(c)/b aN{r/z} is 1/10 (Fig. 11c). When this ratio changes to 1/1 the girdles weaken and the 'opening angle' about X widens until they are effectively 'straight'. As the relative CRSS ratio for baN(c)/b aN{r/z} changes from 1/10 to 10/1 and [c]{m} decreases from 10 to 1, the weak girdles of {10-11} first weaken further and then are replaced by open (i.e. 'straight') girdles about X, possibly with an absolute maximum sub-parallel to the shear direction (Fig. 11c).
In simple shear, the activity of basal slip in baN decreases progressively from more than 80% to~60% with increasing strain, followed by an increase of rhomb slip in b aN from less than 20% to about 35%, if prismatic [c] is virtually inactive (Fig. 12a). The deformation is equally partitioned if basal and rhomb b a N slip systems are dominant (~50% each) (Fig. 12c), whereas rhomb ba N accommodate about 80% of the strain if it is the dominant slip, with smaller activities of basal baN and prism [c] (Fig. 12e). If basal ba N and prism [c] are equally weak and rhomb b a N is hard, the deformation is equally partitioned between basal and rhomb b aN slip systems, with the activity decreasing slightly with increasing strain, followed by an increase of up to 30% activity of prism [c] slip (Fig. 12b). When all the slip systems are equally weak, the deformation is mostly accommodated by prismatic slip in [c] and  10 (a, b), the CRSS for baN{m} is 10 times stronger than the CRSS of baN{r/z} slip. When the relation is 1/1 (c, d), both slip systems are easily activated, whereas 10/1 (e, f) means that rhomb slip in baN is 10 times stronger than prism slip in baN.
its activity increases up to 50% under high strains, followed by a small drop of baN(c) and b aN{r/z} activities (Fig. 12d). Almost all the strain is partitioned between rhomb b aN and prism [c] when basal slip is inactive, with slightly increasing activity of baN{r/z} in relation to [c]{m} (Fig. 12f).

Transitions between baN{r/z} and b aN{m} with constant [c]{m}
The transition baN{r/z} → baN{m} is characterized by CPO patterns varying from single girdles sub-normal to the foliation (when the relative CRSS of baN{r/z} is 1) to broad single girdles parallel to foliation, when the relative CRSS of baN{m} is 1 (Fig. 13a). The transition is more pronounced when the relative CRSS between b aN{r/z}/b aN{m} is larger than 1/1, and passes through a transition stage in which a broad maximum in c-axes occurs about the intermediate (Y) axis. The orientation of the poles of prismatic planes is weak and maximum position varies from parallel to the shear direction to parallel to the foliation with the increasing relative CRSS ratio between rhomb and prismatic ba N slip system (Fig. 13b). The orientation of rhomb planes is also weak   between baN(c) and baN{m}, and different CRSS for baN{r/z}, for an equivalent strain of 1. The scale on the x and y-axes is the same as that in Fig. 9. Distribution based on the orientation of 1000 grains, contour lines are given in the figure.
( Fig. 13c) but is characterized by symmetrical girdles similar to those from Fig. 11c.
The activity of rhomb slip in ba N is dominant and accommodates around 80% of the deformation during the transition of b aN{r/z} → baN{m} with prism [c] and b aN virtually inactive, with minor (but important) activities of prism b a N and [c] slip systems (Fig. 14a). The strain accommodated baN{r/z} and b aN{m} when they are both weak, but the former becomes more active than the later with increasing strain, whereas prism ba N accommodates again around 80% of the strain when this is the dominant slip system (Fig. 14c and e). When the CRSS of prismatic slip in [c] is assumed weak, the deformation is equally partitioned between rhomb b a N and prismatic [c] slip, but under high strain conditions the activity of the former increases, while the later decreases slightly (Fig. 14b). The strain is partitioned between baN{r/z}, b aN{m} and [c]{m} when they are all weak, with a weak predominance of prism slip in [c] and [a] (~40% each) in comparison to rhomb b aN (Fig. 14d). If rhomb b aN slip is hard, then the deformation is mostly accommodated by prism [c], that increases in activity from 50 to 60%, following a similar drop on the activity of baN{m} (Fig. 14f).

Transitions between baN(c) and baN{m} with constant a baN{r/z}
The change on the CRSS ratio between b aN(c)/b aN{m} from 1/10 to 10/1 has little effect on the shape of the c-axis distributions with a constant b aN{r/z} CRSS of 10, exhibiting maxima normal to the foliation connected by weak foliation normal girdles (Fig. 15a). For low relative CRSS for both basal and rhomb ba N slip, the foliation normal girdle is sharper than if prismatic slip ba N dominates. A weak maximum of c-axes parallel to Y is developed when prismatic and rhomb b a N slip have lower CRSS in relation to basal ba N slip (Fig. 15a). The poles of {10-10} planes are distributed either as girdles parallel to the foliation if basal b a N and/or rhomb ba N are the dominant slip systems, or as three symmetrical maxima, if prismatic b aN dominates (Fig. 15b). The preferred orientation of the rhomb planes is weak and is characterized by two maxima orientated at~25°with X but that progressively disappear with the transition to prism b aN dominated slip (Fig. 15c).

CPO development under axial compression and simple shear
Detailed analyses of the CPO evolution during our simulations clearly demonstrate that the preferred orientations are much stronger when developed under axial compression rather than simple shear. The main reason for that is that once the crystals acquire a stable orientation, for a given set of dominant slip systems in relation to the imposed axial compression, the crystals do not rotate anymore, and preferred orientation just becomes more concentrated. On the other hand, simple shear produces strong CPOs just under relatively high strains as we demonstrated in the second part of this paper, where the crystals are continuously rotated and the CPO is continuously regenerated. Similar observations were made by Wenk et al. (2009b), who demonstrate through VPSC modeling that the CPO in halite is continuously regenerated up to shear strains of 35, and also by Keller and Stipp (2011) for quartz. Our observations are also in indirect agreement with fabric studies in calcite deformed under simple shear, where the authors demonstrated that the strength of CPO becomes stronger with the imposed strain even when the aggregates are completely recrystallized (e.g. Barnhoorn et al., 2004).

Transition baN(c) → b aN{r/z} in axial compression
The variations on crystallographic orientations related to the transition from basal baN to rhomb baN slip were first reported in axial compression experiments of flint (a fine-grained, water-rich polycrystalline form of quartz) and quartzites by Green (1967) Green et al. (1970) and Tullis et al. (1973). According to the authors, the modification of the distribution of quartz [0001] axes is directly related to the variations in temperature and strain rates of the experiments. Under relative low temperature or high strain rates, the quartz [0001] axes are distributed along conical girdles with a maximum concentration parallel to the compressional axis. With increasing temperature or low strain rates, the tight small circle girdle becomes broader about the σ 1 (the Z in our pole figures) and the half opening angle of these girdles may vary from 20-45°and tends to increase with increasing temperature. These observations are clearly compatible with out models (Figs. 3, 4). When the relative CRSS of prismatic slip in [c] is high, the increase in the opening angle in relation to the compressional axis is not very evident, The scale on the x-axis is given by the ratio baN(c)/baN{r/z}. When the ratio is 1/10, the CRSS for baN{r/z} is 10 times stronger than the CRSS of baN(c) slip, and then basal slip in baNdominates the development of crystallographic preferred orientation. When the relation is 1/1, both slip systems are easily activated, whereas 10/1 means that basal slip in baN is 10 times stronger than rhomb slip in baN. The same is valid for the y-axis scale showing the relative CRSS for prismatic slip in [c]. Pole figures are plotted on the lower hemisphere, equal area stereonet. Contour of multiples of uniform distribution for 1000 grains, contour interval = 1, or 0.5 for the dashed line. Maximum density is given by the white square in the pole figures, minimum density by the black points. Black line indicates the foliation and x the lineation, whereas the red line indicates the shear plane and SD the shear direction.
becoming more clear when we assumed that basal ba N slip is harder that rhomb b a N slip, following the increasing activity of the rhomb slip (Figs. 3, 4). With the increasing activity of prismatic [c] slip, the quartz [0001] axes in the pole figures distribution do not vary substantially when basal slip is dominant, but has a strong effect on increasing the opening angle, that may reach 90°in relation to σ 1 if rhomb b aN and prismatic [c] slips are equally weak (Fig. 3). Although the pole figures suggest a clear preferred orientation of the quartz basal planes parallel to the flattening plane of our simulations, this observation does not hold in the inverse pole figures, which show that the compressional axes are at least in a very small angle with the pole of (0001) planes in this transition (Fig. 4). In addition, the concentration in our pole and inverse pole figures is much stronger than that observed in the experiments for similar strain conditions. This is probably related to the fact that our numerical models do not take into account any other mechanism rather than dislocation glide. As observed in the experiments of polycrystalline quartz, dynamic recrystallization either by subgrain rotation or grain boundary migration is normally active under relatively high temperatures (e.g. Tullis, 2002, 2006;Hirth and Tullis, 1992). The effect of dynamic recrystallization is quite ambiguous, but experiments on calcite demonstrate that dynamic recrystallization slows down the continuous strengthening of crystallographic  10 (a, b), the CRSS for baN{r/z} is 10 times stronger than the CRSS of baN(c) slip. When the relation is 1/1 (c, d), both slip systems are easily activated, whereas 10/1 (e, f) means that basal slip in baN is 10 times stronger than rhomb slip in baN.
orientation produced in torsion (e.g. Barnhoorn et al., 2004). In quartz, experiments of annealing in pre-deformed polycrystalline aggregates may induce a strengthening of the fabric, possibly by removing high-energy nonfavorable orientated crystals from the aggregates (e.g. Heilbronner and Tullis, 2002). The absence of recovery effects or secondary deformation mechanisms in our simulation also explains why the crystallographic orientation developed in our simulations is very strong in relatively low strains, if one considers the strains applied in experiments.
The predicted distribution of quartz [0001] axes along a continuous girdle parallel to the foliation in the axial compression simulations (Figs. 3, 6 and 9) has never been observed in nature or in experimental samples as far as we are aware of, at least in quartz aggregates deformed by crystal plastic mechanisms. The generation of such a pattern requires a high activity of prismatic [c] and rhomb ba N slips, otherwise the generated patterns are very open small circle girdles around axial compressional axis . The activation of prismatic [c] slip requires high temperatures (above 650°C), relatively high pressures (~200-400 MPa) and hydrous conditions (e.g. Mainprice et al., 1986), similar to those expected in cooling intrusions. In granites deformed in very high temperatures or in some granulites, it is not uncommon to observe chessboard subgrain boundary patterns in quartz, which result from the mutual activation of prismatic [c] and basal baN slip (e.g. Kruhl, 1996;Mainprice et al., 1986;Menegon et al., 2011;Tommasi et al., 1994). The resulting CPO in rocks with this sort of microstructure nevertheless shows a predominance of basal-slip fabric (e.g. Menegon et al., 2011), and the inference of active slip systems purely by pole figure analyses is obviously incorrect. The mutual activation of rhomb b aN and prismatic [c] slip seems to be less common or not studied in detail yet, but may possibly occur under high temperatures and hydrated conditions, above the α → β quartz transition.
The open small circle girdles (N45°) that occur when the CRSS for prismatic [c] and basal b aN slip are equally low (Fig. 3) are also observed in axial compression experiments of flint, a fine-grained, water rich polycrystalline form of quartz (e.g. Green, 1967;Mainprice, 1981;Mainprice and Paterson, 2005) or in some regions of the Saxony granulites (e.g. Lister and Dornsiepen, 1982). In the pole figures, this sort of fabric resembles the crystallographic orientation dominated by rhomb baN slip. In the inverse pole figures however, σ 1 is predominantly parallel to the (10-10) and (01-10), with no concentrations parallel to (0001) and just a weak concentration parallel to the rhomb axes, with the activity of [c]{m} decreasing together with baN(c), followed by an increase on the b aN{r/z} activity (Fig. 5). Although this pattern can be generated by dislocation glide as predicted here, observations on the microstructures of flint led to the conclusion that those patterns are generated by preferential growth of some grains due to strain energy variations related to orientation difference in a deforming environment (Green, 1967). An alternative explanation is where the fabric results from crystal growth of mechanically twinned (Brazil law) crystals, favoring the development of low energy twin grain boundaries with neighbor crystals (Mainprice and Paterson, 2005). On the other hand, the fabrics reported in Lister and Dornsiepen (1982) seem to be compatible with our models. According to the authors, the transition observed in the Saxony granulites results from the basal ba N → prism [c] slip transition (Figs. 3 and 4). Although rhomb baN slip is also considered in these simulations, the relative CRSS for this slip system has to be harder or the same as basal baN slip, otherwise the predicted pattern is a single girdle parallel to the flattening plane in the pole figures.

Transition baN{r/z} → baN{m} in axial compression
The transition from rhomb b a N to prismatic ba N slip under axial compression produces similar pole figures as the transition from basal ba N → rhomb b a N. The main difference is that the later transition is much more evident from the pole figures (Fig. 6). As an example, girdles parallel to the foliation (XY) may be develop even in conditions where prismatic [c] slip is not active, and become relatively strong once this slip system becomes more active (Figs. 6 and 8). The initial transition from dominant rhomb b a N slip, to conditions where rhomb ba N and prismatic b aN are equally weak, is marked by CPO randomization and concentration of compressional axes to intermediary crystallographic axes, such as (11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22) or parallel to acute rhombs (10-12) and (01-12) (Figs. 6 and 7). Although slip systems containing these slip planes are not known for quartz, the concentration of the poles parallel to these unusual axes results from the geometrical solution from the Eqs. (1) to (4), in addition to the statistics applied to calculate the pole figures and inverse pole figures. In fact, the inverse pole figures illustrate that even when rhomb b aN slip is dominant, the concentration of compressional axes parallel to the rhomb forms is rather low, differently to what is observed when prismatic slip baN and/or [c] are dominant, leading to very strong concentrations parallel to the prismatic axes (Fig. 7). This suggests that prismatic slip in baN or [c] geometrically are likely to developed strong CPOs rather than rhomb b aN slip. The scale on the x-axis is given by the ratio baN{r/z}/baN{m} as explained in Fig. 11. Pole figures are plotted on the lower hemisphere, equal area stereonet. Contour lines, shading and reference frame as explained in Fig. 11.

Transition b a N (c) → baN{r/z} in simple shear
The progressive weakening observed in the quartz pole figures (Fig. 11) in the transition from dominant basal rhomb-slip in baN direction possibly means that while basal slip is easily aligned with the imposed simple shear, the same does not happen to the b aN{r/z} slip system. On the other hand, the CPO strength increases again with the increasing activity of prismatic slip in [c], leading to the development of maxima of [0001] axis distributed in a low-angle conjugate girdle with a maxima orientated parallel to the shear direction and at low angle with the when the later is the dominant system. This behavior is the same as that observed in the axial compression simulations, but because the CPO in axial compression is stronger, the effect of increasing activity of b aN{r/z} is not so obvious.

Slip system activity
From these observations, it is clear that basal ba N slip accommodates more strain during the transition b aN(c) → b aN{r/z} during axial compression, assuming that the relative CRSS of b aN{r/z} is harder Fig. 14. Evolution of the slip system activities during the transition baN{r/z} → baN{m} under axial compression with two sets of CRSS for prismatic [c] slip. The CRSS is given in relative values. When the ratio is 1/10 (a, b), the CRSS for baN{m} is 10 times stronger than the CRSS of baN{r/z} slip. When the relation is 1/1 (c, d), both slip systems are easily activated, whereas 10/1 (e, f) means that rhomb slip in baN is 10 times stronger than prism slip in baN.
or the same as b a N (c) and [c]{m} is hard. If then prism [c] slip is as weak as basal and rhomb b aN, then rhomb baN slip is mainly responsible for the strain accommodation, even when its relative CRSS is stronger than that of baN(c), when baN{r/z} becomes progressively more important with increasing deformation (Fig. 5). In the transition baN{r/z} → baN{m}, the slip system that is more likely to accommodate strain is b aN{r/z}, assuming that the relative CRSS of b aN{m} is harder or the same as that of baN{r/z}, no matter what the CRSS of prism [c] slip is (Fig. 8).
In simple shear, we observed a similar behavior of the slip systems related to strain accommodation, and basal b a N slip again accommodates better the deformation during the transition b aN(c) → b aN{r/z}, assuming that the relative CRSS of b aN{r/z} is higher or the same as that of b aN(c) and that [c]{m} is also hard (Fig. 12). If prism [c] is weak, basal and rhomb b a N slips have similar activities when the relative CRSS baN{r/z} is higher or the same asbaN(c), but prism [c] also accommodates a great part of the strain in these conditions (Fig. 12). In the transition b aN{r/z} → baN{m}, the slip system that is more likely to accommodate the deformation is rhomb ba N, assuming that the relative CRSS of b aN{m} is harder or the same as that for b aN{r/z}, and that [c]{m} is hard. If the relative CRSS of [c]{m} is as weak as baN{r/z} and b aN{m}, then it accommodates as much strain as the other active slip systems (Fig. 14). In addition, even slip systems that are theoretically inactive due to the chosen high relative CRSS may present some activity, probably because the other slip systems are geometrically unable to accommodate the imposed deformation. As an example, prism [c] and basal b aN slip in the Fig. 5e are assumed to be hard to activate. Together however, they accommodate more than 20% of the axial compression imposed deformation, possibly because rhomb b aN slip is unable to accommodate larger strains alone. Obviously, in the presence of other mechanisms such as recrystallization, single slip systems might dominate (e.g. Kilian et al., 2011), but such effects are not considered here.

The opening angle of quartz [0001] axes as a proxy for deformation temperature
In recent years, the opening angle of quartz [0001] axes started to be used as a proxy of deformation temperatures in quartz (e.g. Kruhl, 1996;Law et al., 2011;Sarkarinejad et al., 2009Sarkarinejad et al., , 2010. According to these authors, the opening angle increases with increasing deformation temperatures. If we assume that the transition from basal b a N → rhomb ba N → prism [c] slip is directly related to increasing temperatures, then as a first approximation this argument is valid. Nevertheless, detailed measurements of the opening angle between the small circle girdles and the axis of compression (Z) extracted from the quartz [0001] axis pole figures of our models demonstrate that the angle variation according to different sets of CRSS of these three slip systems is considerably high and non-linear (Fig. 16). For example, when prismatic [c] slip is harder than basal or rhomb baN slip, the angle variation produced by the transition b aN(c) → b aN{r/z} is rather stable and actually does not  {m}. The scale on the x-axis is given by the ratio baN(c)/baN{r/z}. When the ratio is 1/10, the CRSS for baN{r/z} is 10 times stronger than the CRSS of baN(c) slip. When the relation is 1/1, both slip systems are easily activated, whereas 10/1 means that basal slip in baN is 10 times stronger than rhomb slip in baN. Note the intense angle increasing when rhomb baN is assumed to be weaker than basal baN slip, and the important effect caused by the increasing importance of [c]{m} slip.
increase considerably until we assume that basal slip is harder than rhomb slip in ba N, when the angle increases quickly. In these conditions, the variation on the opening angle also increases with the progressive decrease of the relative CRSS of prismatic [c] slip, becoming broader with the increasing activity of the later. The direct measurements of the opening angles are also prone to errors, principally related to the trend that small circle girdle has from the center to the border of the stereonet, that might easily result in errors of more than 10°. For these reasons, alternative statistical methods based on calculations of eigenvalues calculated for pole figures, texture indexes and fabric entropy will be described in detail in the third part of this paper. The advantage of the statistical approach is that the indexes are calculated from the orientation matrix and do not depend on external measurements in the pole figures.

The lack of Y maxima in the pole figures
Quartz [0001] axis orientated parallel to the intermediary strain axis (Y-maximum) is a common feature in quartz bearing rocks deformed under relatively high temperature conditions (e.g. Bunge and Wenk, 1977;Law et al., 1992;Lister and Dornsiepen, 1982;Pennacchioni et al., 2010;Schmid and Casey, 1986;Schmid et al., 1981;Toy et al., 2008). In most cases, the Y-maxima fabric is attributed to a predominant activation of prismatic baN slip (e.g. Schmid and Casey, 1986). The simulations under axial compression and simple shear carried out here demonstrate that although a weak maximum parallel to Y can be developed under simple shear if prismatic ba N slip is dominant, this maximum is far from being as strong as the maxima observed in the papers above. The only clear Y-maxima observed in our simulations considering the maximum shear strain of 1.73 was observed when the prismatic ba N slip works alone (Fig. 2), where in addition to the single maxima close to Y, there is a dispersion of c-axes in the YZ plane. Such a pattern is quite similar to the ones reported by Pennacchioni et al. (2010) in non-recrystallized quartz ribbons deformed at 2 b γ b 3 in quartz veins in the Adamelo tonalite. This suggests that such a pattern may result from the activation of a single slip system rather than the activation of multiple slip systems with one dominant over the others (e.g. Kilian et al., 2011) in higher shear strains in comparison to the values used in our simulations. In addition, the single slip on prismatic ba N slip system is geometrically unlikely to produce the "single-crystal" type of orientation reported in the literature at relatively low strain conditions, as the one used here. Therefore, there must be a secondary deformation mechanism or recovery processes that lead to the development of such strong CPOs. Heilbronner and Tullis (2006) demonstrated experimentally that the development of Y-maxima might be caused by progressive recrystallization of quartz under high strain conditions. In their case, the deformation of coarse grain quartz aggregates leads to the complete recrystallization of the aggregates, and the variation of the [0001] axes resulted from the grain size growth and reduction of differently orientated crystallographic domains. The grain size of the recrystallized grains is crystallographically controlled, and the crystals with c-axes parallel to Y are between 1.5 to 2 times bigger than the crystals with other orientations. Similar features were observed by Jessell (1988a,b) also using numerical models, where the author shows the direct relationship between grain boundary migration mobility and the orientation of the quartz c-axes, whether in weak or hard orientation positions. Such a close relationship between Y-maxima and strain conditions were also observed in the high strain simulation in simple shear that are subject of the second part of this paper.

Conclusions
Numerical simulations of the crystallographic preferred orientation transition in quartz via viscoplastic self-consistent approach under axial compression and simple shear in constant strain allowed us to conclude the following: (i) CPOs are strongly developed in axial compression rather than in simple shear. The main cause for that is once the crystals acquire a stable position of "easy-slip" in relation to the imposed axial compression, the crystals do not rotate further, which happens in crystals in simple shear; (ii) A number of the CPOs predicted in the transitions from one dominant slip system to the other are common in naturally and experimentally deformed quartz aggregates. These include (i) small circle girdles of [0001] axes parallel to σ 1 in axial compression; (ii) constant increase of opening angle of these girdles with increasing activity of baN{r/z}, b aN{m} and [c]{m} slip systems in relation to b aN(c); (iii) single girdles of [0001] axes normal to the foliation plane when b aN(c) and baN{r/z} are dominant slip systems; (iv) [0001] maximum parallel to the lineation when [c]{m} slip system has a low relative CRSS. Nevertheless, some classical CPO patterns in quartz such as Y-maxima, or typical crossed girdles are just weakly developed in the simulated conditions. This suggests that they may result either from different deformation/recovery mechanisms, such as dislocation climb, dynamic recrystallization and associated processes or even static annealing, or from different imposed deformation conditions (i.e. higher strains, different imposed strains, such as transpression); (iii) While some of the predicted CPO patterns may be repeated over a large range of different sets of relative CRSS, other patterns are developed in very specific conditions. Considering that all the simulations carried out here were under constant strain and also under constant CRSS of the slip systems (no hardening was evolved during simulation), we suggest that while some combinations or different CRSS do not affect the general aspect of the CPO, other combinations may produce significantly different patterns. In our simulations, the effect of [c]{m} slip system seems to be the most pronounced in producing different CPO patterns at specific combinations of the other studied slip systems in quartz; (iv) Some of the slip systems used in our simulations seem to align more easily with respect to the imposed strain than the others, normally producing stronger CPOs. Considering the CPO transitions and the slip system activities, the baN(c) slip system is the one that produces stronger preferred orientations in quartz, followed by b aN{r/z} and [c]{m} and finally baN{m}. Although the latter does not produce strong CPOs when is the dominant slip system, the combination with [c]{m} results in strong CPOs, particularly in axial compression. For similar reasons, the transition baN(c) → b aN{r/z} in simple shear induces a progressive weakening of quartz CPO, suggesting that the former slip system is more easily aligned with the imposed shear strain than the latter; (v) In certain cases, slip systems that are theoretically inactive due to the chosen high relative CRSS may present some activity, most likely because the combination of weak slip systems is unable to accommodate all the impose deformation; (vi) Although the constant strain simulations cover in great detail the possible CPO patterns developed for the main slip systems in quartz deformed under axial compression, the same is not true for the simple shear simulations.