Constraints on the emplacement of Martian nakhlite igneous rocks and their source volcano from advanced micro-petrofabric analysis

The Martian nakhlite meteorites, ormed from a single magma source region, and emplaced during multiple events spanning at least 93 ± 12 Ma, represent a key opportunity to study the evolution of Martian volcanic petrogenesis. Here 16 of the 26 identiﬁed nakhlite specimens are studied using coupled large area electron backscatter diﬀraction (EBSD) and energy dispersive X-ray spectroscopy (EDS) mapping to determine shape preferred orientation (SPO) textures of contained augite (high Ca-clinopyroxene) phenocrysts by considering crystallographic preferred orientation (CPO). Textural parameters derived from EBSD and EDS analyses were used to calculate maximum and minimum magma body crystallization thicknesses via three endmember emplacement scenarios: thermal diﬀusion, crystal settling, and crystal convection. Results from CPO textural analyses indicate weak to moderate fabric textures that are comparable to those in terrestrial clinopyroxenites. In all samples, a consistent foliation within the { 001 } axis of augite is observed. In all but two of the studied nakhlites this { 001 } foliation is typically coupled with a weaker lineation fabric in one or more of the { 100 } , { 010 } , { 001 } axes. Results from the calculated magma body thicknesses are consistent with an emplacement mechanism for the nakhlites driven by crystal settling. These crystal settling results infer magmatic body thicknesses ranging from < 1 m to several 10’s m, forming two distinguishable groups that appear random when assessed against observed texture, geochemical, and age parameters. Coupled textural and modelling results therefore suggest that the nakhlite source volcano varied in thickness over time yet consistently solidiﬁed via the mechanism of crystal settling.


Introduction
Mars is a planet whose landscape has developed to a large degree by igneous processes (Carr & Head, 2010;Greeley & Spudis, 1981;Grott et al., 2013;Taylor, 2013) and modified by hypervelocity impacts (Carr & Head, 2010;Zuber, 2001).From remote sensing studies of the volcanic provinces present on Mars, it is clear that volcanism on the planet has been long lived spanning several billions of years (Lapen et al., 2017;Lapen et al., 2010;Werner, 2009), although observation of active Martian volcanism has yet to be observed since orbital monitoring of the planet began in 1964.Despite the wealth of knowledge gained from remotely analysing Martian volcanoes, there is a limit to our understanding of the magmatic processes on Mars that can be gained from remote sensing alone (e.g., constraints on interior volcanic processes, magma chamber evolution, quantitative assessment of magma composition, determination of physical parameters such as emplacement textures and compositional variability) that can only be achieved through the analysis of physical specimens.
The Martian meteorites are the only Martian materials that are currently available for laboratory-based examination (McSween & Treiman, 1998;Udry et al., 2020).These samples encompass a range of different Martian rocks (e.g., the various types of shergottites, nakhlites, chassignites, an orthopyroxenite, and polymict breccias) that are interpreted to have been ejected from at least 11 different sites on Mars, probably from the southern highlands (Udry et al., 2020).
However, despite being able to group these samples chemically, isotopically, and texturally, the location of any ejection sites on Mars, their source volcano(es), and their emplacement style (e.g., flow/intrusion, sub-aerial/hypabyssal) are still unknown.Knowing how these meteorites were emplaced (flows versus intrusions) can provide information relating to locating their launch crater on Mars.Identifying patterns in emplacement mechanisms can provide information about whether the Martian crust (as sampled by these meteorites), through volcanic activity, has the potential to act as heat source capable of generating and sustaining liquid water on relevant timescales to create potential habitable regions for the development and sustenance of life (Kurokawa et al., 2014;McSween et al., 2001).
The nakhlites are currently the largest group of Martian meteorites derived from a singular parental magma source interpreted to be a single volcano on Mars (Udry et al., 2020).
Currently, there are 26 known individual specimens (several of which are considered to be paired) representing at least four distinct magmatic events spanning 93 Ma (1416 ± 7 to 1322 ± 10 Ma; Cohen et al., 2017).Originating as basaltic cumulates evolved from the same magma source region, the nakhlites are the only group of igneous rocks, with the exception of Allan Hills (ALH) 84001, to exhibit evidence of Martian fluid/rock interaction (Lee et al., 2018;Treiman, 2005).However, it is unknown whether the nakhlites are sampling a volcanic source which had a regular eruption record or a history consisting of longer dormant phases with condensed bursts of magmatic activity, and whether the meteorites originated solely as surficial flows, only intrusions, or as a combination of both.manuscript submitted to Journal of Geophysical Research: Planets Despite geochemical and isotopic data being excellent at identifying magmatic sources, distinguishing between individual lava flows derived from a single volcano or volcanic province can be difficult to discern in cases where there is little chemical or mineralogical variation between units, unless additional information is obtained through textural analysis, geochronological data, or geological context (Fenton et al., 2004).Textural properties of an igneous rock record the physical processes of a magmatic event from source to emplacement (Jerram et al., 2018).Thus, by analysing textural similarities and/or changes between multiple magmatic units sourced from a single magmatic source, important information on the evolution and history of the magmatic source can be discerned.
Recent EBSD measurements on Governador Valadares, Lafayette, Miller Range (MIL) 03346, and Nakhla by Daly, Piazolo, et al. (2019) showed evidence for gravitational settling in all four samples and weak flow textures in Governador Valadares and Nakhla, indicating potential complexity within the nakhlite launch crater.Daly, Piazolo, et al., (2019) hypothesised at least three distinct magmatic systems encompassing two regimes (subaerial hyperbolic flow and crystal settling) are represented within the nakhlites.Here we use electron backscatter diffraction (EBSD) data to substantially expand on the work of Daly, Piazolo, et al. (2019), assessing petrologic textures formed during emplacement and generating input parameters to model magmatic body thicknesses for 16 known nakhlite meteorite stones.The acquired data were then used to investigate what variation and/or trends can be discerned from observed emplacement mechanism(s) and CPO textural fabric(s) to better understand whether the hypothesised diversity of magmatic systems within the nakhlites proposed by Daly, Piazolo, et al., (2019) are applicable across an expanded nakhlite suite or to just a few individual samples.

