A Multiscale Numerical Modeling Investigation on the Significance of Flow Partitioning for the Development of Quartz c‐Axis Fabrics

Quartz c‐axis fabrics in natural mylonites can vary to such an extent that they apparently give opposite senses of shear in a single thin section. Many hypotheses have been invoked to explain this. Here, we couple our self‐consistent multiscale approach for flow partitioning with the viscoplastic self‐consistent model for c‐axis fabric simulation to investigate quartz c‐axis fabric development. Quartz aggregates are regarded as microscale Eshelby inhomogeneities embedded in a macroscale medium whose effective rheology is represented by a hypothetical homogeneous equivalent medium which is assumed to be rheologically isotropic or has a planar anisotropy. We reproduced the observed quartz c‐axis fabrics. We found that although the microscale flow fields are distinct from one another and from the macroscale flow, the microscale vorticity in every inhomogeneity has the same sense as the macroscale vorticity. This implies that one can use the average of the microscale vorticity axes determined through the crystallographic vorticity axis analysis to obtain the macroscale vorticity axis. However, quartz c‐axis fabrics cannot be used to determine the vorticity number where flow partitioning is significant.

strain increases. But this is not supported by any numerical modeling work with shear strains up to 5 (e.g., Jessell & Lister, 1990;Morales et al., 2011). Little et al. (2016) suggested that plane-strain general flows may produce synthetically-inclined c-axis girdles. The VPSC modeling of Takeshita et al. (1999) and Nie and Shan (2014) considered such flows but did not obtain any synthetically-inclined c-axis girdles. Their models did not consider cases where prism<a> and rhomb<a> slips are more significant than basal<a> slip. Through VPSC modeling, Keller and Stipp (2011) produced synthetically-inclined c-axis girdle with rhomb <a>, prism <a>, and prism <c> all active and all more significant than the basal <a> slip system. However, there are natural samples with synthetically-inclined c-axis girdles (e.g., Kilian et al., 2011;Law et al., 2010;Little et al., 2016) that were clearly produced in the temperature range ∼350°-550°C, much below that required for the activation of prism <c> slip (Toy et al., 2008). Furthermore, c-axis fabrics with both synthetically-and antithetically-inclined c-axis girdles have been observed in the same thin section (Kilian et al., 2011 their Figures 8 and 9 reproduced as Figure 2). Therefore, it is important to clarify if the combination of rhomb<a> and prism<a> slip, without prism <c>, can generate synthetically-inclined c-axis girdles.
In this contribution, we apply a multiscale numerical modeling approach to quantitatively investigate the consequence of flow partitioning on the development of quartz c-axis fabrics and compare our modeling results with natural and experimental observations. Specifically, we use the mylonite thin-section photomicrograph of Killian et al. (2011) as an example and seek to understand if flow partitioning can produce the observed quartz c-axis fabric variation. Figure 2a shows the thin-section photomicrograph. The mylonite is composed of quartz domains whose c-axis fabrics are shown in Figure 2b, feldspar porphyroclasts, and a matrix of fine-grained quartz, feldspar, and mica. According to Kilian et al. (2011), the microstructure was BHANDARI AND JIANG 10.1029/2020JB021040 2 of 17 Figure 1. Current models for interpreting quartz c-axis fabrics at low to medium temperatures ∼350°-550°C presented in the shear zone coordinate system. The sense of shear is dextral. Cross-girdle and single girdle c-axis fabrics are commonly observed. (a) c-axis cross girdle pattern with one girdle normal to the shear plane while the other antithetic to the shear sense. (b) a single c-axis girdle with Y-maxima, inclined either antithetically to the shear sense or normal to the shear plane. (c) a single c-axis girdle inclined either antithetically to the shear sense or normal to the shear plane. The c-axes near the periphery, the center, and in between are interpreted to reflect basal<a>, prism<a>, and rhomb<a> slips respectively. produced in a simple shearing flow, which is consistent with the geometric pattern (Ramsay, 1980;Ramsay & Graham, 1970) of the foliation. In the following, we consider plane-strain general shearing flows as the bulk flow field and investigate the partitioned local flows in the rheologically distinct quartz domains to understand the c-axis fabrics in them. Our method here also applies to any 3D general shearing flows (e.g., Jiang & Williams, 1998), but for the purpose of the present paper, considering plane-strain general shearing is sufficient.

