Compositional Enhancement of Crustal Magnetization on Mars

Martian orbital and lander measurements revealed strong (∼1–2 orders of magnitude greater than Earth) crustal magnetic anomalies and the lack of an active detectable core dynamo. This strong crustal magnetization remains unexplained given that models of an ancient core dynamo on Mars predict surface field strengths comparable to modern Earth. We explored the relationship between Mars' crustal magnetization and its composition in multivariate space. We identified that 530 and 1,000 nm absorptions (from orbital spectrometers) have unique correlations with crustal magnetization in the Terra Sirenum‐Terra Cimmeria region and ∼13% of the variance of the magnetization can be attributed to these correlations. Because we do not expect the topmost material, detectable by the orbital spectrometers, to retain magnetization from an ancient core dynamo, we propose this material is compositionally similar to the bulk rock below it, which is more likely to retain magnetization. Therefore, the observed variance is a lower limit.

Based on the measured magnetic field from orbit, at least 600 kA must be accounted for through the product of crustal volume and its vertically integrated magnetization (Dunlop & Arkani-Hamed, 2005). The areal extent of the magnetized region is known with the approximate precision of the satellite altitude (∼100-200 km). The volume of the magnetized crust can be estimated by placing an upper-bound on the thickness of the magnetized layer based on the calculated depth to the Curie point isotherms (Connerney et al., 2001;Dunlop & Arkani-Hamed, 2005). Given the candidate magnetic phases (magnetite, hematite, maghemite, pyrrhotite, titanomagnetite, titanohematite, and, least likely, kamacite), the amount of magnetic material required ranges from 0.4 to 6 wt.% if the entire layer is magnetized (Dunlop & Arkani-Hamed, 2005). These proportions are consistent with measurements made on the Martian surface (e.g., Vaniman et al., 2014) and of some Martian meteorites (Gattacceca et al., 2014).
Recently, the magnetometer onboard InSight measured a surface magnetic field nearly 10X stronger than modeled by previous satellite-based observations ( Figure 1; Johnson et al., 2020;Langlais et al., 2019;Morschhauser et al., 2014). This observation is not entirely unexpected as smaller (below orbital resolution), local variations are expected. Yet, the potential contribution from local variations does not fully explain the strength and cohesiveness of magnetic fields detected above Mars in the TSTC region.

Data Sets
Predictions of surface magnetic field assessed in this study are derived from a model of the crustal magnetic field of Mars which combines measurements from the Mars Global Surveyor (MGS) and Mars Atmosphere and Volatile EvolutioN (MAVEN) magnetometer (Langlais et al., 2019). The data sets used for comparison included elemental maps from the Gamma Ray Spectrometer (GRS) (Boynton et al., 2004), mineralogical parameters from Compact Reconnaissance Imaging Spectrometer for Mars (CRISM), Observatoire pour la Mineralogie, l'Eau, les Glaces et l'Activité (OMEGA), and the Thermal Emission Spectrometer (TES), as well as TES-derived surface types. A summary of the data sets used, their resolution, and coverage are provided in Table 1. ALHANTOOBI ET AL.