Electron Backscatter Diffraction (EBSD)
Using electron backscatter diffraction (EBSD) a combination of premade and newly made polished nakhlite thick sections (21 sections encompassing 16 individual stones; Table 1) were analysed.The samples were chosen to capture the growing diversity of the nakhlite suite considering the recovery location, known crystallisation age, and model mineralogy.All sections were cut with random orientations i.e., without respect to any known/inferred petrofabrics.For both for consistency and to ensure comparability between the different sections studied, the principal orientations of the samples were defined based on the EBSD scans: X = left-right direction of the thick section, Y = top-bottom direction of the thick section, and Z = direction perpendicular to the plane of the thick section.Data used in our analyses include previously presented EBSD datasets which can be found in Daly, Lee, et al., (2019), Daly, Piazolo, et al. (2019), and Lee et al., (2018).
Prior to EBSD analysis, all sections were mechanically and chemically polished.
Mechanical polishing was achieved using 1 µm and then 0.3 µm aluminium spheres suspended in glycol for 5 minutes each.Subsequent chemical polishing was achieved using 0.1 µm colloidal silica suspended in a NaOH solution for 4 hours.Following polishing, all samples were coated by ~10 nm thick conductive carbon coat using a sputter coater.manuscript submitted to Journal of Geophysical Research: Planets Coupled EBSD and energy dispersive X-ray spectroscopy (EDS) analysis was run on four different instruments: a Zeiss Sigma Field Emission Gun Variable Pressure Scanning Electron Microscope (FEG-VP-SEM) operating Oxford Instruments AZtec analysis software v3.3 (ISAAC imaging centre, University of Glasgow); a Carl Zeiss EVO SEM using a HKL NordlysNano high Sensitivity EBSD detector [Geochemical Analysis Unit (GAU), Macquarie University], a Hitachi SU70 FEG-SEM equipped with a Symmetry CMOS detector and indexed using Aztec analysis software v3.4 (Oxford Instruments Nanoanalysis HQ, High Wycombe), and a Tescan MIRA3 VP-FEG-SEM with the NordlysNano EBSD detector and Aztec EDS/EBSD acquisition system (John de Laeter Centre, Curtin University).All analyses were run at 20 KeV, 4-8 nA beam current, 20 keV, 70° tilt, aperture of 120 µm, under high vacuum (~3.5 x 10 -4 Pa) apart from Lafayette (USNM 1505-1) and MIL 03346 (118), which were run at low vacuum (~ 49 Pa).Selected step sizes for each sample were chosen to maximise the area covered by the EBSD maps and ensure data collection over available timeframes.The detailed analysis settings for each analysed section, including step sizes (ranging from 2 To 15 µm), can be found in Table 1.
All data were processed using Oxford Instruments HKL Channel 5 software.Noise reduction was done using a wildspike correction followed by a consecutive 8-6 point nearest neighbour zero solution reduction to facilitate crystal definition without generating significant artefacts within the data set and to remove erroneous, mis-indexed and non-indexed, data points (Bestmann & Prior, 2003;Daly, Lee, et al., 2019;Forman et al., 2016Forman et al., , 2019;;Tedonkenfack et al., 2021;Watt et al., 2006).The noise-reduced data were then further processed through the one point per grain (OPPG; i.e., one point per crystal) reduction scheme to reduce bias from highly fractured crystals prior to crystal preferred orientation (CPO) and SPO analysis.Crystals <20 µm 2 in surface area were also removed from the dataset as they contained too few pixels (~<10) to robustly sample the mineral's SPO and CPO.The boundary of <20 µm 2 was determined by the scan resolution, step size, and shock levels within the sample (Forman et al., 2019;Watt et al., 2006).Crystal grain boundaries (defined as >10° internal crystallographic misorientation from the nearest-neighbour pixel) within the datasets were determined using Channel 5's automated "grain detect" algorithm.These crystal boundaries were further processed to exclude twin boundaries (180° rotation) identified to occur along augite's (100), ( 001), (204), or (104) axes.Thus, Channel 5 derived crystal lists can be used to interpret and quantify crystal textures within the different scans.To assess textural fabric there are three different quantitative metrics J-Index, M-Index, and eigenvalues (PGR) which all assess the orientation distribution function (ODF) of the crystals in a slightly different manner.J-Index uses the euler angles (Bunge, 1982), M-index assesses the crystal axes with equal weighting (Skemer et al., 2005), and eigenvalues assess the ODF of each crystallographic axis individually between three endmembers (Vollmer, 1990).For calculating the quantitative metrics and providing crystal statistics for thickness calculations raw CTF files were imported into MTEX.In MTEX the data underwent denoising (equivalent to Channel 5's wildspike), before establishing 10° misorientation as the crystal boundary and identifying twin crystals.The detected crystals were then separated into two groups; crystals <200 µm were classed as mesostasis crystals, and crystals >200 µm were classed as phenocryst crystals.
manuscript submitted to Journal of Geophysical Research: Planets 2.2 Shape Preferred Orientation (SPO) SPO represents the orientation relationships of the crystals (i.e., texture) within a given sample.The study of Daly, Piazolo, et al. (2019) showed that a coupling exists in augite between the long crystal-shape axis and the {001} also referred to as the <c> crystallographic axis within the nakhlites Governador Valadares, Lafayette, MIL 03346, Nakhla, Y 000593, Y 000749, and Y 000802.Thus, for SPO analysis, the slope of the crystallographic {001} axis of the OPPG subset is used as a proxy for the long crystal-shape axis within the sections to assess the 3D crystal textures.The oriented crystals from the OPPG datasets provides consistency across the multiple specimens analysed, removes bias against the lack of geological context or known sample orientation (e.g., Fig. 1), as well as bias from crystal fracturing within each of the sections.The slope angle was calculated using Oxford Instruments Channel 5 Tango module best-fit ellipse algorithm, which was also used to calculate crystal area and length of the long crystal-shape axis.