Approach
Our approach is illustrated in Figure 3. We consider two different scales. The macroscale shear zone flow is a plane-strain general shearing flow, which in the coordinate system as shown in Figure 3a and is defined by the following Eulerian velocity gradient tensor: where  is the shear strain rate for the simple shearing component and  is the strain rate parallel to the shear zone boundary (x-axis). The flow in Equation 1 corresponds to a kinematic vorticity number Jiang & White, 1995;Li & Jiang, 2011;Truesdell, 1953). We consider the variation of the macroscale flow by varying k W from 0 to 1.
The shear zone contains numerous quartz domains which are rheologically distinct. The microscale flow refers to the local partitioned flow in an individual quartz domain which is distinct from the macroscale flow because each quartz domain is unique in its instantaneous rheology, shape, and orientation (Eshelby, 1957;Jiang 2014Jiang , 2016Mura, 1987). Clearly, it is microscale partitioned flows in quartz domains that are relevant for the quartz c-axis fabrics ( Figure 2b). We use our self-consistent Multi Order Power Law Approach (MOPLA) (Jiang, 2014(Jiang, , 2016Jiang & Bentley, 2012;Lu, 2020 andQu et al., 2016) to obtain the partitioned flow fields in quartz domains first and then the Viscoplastic Self-Consistent (VPSC) model of Lebensohn and Tomé (1993) to get the quartz c-axis fabrics in quartz domains. The procedure is described briefly below.
Each quartz domain is regarded as an ellipsoidal Eshelby inhomogeneity embedded in and interacting with the macroscale material ( Figure 3c). The rheology of the latter is approximated by a homogeneous BHANDARI AND JIANG 10.1029/2020JB021040 3 of 17   Jiang, 2014Jiang, , 2016. The partitioning and homogenization equations are solved simultaneously so that self-consistent solutions to the partitioned flow fields and the macroscale rheology are obtained. As the relevant rheology for mylonites are power-law viscous (Kohlstedt et al., 1995), the MOPLA formulation adopts a linearization approach (e.g., Lebensohn & Tomé, 1993;Molinari et al., 1987) where linearized viscosities such as tangent viscosities are used in the formulation.
We consider two situations for the rheology of the quartz domains and HEM. In the first, both the quartz domains and the HEM are assumed to be isotropic ignoring the rheological anisotropy development in HEM as a result of fabric buildup with strain. This simplification is because anisotropic rheological response of the constituent quartz domains is not available. Because of this simplification, the rheological contrast between a quartz domain and the HEM is reduced to an effective viscosity ratio r between the domain to the HEM. Owing to power-law rheology, r varies with time. We thus consider a range of constant r values. The quartz domains in Figure 2a all have convex shapes with surrounding matrix material wrapping around them, suggesting that quartz domains were rheologically stronger (r > 1) than the ambient HEM. But r cannot be too high (e.g., r > 10) because of power-law rheology or the quartz domains would behave like rigid clasts (Jiang, 2007a;Xiang & Jiang, 2013)   tion. We consider the situations of r being 2, 5, and 10. The situation of r = 0.5 is also considered here for comparison to show what the c-axis fabrics might be like if quartz domains were mechanically weaker than the ambient medium.
In the second situation, we consider isotropic quartz domains in a HEM of simple planar anisotropy which approximates a foliation and/or layering common in natural mylonites. The rheology of such a HEM can be characterized by two distinct viscosities:  n, the normal viscosity for the resistance to pure shearing along and perpendicular to the layering, and  S the shear viscosity measuring the resistance to shearing parallel to the layering (e.g., Fletcher, 2009;Jiang, 2016;Johnson & Fletcher, 1994). The strength of anisotropy is measured by the ratio m of  n to  S. For foliated and layered rocks, m > 1 (Treagus, 2003). The rheological contrasts between the isotropic quartz domain and HEM can be defined by the following two parameters: the ratio, r eff , between the viscosity of the domains to  n and m. In such case, the effective viscosity of the domains is simply given by r eff  n. Similar to the isotropic cases, we consider r eff being 0.5, 2, 5, and 10. Figure 4 shows the geometric relation between the flow field and the plane of anisotropy. Macroscale flow is simple shear with shear plane parallel to x-z plane.
To cover the shape variation of quartz domains, we considered three reference initial shapes: prolate (5:1:1), oblate (5:5:1), sphere (1:1:1), and initial triaxial domains with long and short semi-axial length fixed to 5 and 1, respectively, and intermediate semi-axial length ranging from 2 to 4. These initial domains will deform into various possible triaxial shapes in nature. The orientations of these domains are defined by spherical angles and the initial orientations are generated using the method described in Jiang (2007aJiang ( , 2007b) so that they exhibit uniform random distribution in 3D space.
We use the Viscoplastic Self-Consistent (VPSC) model, originally due to Molinari et al. (1987) and further developed by Lebensohn and Tomé (1993), to simulate the quartz c-axis fabric development in a given flow field. Specifically, VPSC7 (Lebensohn & Tomé, 2009) for Windows is used in this investigation. We track the c-axis evolution of 500 quartz crystals with initial orientations being uniform randomly distributed in 3D space (Jiang, 2007b). The initial shapes of quartz grains are equant. The deformation of the grains is achieved by intracrystalline slip on slip systems that follow a rate-sensitive viscoplastic behavior (relevant equation and model parameters are provided in the supplementary material). The crystal shapes evolve with strain and the grain fragmentation scheme of Beyerlein et al. (2003) is used to limit the aspect ratio of quartz grains to mimic some effect of dynamic recrystallization. While other mechanisms such as recrystallization may affect the efficiency of c-axis fabric development, the c-axis fabric patterns are mainly determined by dislocation glide (Llorens et al., 2016). As the objective of this work is on effect of flow partitioning on quartz c-axis fabric, we consider dislocation creep only. VPSC output of c-axis data are plotted using the MTEX toolbox (Bachmann et al., 2010). Since we are concerned with variation in quartz c-axis fabrics, we present pole figures of the (0001) directions here. Pole figures of other crystal directions are provided in the supplementary material.
The relative activity of a slip system is modulated in VPSC by their relative critical resolved shear stress (CRSS). A lower CRSS corresponds to higher activity. For quartz, the CRSS for a slip system is largely temperature dependent with basal <a> slip occurring at about 350°C-500°C, rhomb<a> and prism<a> slip at ∼500°-600°C, and prism<c> slip at >650°C (e.g., Toy et al., 2008). Table 1 summarizes the slip system combinations used in this study based on previous work (e.g., Lister & Paterson, 1979;Morales et al., 2014;Wenk et al., 1989). They are labeled as Model-A to F. Morales et al. (2014) and Keller and Stipp (2011) made a distinction between rhomb (+) slip{r} in <a> direction and rhomb (−) slip{z} in <a> direction in their modeling. This distinction may be important in regime 1 dislocation creep (Hirth & Tullis, 1992) but it is not necessary for our work here on the effect of flow partitioning on quartz c-axis fabrics.
If one considers quartz c-axis development in a single-scale plane-strain general shear, then the macroscale flow such as the one defined in Equation 1 is used directly with VPSC. To consider the development of quartz c-axis fabrics due to partitioned flows in quartz domains, calculated from MOPLA, a coupled computation between MOPLA and VPSC is necessary. This is because the partitioned flow field in any given quartz domain is non-steady as the domain continuously changes shape and orientation during defor-mation. The coupling is carried out as follows: With a given macroscale flow field defined in Equation 1, we use MOPLA algorithm implemented in MATLAB (Jiang, 2007a(Jiang, , 2014(Jiang, , 2016Jiang & Bentley, 2012;Lu, 2020;Qu et al., 2016) to calculate the partitioned flow in every quartz domain, which is expressed as a velocity gradient tensor. We export to a data file the domain's velocity gradient tensor for every prescribed macroscale strain increment until a given macroscale finite strain is reached. This data file is then used as input flows into the VPSC code to calculate the c-axis fabric evolution within the quartz domains as macroscale strain increases until the set magnitude. The velocity gradient tensor files for all the domains are available in the supplementary material.
The progression of the finite strain is measured by a strain intensity de- In the event of simple shearing ( k W = 1), the strain can also be measured by the shear strain γ parallel to the zone boundary. The relation between γ and ρ in simple shearing situation is shown in Figure 5.

Results
As mentioned above, we have simulated quartz c-axis fabric development in three situations. The first is in homogeneous macroscale plane-strain general shearing flows, without flow partitioning, with quartz slip system combinations that have not been covered by previous studies. The second situation is when rheologically distinct quartz domains are within an isotropic HEM, and the third is when the quartz domains are within a HEM of planar anisotropy.

Quartz C-axis Fabric Development in Homogeneous Plane-Strain General Shearing Flows
Figures 6 and 7 present quartz c-axis fabrics produced for models-A-F (rows) under uniform macroscale flows of plane-strain general shearing from k W = 0 to k W = 1 (columns), at macroscale strain states ρ = 2 and 6.
In pure shearing ( k W = 0), for models A, C, and D, peripheral c-axis maxima form and remain at the maximum shortening direction regardless of strains (Figures 6a, 6d, 6m, 6p, and 7a, 7d). For model-B, a cross-girdle ( Figure 6g) form at ρ = 2, with its central segment lying along the maximum shortening direction. A single girdle is produced lying along the maximum shortening direction (Figure 6j) at ρ = 6. For models-E and F, peripheral maxima develop and remains at the maximum stretching direction at both strain states (Figures 7g, 7j, 7m, and 7p).
When 0 < k W ≤ 1, for models-A, C, and D and at ρ = 2, peripheral c-axis maxima are developed inclining antithetically to the shear sense (Figures 6b, 6c, 6n, 6o, 7b, 7c). As k W increases, the angle between the peripheral c-axis maxima and the shear plane normal increases. The peripheral maxima rotate with macroscale vorticity as finite strain increases but do not pass the shear plane normal (Figures 6f, 6l, 6q, 6r, and 7f). For model-B, a cross-girdle, with its central segment lying normal to the shear plane, is formed at ρ = 2 (Figure 6h, 6i). At ρ = 6, this cross-girdle becomes a single girdle lying normal to the shear plane (Figures 6k,  6l). In some general shear (0 < k W < 1) cases, the c-axis girdles can be slightly synthetically-inclined at ρ = 6 (Figures 6e and 6k,7e), but the angle between the peripheral c-axis maxima and the shear plane normal is small (<10°). For models-E and F, the c-axis peripheral maxima are synthetically inclined near the shear direction at ρ = 2. The angle between the peripheral c-axis maxima and the shear direction increases with k W (Figures 7h, 7i, 7n, 7o). The c-axis peripheral maxima rotate with vorticity toward the shear direction as the finite strain increases (Figures 7k, 7l, 7q, 7r).

Quartz C-axis Fabric Development in Quartz Domains Embedded in an Isotropic HEM
Since models-A, B, C and D all produce similar c-axis fabrics with c-axis girdles antithetically-inclined or nearly normal to the shear plane, we only present results for model-A for multiscale deformation. Figures 8 and 9 report 10 results of c-axis fabrics produced in quartz domains of varying initial shapes, orientations, and viscosity ratio r under a range of macroscale flow fields. These results can be summarized as follows: When 0 < k W ≤ 1, c-axis girdles are always antithetically inclined (Figures 8a, 8b, 8d, 8e, 8g, 8h, 8j, 8k, 8m, 8n, and 9a, 9g, 9j, 9m) regardless of initial conditions of quartz domains, up to the finite strain ρ ∼ 4. With increase in finite strain, the girdles rotate with bulk vorticity (rows of Figures 8 and 9) but do not pass the shear plane normal unless r >= 5 (Figure 8 -i, 8l, 8o, 9-b, 9c, 9i, 9l, 9o). If the quartz domains were weaker than the HEM (r = 0.5), c-axis girdles remain close to normal (Figure 8c) to the shear plane even at very high finite strains (ρ ∼6-7). In pure shearing, a cross girdle is produced at ρ ∼2 which becomes a single girdle at ρ ∼6 with peripheral c-axis maxima always parallel to the maximum shortening direction (Figures 9d, 9e, 9f). Figure 10a-o reports c-axis fabrics developed in quartz domains of varying initial shapes, orientations, and viscosity ratio r eff , embedded in a planar anisotropic HEM of anisotropic strength m as described in Section 2. These results can be summarized as follows: The c-axis girdles are always antithetically-inclined (Figures 10a, 10d, 10g, 10j, and 10m) at ρ ∼2 regardless of k W , r eff , m, initial orientations, and shapes of domains. The girdles rotate with bulk vorticity as ρ increases (rows of Figure 10) but do not pass the shear plane normal unless r eff >= 2 (Figures 10i, 10l, 10o). If quartz domains were weaker (r eff = 0.5), c-axis girdles are close to normal (Figures 10c and 10f) to the shear plane even at high finite strains (ρ ∼6-7).

Discussion
In Section 2, we argued that r must be greater than 1 from microstructures ( Figure 2a) but not so high that the quartz domains do not develop enough internal strain for c-axis fabric formation. In our modeling, we considered the range r between 2 and 10. Our modeling results show that, with basal<a>, rhomb<a>, and prism<a> slips, c-axis girdles in a quartz domain always develop at an antithetical orientation initially. But the girdles rotate with vorticity as the macroscale finite strain increases. If the domain is sufficiently strong (  5 r in isotropic HEM case and r eff ≥ 2 in planar anisotropic HEM case), the girdles will rotate pass the shear zone normal and lie in the synthetical sector at high strains (ρ ∼6, Figure 8i , 8l, 8o, 9-b, 9c, 9i, 9l, 9o, and 10-i, 10l, 10o). Our modeling results show that for synthetically-inclined girdles, r should be between 5 and 10. These results are consistent with the strain gradient of the thin section (Figure 11, based on Killian et al., 2011, Figure 9 there). The antithetically inclined c-axis peripheral maxima (yellow in Figure 11) correspond to a lower strain and synthetically-inclined c-axis peripheral maxima (red in Figure 11) to a higher strain.
Despite the variability of the partitioned flow field in quartz domains, we found out that for the viscosity ratio range used in this investigation (r ∼ 2-10), which is appropriate based on microstructures, the microscale vorticity in every quartz domain still has the same sense as the macroscale vorticity. Figure 12 shows the dot product of the unit vector ω parallel to the vorticity in a quartz domain and the unit vector Ω parallel to the macroscale vorticity.  ωΩ is positive for all quartz domains (a few selected ones are shown in Figure 12). In other words, there is no vorticity sense reversal in any quartz domains. This implies that the volume-weighted average of microscale vorticity vectors is parallel to the macroscale vorticity vector. Therefore, the crystallographic vorticity axis (CVA) analysis (Giorgis et al., 2017;Michels et al., 2015) can still be used in every quartz domain and then the averaged microscale vorticity axes represent the macroscale vorticity axis.

Model-F
Synthetically inclined c-axis girdles have also been reported in some creep experiments on pure quartz aggregates with similar slip systems (basal<a>, rhomb<a>, and prism<a>) at high finite strains (Heilbronner & Tullis 2002;2006). Although Keller and Stipp (2011) obtained synthetically-inclined c-axis girdles by including prism<c> slip in their VPSC models, we have further confirmed that without prism <c> slip (our model-D), synthetically-inclined c-axis girdles cannot be produced (Figures 7a-7f), unless flow partitioning is considered. We suspect that some degree of heterogeneous strain and therefore partitioned flow was responsible for such c-axis girdle orientations. Heilbronner and Tullis (2006)  a a a ), and initial orientation given by spherical angles (θ 1 , Φ 1 , θ 2 ) for general domains or (θ, Φ) for spheroidal domains]. Each row presents the results for the domain as the macroscale strain increases. The c-axis fabric strength is evident from the pole densities expressed in multiples of uniform distribution (m.u.d.) and shown at the right of each plot. At any given finite strain intensity, the m.u.d. values are higher in the weaker quartz domains (r < 1) than the stronger quartz domains (r > 1). This is consistent with the fact that weaker quartz domains are deformed more intensely and develop a stronger c-axis fabric. With increase in finite strain intensity, the m.u.d. values increases.  strated by our modeling based on microstructures of natural mylonites raises further issues with using quartz c-axis fabrics to estimate the (macroscale) vorticity (e.g., Law, 2010;Vissers, 1989;S. R. Wallis, 1992S. R. Wallis, , 1995Xypolias, 2009  a a a ), initial orientation given by spherical angles (θ 1 , Φ 1 , θ 2 ) for general domains or (θ, Φ) for spheroidal domains, and anisotropic strength m]. Each row presents the results for the domain as the macroscale strain increases. The c-axis fabric strength is evident from the pole densities expressed in multiples of uniform distribution (m.u.d.) and shown at the right of each plot. At any given finite strain intensity, the m.u.d. values are higher in the weaker quartz domains (r < 1) than the stronger quartz domains (r > 1). With increase in finite strain intensity, the m.u.d values increases in all models. that the assumption commonly used in vorticity determination that the dominant c-axis girdle is perpendicular to the shear plane is not always valid.
Peripheral c-axis maxima close to the shear direction have commonly been taken to reflect prism <c> slip (Passchier & Trouw, 2005). Our modeling suggests that they can also be significantly rotated peripheral basal <a> maxima from certain quartz domains (Figure 9c). Larson et al. (2014) have reported a possible example of this. They presented peripheral c-axis maxima close to the shear direction in a temperature condition much below that required for the activation of prism <c> and used the concept of flow partitioning to explain their observation. Our modeling lends support to this explanation.

Conclusions
The co-existence of both synthetically inclined and antithetically-inclined quartz c-axis girdles in a single thin section can be explained by flow partitioning at the thin section scale. The antithetically inclined girdles correspond to relatively low finite strains and the synthetically inclined girdles to high finite strains.
Although the microscale flow fields vary from one quartz domain to another and are distinct from the macroscale flow, the sense of vorticity in every quartz domain remains the same as the macroscale vorticity.
Because of flow partitioning, it is not possible to estimate the vorticity number of the macroscale flow from quartz c-axis fabrics. But it is still possible to obtain the macroscale vorticity axis by averaging the microscale vorticity axes from quartz domains. The latter can be obtained through the crystallographic vorticity axis analysis.
As a result of partitioned flows, the dominant quartz c-axis girdle can lie antithetical, normal, synthetical to the shear plane. The basal <a> peripheral maxima may end up lying close to the shear direction at high macroscale strains. Caution should be taken not to misinterpret these peripheral maxima as reflecting prism<c> slip.  Figure 9 there). C-axis fabrics comprise of peripheral c-axis maxima that are antithetically-inclined (yellow), nearly normal (blue), and synthetically-inclined (red) to the shear plane. The sense of shear is dextral. Slightly curved lines around the quartz domains traces the matrix foliation. Dotted green arrow at the left shows the direction of strain gradient. Figure 12. Plot of dot product  ωΩ with increasing macroscale strain ρ, for 50 quartz domains with random initial shapes and orientations, r = 0.5, 2 and 5 (rows) and W k = 1 and 0.75 (columns). ω and    are, respectively, the unit vectors parallel to the microscale vorticity vector in a quartz domain and the macroscale vorticity vector.  ωΩ is close to unity for the domains. Therefore, the microscale vorticity vectors are always of the same sense as and nearly parallel to the macroscale vorticity vector.

Data Availability Statement
Numerical modeling in this paper used our own MOPLA code (Jiang 2007a(Jiang , 2014(Jiang , 2016Jiang & Bentley, 2012;Qu et al., 2016) coupled with the VPSC code kindly provided by R. Lebensohn (Lebensohn & Tomé, 2009). Visualization used the MTEX toolbox of Bachmann et al. (2010). All model parameters are listed in Tables 1 and Section 2 of the text. VPSC input files and their description are provided as supplementary information.