Statistical Analysis
There are many covarying relationships between crustal properties. Considering two variables at a time in bivariate analysis can create misleading apparent correlations (Karunatillake et al., 2012). To accurately compare these disparate data sets, we applied a multivariate correlation and regression methodology similar to that of Karunatillake et al. (2012) with magnetic field model result as the response variable and compositional data sets as the predictor variables. We used two distinct yet complementary multivariate analysis methods of Karunatillake et al. (2012)-ordinary linear regression (OLR) and spatially weighted linear regression (SWLR)-to estimate the correlation between the magnetic field model and our predictor variables (all other data sets). OLR estimates the partial correlation of the magnetic field with each predictor variable. SWLR modifies the OLR method to account for spatial autocorrelation. The multivariate analysis methods help to verify the significance of observed correlations in the bivariate analysis and can reveal correlations that were not initially apparent.
We calculated the Langlais crustal magnetic field model (Langlais et al., 2019) at a resolution matching that of TES. We then trimmed the geographic extent of each data set to match the most limited footprint and sampled each data set to match the resolution of the lowest-resolution data set. Finally, we excluded regions where the predicted surface magnetic field was below 100 nT ( Figure 1) and where TES lambert albedo was above 0.2 (to avoid extensive dust cover).
Full equations and methods are available in Karunatillake et al. (2012) and references therein with differences in Section S1 (supporting information). The variance contribution, R 2 , is a numeric value between 0 and 1 that represents the fraction of the response variance attributable to the predictors in each method. The partial correlation coefficient, r i , indicates the unique correlation between the response variable and the individual data set. To verify if the methods are sufficiently robust for the data, we use the expected R 2 if the response were uncorrelated with any of the predictors (R 0 2 ) and R 2 adjusted for the degrees of freedom (R df 2 ). If R 0 2 exceeds 0.1 or R df 2 is less than 0.1, then any apparent correlation can be considered insignificant. To assess the significance of the correlation estimates, we used the p-value from the Student's t-test distribution for OLR and standard normal probability distribution for SWLR to eliminate predictors that fail to meet the 95% confidence threshold. Univariate (OLR and SWLR) and Bonferroni (OLR) confidence intervals were computed. Predictor variables with confidence intervals that include zero (indicating very minimal/no significant correlation between the response and the predictor) were eliminated. Any analysis that revealed key predictor variables through OLR was then subjected to SWLR. Predictor variables with contributions >0.05 were considered key predictors.
SWLR analysis requires quantifying the degree to which individual bins are distinct from neighboring bins. Additionally, spatial weighting due to map projections and sampling differences at the poles and equator is incorporated. To define a spatial weighting matrix, we plotted the similarity function for all our data sets. This function allowed us to determine the similarity between values of bin pairs at given angular separation. Beyond a threshold angular separation, the similarity value remains constant. This threshold angular separation can then be used to find the constant of decay for bin similarity. We obtained the constant of decay for the highest threshold angular separation to ensure that spatial autocorrelation has been accounted for across all predictors. The "Model σ" is the uncertainty of the overall spatial autocorrelation fit, "ρ" is a scalar of the spatial weights matrix determined by the similarity function and "C" is the constant of regression.
We first attempted multivariate analysis at the lowest resolution data set (GRS), which enabled the inclusion of all data sets. The geographic extent and large footprint of GRS resulted in too few datapoints for analysis. These results are further discussed in Section S2 (supporting information). OLR and SWLR were then applied at the second lowest resolution (TES) with the goal of reducing the number of predictors and simultaneously avoiding degradation of the potential variance contribution as a result of binning. Analyses were conducted both on global and regional scales.
For regional analyses, we used the compositional provinces identified by Gasnault et al. (2010) and Rogers and Hamilton (2015) and geologic units mapped by Tanaka et al. (2014). Individual geologic units are quite small for robust statistics, so geologic units of similar age were also grouped and analyzed. Results of these analyses did not yield meaningful correlations and are discussed in Section S2 (supporting information).

Analysis at TSTC
We conducted additional analyses on the region of maximum crustal magnetization, TSCS. The maximum field in this region is > 11 μT at the surface as predicted by the global model (Langlais et al., 2019). To define the exact region of interest, we generated scatter plots between the predicted magnetic field and the ferric mineral indicators. Those that show a positively correlated trend within the global distribution were chosen ( Figure S1), mainly the 530 nm band depth (BD530) and nanophase iron spectral parameter (NNPH). The data sets were normalized and each point on the graph was assigned a slope value relative to the origin. We then selected a range of slope values and captured all the points that lie within the range. The slopes selected for BD530 ranged between 0.4 and 1.2 and those of NNPH between 0.6 and 1.6. Furthermore, we excluded parameter values that were within 10% of the minimum value and magnetic field values that were below 400 nT (to exclude transitional areas between positively and negatively magnetized regions). The remaining points from both scatter plots were then combined into one mask and points outside of TSTC were discarded (Figure 2).