Crystallographic Preferred Orientation (CPO)
CPO refers to the preferred orientation of the crystallographic lattice axis for a selected phase within a sample.When a given crystallographic axis is linked to a given shape axis within a given crystal (e.g., long-shape axis), it can also be used as a proxy for SPO analyses.CPO was visually assessed from the OPPG dataset (one pixel/crystallographic orientation defining a crystal), as the OPPG reduction removes bias from fractured crystals within the dataset.Lower hemisphere equal area projections contoured using a maximum multiple uniform density (MUD.; representing the density of data points) of 5 with 5° clustering and a half width of 15° were used to plot {100}, {010}, and {001} crystallographic axes.CPO petrofabric intensity and/or presence was quantified from the phenocryst group of crystals using an adapted MATLAB MTEX code from Daly, Piazolo, et al. (2019), adapted from Vollmer (1990) (Bunge, 1982;Skemer et al., 2005;Vollmer, 1990).For each quantitative metric CPO strength was determined using the same comparative terrestrial equivalent parameters as stated in Daly, Piazolo, et al. (2019).Briefly, Jindex values and their associated fabric strength were assessed as: low fabric = 1.40-1.80,lowmedium fabric = 1.80-3.00,medium fabric 2.40-5.00,medium-strong fabric = 4.00-12.00and strong fabric <12; M-index: 0 = random orientation and 1 = single crystal (note that although M-Index calculations work well for many different minerals, clinopyroxene has been shown to produce weak M-Index textures irrespective of the overall classified fabric strength of the rock (Hidas et al., 2019), thus the clinopyroxene calculated values will be compared relative to each other rather than the standard index range); and Eigenvalues: random (R) = >90%, weak point (P) and girdle (G) fabric = 10-30%, moderate P and G fabric = 30-50 %, and strong P and G fabric = >50%.Where possible for J-Index and M-Index calculations only the largest 300 crystals were used to enable more accurate comparison between the different sized scans.For J-Index results datasets of 100-150 crystals are considered to be statistically stable (Ismaïl & Mainprice, 1998), while M-Index calculations are considered to have a stabilised uncertainty of 2% in datasets >150 crystals (Skemer et al., 2005), eigenvalue assesses the textural symmetry and axis orientation as a fraction between three endmembers where the combined PGR value for each crystallographic axis equals 1 (Mainprice et al., 2015;Vollmer, 1990).Common eigen index parameters, BA and LS were calculated using the following equations: (1) These index parameters assess the shape of the eigen fabric between perpendicular crystallographic axes giving a range between 0 and 1 [i.e., {010} (B or L) lineation = 1 against either {100} (A) or {001} (S) foliation = 0].

Calculating magma body thicknesses
Rough magma body thicknesses (L) for each dataset were calculated using MATLAB with MTEX toolbox (Bachmann et al., 2011), for three endmember mechanisms.These three equations of thermal diffusion (3), crystal settling (4), and crystal convection (5) were selected to provide the maximum possible thickness, intermediary thickness, and minimum thickness, respectively for each meteorite: In equation ( 3) κ is the thermal diffusivity of basalt (7.5 x 10 -7 m/s; Durham et al., 1987).In equations ( 3), (4), and (5), T is the residence time taken from (Krämer Ruggiu et al., 2020;Udry & Day, 2018).In equation ( 4) and ( 5), θ is the crystal volume fraction (calculated from the EBSD detected mineralogy), Δρ is the density difference between the melt and solid lava manuscript submitted to Journal of Geophysical Research: Planets (calculated from CIPW normative EDS spectra), g is the Martian gravity, a is the crystal diameter (from MTEX calculated phenocryst grouped crystal), and µ is the viscosity (calculated from CIPW normative EDS spectra).In general, thermal diffusion and crystal convection estimate the possible maximum and minimum thickness of the magmatic body.
To calculate CIPW normative spectra, large area spectra encompassing the EDS map areas were collected.The collected spectra was then processed following standard CIPW normative procedures (Bickel, 1979;Pruseth, 2009).The CIPW normative results were then used to calculate different model input parameters such as density and viscosity.Liquid magma temperatures used for density and viscosity estimates were calculated using the linear equation of McBirney (1993), which bases temperature from silica content within the sample.These estimated temperatures ranged between 1043 and 1121°C (For CIPW results, full code, and input parameters see supplementary materials).Among the set of analysed nakhlite meteorites, there is a range of SPO from moderate (one dominant direction of alignment with a smaller less dominant direction e.g., NWA 12542) to random (no alignment e.g., NWA 11013; Fig. 4).Moreover, the data show an average crystal orientation spread (GOS) difference of 113° ± 7°, ranging from 0-101° in MIL 090032, 108 to 0-128° in NWA 10153, SH65 T2-2 (Fig. 5, Table 2).Within the 21 individual datasets a low percentage of augite crystals were observed to have any of the three crystallographic axes aligned parallel (±5°) to the cut surface of the thick sections (Fig. 6,      2).While the CPO metrics for NWA 817, N8-1 are included, the results are not discussed and interpreted owing to the large associated error (e.g.Fig. 7).
Both replicate sections of Y 000593 exhibited low intensity (~0.3) weak point fabrics (P = 0.13) within {010} instead of {100} (Fig. 8).indicating S type to LS type fabrics for all analysed samples (Fig. 9, Table 2).For a thermal diffusion emplacement mechanism (equation 3), calculated magma body thicknesses were observed to sit within the tens of metres scale, where all calculated values were within error (Table 3).For a crystal settling emplacement mechanism (equation 4), magma body thicknesses were observed to span several different orders of magnitude (  NWA 817, NWA 12542).From this emplacement mechanism two distinct groups outside of error can be identified: group 1 being <10 m and group 2 being > 10 m.Lastly, the calculations for the crystal convection emplacement mechanism (equation 5), provided unrealistic thickness values in the range of tens to hundreds of nanometres (Table 3).The thicknesses as calculated via equations 3-5 are very different to each other, the reasons for these differences are further evaluated in section 4.2 below.No trends were able to be discerned when assessing calculated magma body thicknesses against modelling input parameters (e.g.,

Nakhlite textural observations
Textural observations, here defined as the geometric relationships between the sample's constituent crystals, provide insights into important petrological processes such as crystallisation growth rates, magma mixing, emplacement mechanisms etc.. Specifically, EBSD provides third dimension orientation information for crystals analysed within the 2D sample analysis plane (thin or thick petrologic section) through the relationships between crystallographic and crystalshape axes.This information is important for the study of meteorites derived from large achondritic bodies, such as the Martian nakhlites, where gravitational forces undoubtedly influenced crystallisation but the original orientation of the meteorites is unknown (Daly, Piazolo, et al., 2019).
Previous EBSD studies of a limited set of Martian nakhlite meteorites provide a baseline for understanding crystallisation processes in their origin on basis of detected SPO and CPO of earliest crystallisedand gravitationally settledmineral phases, particularly augite (Daly, Piazolo, et al., 2019).The data in this study, acquired on a much larger set of nakhlites (21 specimens (Fig. 4), confirm these observations; i.e., augite in each of the analysed materials shows a coupling of the {001} axis -long crystal-shape axis and {100} + {010} axes (i.e., <a> + <b> axes)short crystal-shape axes.
In general, clinopyroxenes (i.e., augite investigated in this study) will typically exhibit weaker fabric strengths than other associated minerals (e.g., olivine and orthopyroxene) as a function of its lower crystallographic symmetry and euhedral shape (Hidas et al., 2019;Van Der Werf et al., 2017).For augite within the nakhlites, the observed coupling of the crystal's euhedral shape to with the monoclinic crystal symmetry will contribute to the observed fabric strengths.Thus, when using quantitative metrics to assess fabric strength, the overall crystal's shape, crystal symmetry, and their relationship to each other needs to be considered to properly contextualise the results.
The M-Index values, which span a range of 0.03 (0.01 ± 0.0002 -0.04 ± 0.0007; Figs.7     and 11A), indicate no fabric within any of the analysed nakhlites.Note, however, that M-Index is commonly applied in the study of olivines and other orthorhombic minerals but not clinopyroxenes mainly because consideration of the latter alone can provide either erroneously low fabric strength estimates or even nonsensical values even in rocks with moderate fabric strengths (Hidas et al., 2019).The M-Index determines fabric strength through assessing the degree of rotation required for two crystals to become crystallographically aligned.To achieve this, the M-Index utilizes a subset of the overall orientation distribution function (ODF) spherical harmonic (Mainprice et al., 2015;Skemer et al., 2005).Among all the nakhlites studied here, a significant proportion of the augite crystals exhibit simple twins (180° rotations; Fig. 6).Note, however, that a current limitation of the used MTEX code is the inability to remove twinning orientation discrepancy within crystals exhibiting twinning of lower order symmetry crystals such as monoclinic augite.In other words, simple twins impact the M-Index calculations by underestimating true fabric strengths.The values that are observed do sit in a similar range to reported ultramafic rocks on Earth (Hidas et al., 2019;Van Der Werf et al., 2017), however in these studies clinopyroxene is not the dominant phase.For the analysed sections, the fabric manuscript submitted to Journal of Geophysical Research: Planets strength from M-Index calculations is at odds with the estimates from the CPO plots, but does agree with the observed SPO (Fig. 4).
The  2).These results, which are consistent with data from finegrained clinopyroxenites and coarser grained websterites on Earth (Frets et al., 2012;Henry et al., 2017), indicate significant textural diversity within the nakhlites emplacement environment.
Variation in textural fabric within a singular magmatic body (flow and/or intrusion) is common within igneous rocks, particularly those that formed in convecting flow regimes (Lofgren, 1983;Marsh, 2013;Perugini et al., 2003).
Girdle fabric CPO isdefined here as the alignment of one or two crystal axes along a singular plane, also known as magmatic foliation (Paterson et al., 1998).In addition to the more dominant {001} girdle CPO, almost all sections (exempting Caleta el Cobre 022 and Lafayette) exhibit very weak to weak intensity point fabrics within either the {100}, {010}, or {001} i.e., <a>, <b>, and <c> axes, respectively (Figs. 4,8,9 and 11), CPO point fabrics being defined as the consistent alignment of a given crystallographic axis along a singular plane also known as magmatic lineation (Paterson et al., 1998).Finally, some sections show PGR values within the {100} and {010} axes which indicate a combination of weak point with weaker {100} girdle fabrics (Figs. 8,9 and 11C).
The majority of the analysed samples exhibit a moderate {001} foliation (Fig. 4, 8, 9, and 11, Table 2) indicating a dominant shear or flattening flow regime, which is commonly observed within cumulate-type igneous rocks (Hunter, 1996;Iezzi & Ventura, 2002;Merle, 1998).Only Caleta el Cobre 022, Lafayette, MIL 090032, Y 000593 (section 127-A only), and Y 000802 exhibit weak {001} axis foliation (Fig. 4, Table 2), indicating a lower influence of pure shear (σ, i.e., where σ1>σ2=σ3).Out of the identified nakhlites with weaker strength foliations, Caleta el Cobre 022 is the only sample to exhibit no other form of textural fabric (foliation or lineation) within its {100} and {010} axes, which is consistent with its low J-Index value.Even Lafayette, the section with the next lowest foliation strength, exhibits a weak {010} girdle fabric that is weaker than the {001} axis girdle fabric (Fig. 4, Table 2).2).The general preference within the nakhlites towards the development of a weaker secondary fabric within the {100} axis rather than the {010} axis is inconsistent with most common CPO fabrics observed within clinopyroxene-rich igneous rocks on Earth despite plastic deformation within the crystal favouring dislocation along the {001} and {100} planes (Bascou et al., 2002(Bascou et al., , 2009;;Godard & van Roermund, 1995).Even though the values for lineation within Y 000593 are classified as weak, the alignment of the {010} axis over the {100} axis could indicate different external strain conditions during emplacement compared manuscript submitted to Journal of Geophysical Research: Planets to the other nakhlites (Bascou et al., 2009).For all but two of the studied nakhlites, the combination of stronger foliation with a weaker lineation suggests a common emplacement mechanism/environment, regardless of whether or not each sample represents a individual flow/intrusion (Corrigan et al., 2015;Daly, Piazolo, et al., 2019;Udry & Day, 2018).
Dominant pure shear (flattening) with a lower component of simple shear (rotational) is commonly associated with gravitationally driven crystal settling, the dominant emplacement mechanism in stagnant magma body regimes.Environments where such textures can develop include intrusions, lava pools and ponds (Iezzi & Ventura, 2002;Merle, 1998), low viscosity slow moving basaltic lavas, or as has been previously suggested for the nakhlites thick lava flows (Friedman Lentz et al., 2011;Treiman, 2005).However, variation (heterogeneity) in CPO fabric strength and the formation of weaker lineation textures is a common occurrence within magmatic bodies where pockets of texture can be highly localised, particularly for bodies associated with a higher discharge volume.In general, stronger fabric strengths which tend towards lineation are often observed along the edges of the magmatic body, gravitational settling fabrics towards the lower portions (but not the direct base), and random fabrics towards the central regions (Nicolas, 1992;Shelley, 1985).Previous CPO analyses of several nakhlites resolved microfabrics that are consistent with both subaerial flow and stagnant lava lake emplacement environments (Daly, Piazolo, et al., 2019).Presented observations among a much larger set of nakhlites are consistent with variable emplacement mechanisms and environments.Importantly, they also indicate that formation in pure crystal settling environments (e.g., stagnant lava lakes/ponds) may have been less common formation in dynamic flow regimes, be it intrusive or subaerial.In the case of an intrusive emplacement environment, such as a sill or a dyke, stronger CPO with a higher simple shear component is common along the margins but not the central portions of the magmatic flows (Gibb & Henderson, 1992;Komar, 1972;Piazolo et al., 2002;Shelley, 1985).In such a scenario, the here observed, moderate {001} foliations among augite crystals suggest crystal settling via dynamic sorting (Bagnold effect), where each of the nakhlites represents a different axial region from a series of interconnecting dykes or inclined sills (Corrigan et al., 2015;Marsh, 1996) which span at least four separate magmatic events (Cohen et al., 2017).The identification of secondary multiple axis lineations plus the random SPO in most of the analysed samples suggest that the observed {001} foliations have overprinted the initial primary flow textures within each of the putative nakhlite igneous bodies.During emplacement, such a switch from flow to gravitational crystal settling would account for random SPO, conform with the textural interpretations of Udry and Day (2018) whilst also maintaining regions of stronger texture particularly along the margins of intrusion (for comparison see Berkley et al., 1980).
To assess variability in individual analyses for each meteorite two replicate sections from five different meteorites (Governador Valadares, Nakhla, NWA 998, Y 000593, and Y 000749) were analysed.Previous work assessing modal mineralogy variation between replicate nakhlite sections (i.e., sections sourced from the same meteorite) showed up to 40% variation in the nakhlites (Corrigan et al., 2015).Significant changes in modal composition will impact the expression of texture, particularly within igneous rocks.Higher levels of mesostasis material will, in general, enable easier mobility of phenocryst crystals, whilst highly crystalline solids will have less movement (Nicolas, 1992;Piazolo et al., 2002).The quantitative CPO metrics calculated standard deviation between replicate sections were highest for J-Index (0.75 ± 0.13), followed by PGR (0.022 ± 0.002), then M-Index (0.002 ± 0.001; Table 2).The M-Index and J-Index standard deviations, sit within the same fabric strength groupings (weak, moderate etc.)  2).Both M-Index and J-Index calculations are highly dependent on the orientation distribution function (ODF).A sample's ODF represents the variability in crystal orientation relative to the analysis plane.Crystal orientation used in ODF is calculated using the Euler density, which is derived from the crystal's crystal symmetry.Calculation of an accurate ODF is highly dependent on the associated halfwidth (selected as 10° for all samples in this study), which if not properly considered can distort the Euler space and bias the results (Gerth & Schwarzer, 1993;Kalidindi et al., 2009).Variable sized datasets and failure to properly accommodate fractured crystals can further influence CPO calculations.However, the third CPO metric, eigenvalue analysis (PGR), does not rely on ODF.What the PGR metric does instead is utilise the direct Euler angle of the crystal with respect to the crystal's Miller reference frame (which for augite is {100}, {010}, {001}).The orientation deviation of each crystal relative to the overall sample is then normalized.This normalised orientation is then assessed as a fractional component (i.e., using one Miller reference frame) between three endmembers (P, G, and R).This is distinctly different to the 2 end-member scale used for both M-Index and J-Index, which both consider a crystal in its' entirety.
Across the replicate sections the major variable that could have affected the resulting CPO values is the discrepancy in the number of detected crystals.Assessment of the compositional breakdown between phenocryst, mesostasis material, and melt/void space within the samples (Fig. 2) showed the largest deviations between all three components in Y 000593 importantly, it provides an increased understanding on how the CPO metric is influenced by the inherent heterogeneity of the sample.Variability in the number of augite crystals for CPO calculations was minimised through setting an upper limit of ≤300 crystals for the larger datasets (affecting 12 datasets including all but Y 000593 in the replicate sections).Despite 100-150 crystals being shown to produce statistically stable results for fabric calculations (Ismaïl & Mainprice, 1998;Skemer et al., 2005), ≥300-≤600 crystals over several different sections is currently considered the gold standard for assessing stable textural results within terrestrial geological samples.This gold standard can be impractical and potentially problematic for many planetary studies, where often the material available to study is limited and the amount of sample per section is generally a fraction of the size typically found in a traditional petrological thin section.If this ideal range were to be applied to the presented nakhlite datasets, only 12 of the 21 sets using nMTEX (4 datasets using nOPPG; Table 2) would be considered viable for assessing texture.Thus, textural CPO metrics should only be used crude indicators for texture particularly in the study of scare meteorites.
manuscript submitted to Journal of Geophysical Research: Planets 4.2 Textural implications for the nakhlite source volcano and magmatic edifice Igneous rocks, whether formed on Earth or elsewhere in the Solar System, usually show significant mineralogical heterogeneities down to the scale of a few tens of micrometers (Hammer et al., 2010;Jankovics et al., 2012;Kouchi et al., 1986;Lofgren, 1983;Marsh, 2013;Perugini et al., 2003).This characteristic hampers the extrapolation across scales, particularly for meteoritic materials like the nakhlites, of which only millimetre to centimetre-sized fragments are available for in-depth microtextural analysis.To date, however, significant textural variation within an individual Martian meteorite specimen (303 samples as of 2021) has only been reported for the shergottite Zagami (Becker, 2011).Different conditions on Mars e.g., lower gravity, would undoubtedly affect the textural and petrographic expression of Martian magmatic rocks when compared to those formed on Earth.For example, a relatively stronger strain field is required to produce lineation textures within equal-sized magmatic rocks on Mars compared to Earth, mainly because of the lower gravity on Mars and its subsequent effects on parameters like viscosity, strain, apparent density, and general flow dynamics (Niu & Pang, 2020;Vetere et al., 2019).However, despite the differences between Earth and Mars creating variation in the parameters required to form particular textures, the texture's significance in terms of the general environment of formation are directly comparable.For igneous rocks, increased M-Index and J-Index metric values signal a flow-dominated environment, accounting for increased lineation/simple shear; in PGR analysis, the same is indicated by higher P values.Crystal alignment, which influences texture development, has been shown to have a number of controls outwith the external strain field (simple shear).The alignment of crystals within an igneous body will be influenced by factors such as the elongate nature of the individual crystals, density and viscosity of the melt, the mineralogical composition and distribution, as well as the surface: volume ratio of the igneous body (Piazolo et al., 2002).Textural heterogeneity between different igneous specimens may be even more pronounced in cases that they formed in multiple events of emplacement (Fenton et al., 2004).Despite the lack of (geological emplacment) context for the nakhlites, the size of each meteorite constrains the maximum relative distance between the replicate sections to centimetre-scale, and the vicinity of the individual meteorites relative to one another to within the area of the currently unknown ejection crater, whose diameter is estimated to be on the scale of kilometres (Artemieva & Ivanov, 2004;Kereszturi & Chatzitheodoridis, 2016).
The bulk volume of an igneous body can have a significant impact on the external strain field's textural expression within an igneous sample (Marsh, 2013).Specifically, a specimen that forms in a large crystal convecting flow/intrusion (which produces localised strain fields with at least two major strains directions within a larger strain field) would be expected to exhibit a stronger textural fabric.The strength of this textural fabric would also vary more throughout the igneous body than a sample formed from either crystal settling (one dominant strain field a.k.a.gravity) or thermal diffusion [no strain field; Callot & Guichet (2003), Marsh (2013), and Nicolas (1992)].However, the model calculation of equation 5, which was intended to provide the absolute minimum potential magma body thickness, results in unit thicknesses for each specimen that are significantly below the average crystal size (Tables 2 and 3).Thus equation 5's results do not support a mechanism of emplacement that was purely controlled by crystal convection.However, crystal convection would be expected for magma body unit thicknesses calculated to be on the 10's m scale in order to satisfy fundamental physics regardless of the main mechanism of emplacement.If crystal convection was the overriding emplacement manuscript submitted to Journal of Geophysical Research: Planets mechanism within the nakhlites, then based off the calculated magma body thicknesses the residence times used in our calculations would have to represent a single round of convection rather than representing the entire crystallisation time period.Out of the three emplacement endmembers modelled for the nakhlites, crystal settling agrees best with the observed textural data (Figs. 4 and 8) and moreover, produces the largest range of magma body thicknesses (<1 m to several 10's m) which can be distinguished as two separate groups (Figs.10-11, Table 3).
However, because geological samples rarely represent pure endmembers components both crystal convection and/or thermal diffusion cannot be discounted for the nakhlites, particularly in identified low fabric strength samples, such as Caleta el Cobre 022.
Coupled modelling and CPO textural analysis results demonstrate that the nakhlite specimens analysed formed via gravity settling dominated crystallisation with variable minor components of flow.These results are in agreement with, and expand upon the findings in Daly, Piazolo, et al. (2019).Whereas these results are in striking contrast to the findings in Udry and Day (2018), who reported random textures for the thirteen nakhlite specimens they studied.The conflicting results in Udry and Day (2018) may be due to their use of more traditional SPO measurements corrected for slope via measured crystal shape dimensions (rather than proxied SPO via CPO used here), or from the nakhlites forming under a dynamic sorting regime (Corrigan et al., 2015).An intrusive dyke field could produce the observed random SPO with an aligned CPO via a dynamic sorting regime (Horsman et al., 2005), as reported in Figure 4.In all studied specimens, the percentage of crystallographically aligned phenocrysts within the analysis plane is low (2.53-10.04%;Table 2).This results in the observable SPO, even when accounting for the phenocryst's slope relative to the analysis plane, to appear random in almost all analysed nakhlite sections (Fig. 4).If each section's CPO was not considered, then the presented data from this study would agree with the interpretations of Udry and Day (2018).
Lastly, it is important to evaluate whether the calculated magma body thicknesses for each specimen can aid in localisation of the nakhlites' origin in the Martian crust.Magma body thicknesses modelled using crystal settling for the individual nakhlites show significant variation in thickness from <1 to several 10's m, enabling (with the models current resolution) two distinct groups to be identified (group 1 being <10 m and group 2 being >10 m; Table 2).These two groups show no distinct trends or correlations to textural data, published geochemical data [e.g., Day et al. (2018) andKrämer Ruggiu et al. (2020)], or age [e.g., (Cartwright et al., 2013;Cohen et al., 2017;Krämer Ruggiu et al., 2020;Mikouchi et al., 2016;Park et al., 2009Park et al., , 2016))] inferring the overall dataset to be inherently random in terms of both textural fabric variation and magmatic body thicknesses through time (Figs 10,11).This heterogeneity of the larger dataset could be very easily overlooked when only assessing a subset of the samples.On Earth, individual magmatic events from the same volcano or even multiple lobed flows from a singular event have been observed to shift and evolve in an inconsistent manner, particularly when analysing relatively short periods of a volcano/magma bodies lifespan (Gamble et al., 1999).
This intrinsic variability is why field-mapping and high-resolution age dating are critical tools used to understand the evolution of volcanic systems over shorter time periods.Even though the nakhlites are shown to have a formation window of at least 93 ± 12 Ma (Cohen et al., 2017) magmatic sources on Mars have been shown to span several billion years (Lapen et al., 2017).In terms of the positioning of the nakhlites, the distial relationship between the samples is limited to their ejection crater on Mars (Cohen et al., 2017;Day et al., 2018;Nyquist et al., 2001;Udry & Day, 2018).Comparing the calculated magma body thicknesses to published age data showed manuscript submitted to Journal of Geophysical Research: Planets that thicknesses varied over time in an inconsistent manner, which could indicate variable discharge from the nakhlite magmatic source, or be the result of other factors such as magma viscosity, topography, restriction/spread of the flows/intrusions etc. (Baumgartner et al., 2017).
The resolution of the here presented magma body thicknesses do not yet enable localisation parameters to be fully discerned at a reasonable resolution.However, they do provide a starting point for building more complex models to start investigating other influential parameters such as, mixing of the here calculated endmember emplacement mechanisms, intrusion versus flow, depth of burial, cooling rate, distance travelled from source etc., which in lieu of field-mapping should assist in narrowing the localised criteria for the nakhlites.