Results
A scatter plot matrix of points within the region of interest (n = 19,133) was generated ( Figure S1). This initial bivariate analysis reveals that the 530 nm band depth (BD530), NNPH, lambert albedo, and thermal inertia (TI) all show potential positive correlation while high calcium pyroxene (HCP), 2,000-nm band depth (BD2000), 1,000-nm band depth (BDI1000), and surface type 1 abundance all show potential negative correlation (see Table 2). Importantly, the bivariate analysis results are not considered of significance until they are verified using the multivariate OLR and SWLR methods.
10.1029/2020GL090379 5 of 10  Using OLR analysis, many predictor variables indicate unique correlations with magnetic field (Figure 3a). Several predictor variables have confidence intervals that do not include zero and therefore pass both the univariate and Bonferroni confidence intervals (Figure 3b). The variance contribution from this analysis indicates that ∼17% of change in the magnetic field can be linked to changes in the predictor variables.  The maximum threshold angular separation, after which there is insignificant contribution from neighboring bins, is 72°. Using these threshold angular separations, we were able to determine the spatial weighting matrix for SWLR analysis. As expected, fewer predictor variables pass this analysis with significance (see Figure 3c, Table 2). Of the predictor variables that pass, only two do so with significant contribution (>0.05; Figure 3d). The variance contribution from this analysis indicates that ∼13% of the change in magnetic field can be linked to changes within these variables, primarily from BD530 and BDI1000VIS.

Meaning of Key Predictors
Although several predictor variables pass our SWLR analysis, only BD530 and BDI1000VIS contribute more than 0.05% to the variance contribution. Both indicators are related to iron bearing mineral phases (Ody et al., 2012;Viviano-Beck et al., 2014) and therefore potentially sensitive to magnetic phases. BD530, associated with absorptions at 530-nm, has been attributed to multiple spectroscopic effects from Fe 3+ while BDI1000VIS, a broad absorption at 1-μm, is typically associated with electronic transitions and charge transfer absorptions from iron in the crystal structure of olivine and pyroxene or from Fe-bearing glass. Vis-NIR reflectance for four candidate magnetic phases from the USGS spectral library (Clark et al., 2007) and the calculation of their parameters is shown in Figure S2. Our statistical analysis does not indicate which magnetic phases are most likely.
Different iron-rich minerals have distinct susceptibilities to remanent magnetization. Magnetite, pyrrhotite, and hematite have all been proposed as magnetic carriers on Mars (Dunlop & Arkani-Hamed, 2005;Gattacceca et al., 2014;Kletetschka et al., 2000;McEnroe et al., 2004a;McEnroe et al., 2004b;Nimmo, 2000). Previous studies assumed that the magnetized layer in TSTC began at ∼10 km depth and continued downwards until the Curie point(s) of the magnetic mineral(s) were reached (Dunlop & Arkani-Hamed, 2005;McEnroe et al., 2018). For example, ∼0.8 wt.% of single-domain magnetite in a ∼30-40-km thick layer could produce the necessary vertically integrated magnetization of ∼600 kA. However, recent results from the NASA In-Sight Mission  and new MAVEN data  show that some magnetized layers on Mars could be thin and/or shallowly buried (∼200 m-10 km). Figure S5 illustrates that thin layers of magnetized crust (<2-km thick) can produce the predicted magnetic fields at TSTC if roughly one third of the iron (>10 wt.% from GRS data) is in the form of single-domain magnetite and/or pyrrhotite (using an average density for the bulk crust of 2,900 kg/m 3 ; Zuber, 2001). Other minerals (hematite, maghemite, and/or multidomain magnetite) are only credible candidates for the primary magnetic carrier(s) if the magnetized layer is >10-20-km thick (Section S3, Table S1, supporting information).