Conclusions
Nakhlite CPO indicates the presence of a low intensity weak-medium foliation fabric within the {001} (<c>) axis.In addition to the identified foliation fabric all analysed specimens, all but two of the analysed specimens, also exhibit a weak form of lineation fabric within either one, two, or all three crystallographic axes.The dominant foliation fabric is common for many cumulate terrestrial rocks, and agrees most with an emplacement mechanism based on crystal settling such as lava ponds, lava lakes, intrusive sills, and dykes.Disparity between observed SPO and CPO across the nakhlites could either be the result of the sections' orientation or an overprinting of initial fabrics by the now more dominant {001} foliation.Calculating end member magma body thicknesses showed crystal settling to be the mechanism that best agrees with reported textural data.From the crystal settling results two distinct groups (outwith error) were identified for the nakhlites; those with calculated thicknesses <10 m and those with calculated thicknesses >10 m.The relationship of these two groups to the presented CPO, SPO, as well as published age and geochemistry show no correlations or trends highlighting the inherent randomness of the large dataset despite all analysed samples exhibiting a common dominant emplacement mechanism and similar CPO.Overall, this study highlights the importance of assessing larger datasets, particularly for samples (such as meteorites) where due to the lack of external context, the assessment of limited data can lead to the identification of non-existing trends and relationships.