Local Compositional Enhancement of Magnetization
The variance contribution from BD530 and BD1000VIS indicates that composition of the crust influences its magnetization. Using the data sets available for this study, we can account for ∼13% of the variance in magnetic field. It is important to consider the sensitivity of the instruments we used for comparison with magnetic field are only sensitive to the top ∼50 μm (TES) or top ∼ few μm (CRISM and OMEGA). This uppermost material has likely undergone some combination of erosion, alteration, impact cratering, or other process which removed any magnetization from an ancient dynamo. Instead, the magnetization detected from orbit likely originates in either the bulk rock or within deeper layers of similar composition, which are more likely to preserve remanence from an ancient dynamo. This interpretation is consistent with the measurements observed at the InSight landing site  and Lucus Planum  in which the magnetization is carried in comparatively shallow layers, contrary to previous magnetic source depth models (Dunlop & Arkani-Hamed, 2005;Lewis & Simons, 2012). We have applied our statistical analysis at Lucus Planum and found that at least 7% of the variance contribution to predicted magnetic field can be explained through mineralogy ( Figure S3).
Geological mapping in the TSTC region (Tanaka et al., 2014) indicates that the majority of the exposed material is a Noachian highland unit, which in our analysis is exposed within a ∼6-km deposit ( Figure S4). In order to explain the correlation between surface composition and bulk or buried magnetization, there must be compositional similarity between the uppermost material in TSTC and the layers that retain remanent magnetization. Our preferred explanation is that the uppermost material is derived from the bulk rock, which carries the magnetization. This assumption is further validated by the recent crustal magnetism analysis at Lucus Planum which has shown magnetization carried within an upper unit . Using estimates for the lifetime of the ancient dynamo (Moore & Bloxham, 2017;Roberts et al., 2009), we can expect that the majority of the remanence is carried in the early and middle Noachian units mapped by Tanaka et al. (2014) in TSTC. Other possible explanations are iron-enrichments resulting from hydrothermal alteration or dike emplacement. However, these later deposits/intrusions must be large enough to produce an observable signal in the remote sensing data taken of the surface, and therefore we consider them less likely sources of the predicted magnetic field. In all scenarios, the observation of mineralogy correlated with magnetization remains.
The variance contribution is likely an underestimate of the compositional enhancement in TSTC. Any portion of our region of interest which lost its magnetic remanence, such as the uppermost material discussed above, would have the effect of reducing the variance contribution we observe because it would be uncorrelated but interspersed within the correlated data. Therefore, the positive correlations we observe indicate that without this compositional influence the intensity of the crustal magnetization would be reduced. As previously discussed, the magnitude of Mars' crustal magnetization has previously conflicted with current models for the strength of potential ancient core dynamos (e.g., Stevenson, 2001). This result reduces the discrepancy between models of ancient field strength and crustal magnetization.
In addition to the compositional enhancements observed through this analysis, there are broad implications for the history of the TSTC region. The correlation we observe between composition and magnetization suggests that the surface material in TSTC is locally sourced and has not mixed significantly with global regolith. Tanaka et al. (2014) indicate that this region primarily consists of Noachian-aged units ranging in thickness between hundreds of meters to a few kilometers. Therefore, TSTC as a whole has remained largely unaffected by demagnetizing processes since the Noachian.

Conclusion
Our analysis indicates that in the TSTC region on Mars, the strong magnetization is at least partially due to compositional (iron) enhancements and not necessarily solely from a strong paleofield. This finding significantly reduces the discrepancy between models of ancient field strength and crustal magnetization and suggests that Mars' ancient paleofield does not need to be of exceptionally high strength. This work places additional constraints on the depth and thickness of the magnetized layer and/or the amount of magnetized minerals and is consistent with previous findings related to Mars crustal magnetic mineralogy (Dunlop & Arkani-Hamed, 2005; e.g., Kletetschka et al., 2000). For example, two possible scenarios are (1) <2-km thick layer containing minerals with high saturation magnetization (such as single domain magnetite and/or pyrrhotite), or (2) a > 10-20-km thick layer with minerals of weaker saturation magnetization (such as hematite and/or maghemite). When considered with other recent results (e.g., Johnson et al., 2020;Mittelholz et al., 2020), much of Mars' crustal magnetization can be explained by shallow and/or thin layers of magnetized material, dramatically different from early models (e.g., Connerney et al., 1999). Finally, heating or other demagnetizing processes must be absent at large scales in TSTC for the magnetization to be coherent (consistent with Bouley et al. (2020)). The uppermost material (<10 μm) in TSTC is compositionally similar to the bulk rock, which retains the magnetization, indicating it is recently eroded or locally sourced.

Data Availability Statement
All data sets used in this work are available through JMARS (https://jmars.mars.asu.edu/) or the NASA Planetary Data System (PDS) (https://pds-geosciences.wustl.edu/dataserv/mars.html). The Mars Odyssey Gamma Ray Spectrometer (GRS) derived chemical maps were derived from the spectral data archived on the PDS by applying the spectra-to-chemistry modeling methods described by Boynton et al. (2007), Evans et al. (2006), and Karunatillake et al. (2007). The derived chemical maps have also been described and ana-