Figure 1 .
Figure 1.Examples of SPO bias and resective CPO relative to the section analysis surface orientation [red (perpendicular), yellow, (oblique) and green (parallel)] with respect to the fabric orientation for a crystal settling regime producing girdle fabrics (i.e., foliation).

Figure 2 .
Figure 2. Nakhlite EBSD phase maps showing the typical mineralogy of the nakhlite meteorites, consisting of augite (high Ca-clinopyroxene; green) and lesser olivine (blue; indexed as forsterite) together with minor orthopyroxene (indexed as enstatite), iron-titanium oxide (indexed as titanomagnetite), and plagioclase feldspar (indexed as albite).Note the high proportion of plagioclase mesostasis material in the Caleta el Cobre 022 meteorite compared to the other nakhlites studied here.3.1 Modal mineralogy and augite morphology

Figure 3 .
Figure 3. Modal mineralogy of the nakhlites from EBSD data.A) Volume percent of phenocrysts (>0.2 mm crystals), mesostasis material (<0.2 mm crystals), and glass-non indexed material.B) volume percent of indexed phases in EBSD data i.e., this is for both phenocryst and mesostasis phases.

Figure 7 .Figure 8 .Figure 9 .
Figure 7. M-Index and J-Index values for the nakhlites.Caleta el Cobre 022 exhibits the lowest J-Index and M-Index values, while the highest observed J-Index and M-Index values for NWA 817 are considered to be a function of the smaller dataset size (nMTEX = 43 crystals).
J-Index values indicates different fabric strengths across the across the nakhlite suite, ranging from borderline no fabric to weak fabric in Caleta el Cobre 022, to medium strength fabrics in Governador Valadares, MIL 090136, NWA 998, NWA 11013, Y 000593, Y 000749, and Y000802 (Figs. 7 and 11A, Table manuscript submitted to Journal of Geophysical Research: Planets but are often outwith the metrics associated calculation error.The PGR values sit qualitatively within the same regions (Fig 8), however, quantitatively they do not overlap (Table

[
±21 (c), ±16 (mes), ±5 (ni)] and the smallest deviation within in Y 000749 [±2 (c), ±1.5 (mes),±3.5(ni)].In terms of the vol.% augite within the samples, deviation within the replicate samples is observed to range from ±0.7 (Y 000593) to ±19.3 (Governador Valadares).Placing the replicate section CPO metric deviations within the context of augite's variable abundance within each section helps to explain the observed differences in metric results outside of analytical uncertainty.It also reinforces the importance of crystal distribution heterogeneity to the development and variation of observable textural fabric within a given sample.More with corrections from Mainprice manuscript submitted to Journal of Geophysical Research: Planets et al. (2015), generating J-index, M-index, and Eigenvalue values

Table 2
manuscript submitted to Journal of Geophysical Research: Planets