New and Extended Data Processing of Mars Odyssey Neutron Spectrometer Data

TL;DR

This paper introduces a new extended analysis of Mars Odyssey Neutron Spectrometer data, providing updated neutron maps and preliminary insights into inter-annual polar variability on Mars.

Contribution

It presents a novel extended dataset and analysis methodology for MONS data, improving temporal coverage and enabling new studies of Martian polar regions.

Findings

Extended MONS neutron data maps created

Preliminary inter-annual variability results in polar regions

Comparison with previous MONS data analyses

Abstract

The Los Alamos National Laboratory designed and built Mars Odyssey Neutron Spectrometer (MONS) has been in excellent health operating from February 2002 to the present. MONS measures the neutron leakage albedo from galactic cosmic ray bombardment of Mars. These signals can indicate the presence of near-surface water deposits on Mars, and can also be used to study properties of the seasonal polar CO$_2$ ice caps. This work outlines a new analysis of the MONS data that results in new and extended time-series maps of MONS thermal and epithermal neutron data. The new data are compared to previous publications on the MONS instrument. We then present preliminary results studying the inter-annual variability in the polar regions of Mars based on 8 Mars-Years of MONS data from the new dataset.

Figures

Tables

Table 1. Table 2: ADC Nonlinearty correction factor, repeated for each 16 channel group.

Channel	Correction Factor
0x	1.00864
1x	0.97306
2x	1.04867
3x	1.03395
4x	0.95991
5x	1.05862
6x	0.92326
7x	1.05427
8x	0.97688
9x	1.02267
10x	0.99594
11x	0.99265
12x	0.89037
13x	1.14551
14x	1.00161
15x	0.91398

Table 2. Table 3: Elemental weight fractions in the base soil composition based on S21 from [ 20 ]

O	0.4223		K	0.0022
Na	0.0147		Ca	0.0419
Mg	0.0478		Ti	0.0052
Al	0.0464		Cr	0.0016
Si	0.2075		Mn	0.0033
S	0.0278		Fe	0.1466

Table 3. Table 4: Fit parameters (squared term, linear term, and constant term) for determining atmospheric scaling factor relative to 16 g/cm 2 .

Prism/Term	1% WEH	3% WEH	4% WEH	5% WEH	6% WEH	8% WEH
Prism 1 Sqr	3.169e-4	3.133e-4	3.752e-4	4.160e-4	4.265e-4	4.237e-4
Prism 1 Lin	-2.227e-2	-1.990e-2	-2.077e-2	-2.143e-2	-2.060e-2	-1.880e-2
Prism 1 Cns	1.275	1.241	1.239	1.235	1.221	1.192
Prism 2 Sqr	2.831e-4	1.827e-4	2.405e-4	2.956e-4	2.545e-4	2.088e-4
Prism 2 Lin	-1.145e-2	-1.114e-2	-1.393e-2	-1.608e-2	-1.485e-2	-1.418e-2
Prism 2 Cns	1.110	1.137	1.165	1.179	1.174	1.177
Prism 3 Sqr	3.148e-4	3.132e-4	3.751e04	4.158e-4	4.289e-4	4.258e-4
Prism 3 Lin	-2.216e-2	-1.981e-2	-2.067e-2	-2.132e-2	-2.056e-2	-1.873e-2
Prism 3 Cts	1.274	1.240	1.238	1.234	1.219	1.190
Prism 4 Sqr	2.753e-4	3.189e-4	3.701e-4	4.173e-4	4.571e-4	4.629e-4
Prism 4 Lin	-2.073e-2	-2.097e-2	-2.183e-2	-2.298e-2	-2.321e-2	-2.157e-2
Prism 4 Cts	1.261	1.257	1.257	1.259	1.255	1.226

Equations

$$\Omega (h)=2\pi \left[1-\sqrt{1-R^2/{\left(R+h\right)}^2}\right].$$

$$\lambda =\text{atan2}\left(\sin \left({\phi }_2-{\phi }_1\right)\cos {\theta }_2,\cos {\theta }_1\sin {\theta }_2-\sin {\theta }_1\cos {\theta }_2\cos \left({\phi }_2-{\phi }_1\right)\right).$$

$${\theta }_{\text{new}}^i$$

$${\phi }_{\text{new}}^i$$

$$\Delta t_r=\frac{R_M(m/2V{)}^{1/2}}{2{\left[1-K/V\right]}^{3/2}}\times$$

$$\left\{2\mu {\left(1-\frac{K}{V}\right)}^{1/2}{\left(\frac{K}{V}\right)}^{1/2}\left[1-{\left(\frac{{\tan }^2\theta }{{\tan }^2{\theta }_R}\right)}^{1/2}\right]+{\sin }^{-1}\left(\frac{B}{[A^2+B^2{]}^{1/2}}\right)+{\sin }^{-1}\left(\frac{1-2K_R/V_R}{[A^2+B^2{]}^{1/2}}\right)\right\},$$

Full text

1 Introduction

The Los Alamos National Laboratory (LANL) designed and built Mars Odyssey Neutron Spectrometer (MONS) has been continuously measuring the leakage flux of neutrons from a Mars polar orbit since February 2002. These data have been used to map the hydrogen content of Mars up through July 2009 [1, 2]. However, both the Odyssey spacecraft and the MONS instrument have been operational to the present time; here we present a new and extended analysis of the MONS dataset through December 2017 that can be used to search for long-term climate variations, particularly in the polar regions of Mars. Processing of the integrated total dataset is an important step before such interpretations can be drawn. Throughout this process, many choices and alternatives must be clearly documented to establish the accuracy, precision, and robustness of these data.

The MONS instrument collects neutron fluxes continuously from the Mars surface in three energy bands: thermal (0–0.4 eV), epithermal (0.4 eV–700 keV), and fast (0.7–5 MeV). The purpose of this work is to transform time-tagged measurements through the present time into relevant neutron maps, which are time dependent because of seasonally changing CO2- and water-ice precipitation at high latitudes. MONS has continued to collect data over 17 years of operation to present, well over its initial science phase duration. However, over that time period, there have been several unresolved issues regarding our understanding of the MONS systematic biases. At present, most, if not all, of these biases have been removed. Each generation of MONS data processing was built independently of its predecessors, often by different people, to limit the propagation of erroneous assumptions. Efforts were also devoted to compare the results of each approach.

The initial processing of MONS data was performed by Tokar et al. [3] and was used for early discovery results. Subsequently, Prettyman et al. developed an independent approach [4] that has been the reference for publications between 2004 to present. This code is currently used to deliver level‐1 derived neutron data (DND) to the Planetary Data System (PDS). These products are time series of corrected neutron counting rate data that can be used for scientific investigations. The level‐1 dataset includes averaged neutron data (AND), which consists of neutron maps built from the neutron time series data. The most recent processing and analysis of MONS data was performed by Maurice et al. [1] and covered data through July 2009. This led to work studying the depth-distribution of water on Mars [2, 5].

On the way to any science interpretation of inter-annual variability in the MONS dataset, this paper intends to document and provide the necessary elements for understanding the new data processing method and resulting dataset. We then build on this stage by presenting averaged counting rate maps and a preliminary comparison of the inter-annual variability in the Mars polar regions over 8 Mars Years by means of the neutron counting rates. These can be compared to previous results presented in [4, 6, 1].

2 MONS Instrument

The MONS instrument consists of an 11$\times$11$\times$10 cm3 BC454 plastic scintillator separated into four optically isolated segments, or “prisms.” This plastic is loaded with 5% natural boron by weight, which provides sensitivity to thermal and epithermal neutrons through neutron capture on 10B. The predominant interaction that occurs is 10B(n,$\alpha$)7Li*∗* with a Q value of 2.8 MeV. Due to inefficiency in light production from the heavy isotopes produced in this reaction, this energy is quenched and detected at 98 keV electron equivalent (keVee). A schematic of the MONS instrument is shown in Fig. 1.

There are two types of primary data products produced by the MONS instrument that define what type of neutron was detected [1]. Category 1 events (thermal, epithermal) are defined by a prompt interaction with an energy between 40 keVee and 630 keVee that is not followed by a delayed interaction within 25.6 $\mu$s. Category 2 events (fast) are defined as a similar prompt pulse with an expanded energy range of 40 keVee to 2.55 MeVee followed by a delayed pulse with an energy between 40 keVee and 630 keVee within a 25.6 $\mu$s window. For both Category 1 and Category 2 events, only one or two prisms can detect the event, otherwise the event is thrown away. Prompt events with an energy greater than 2.55 MeVee are categorized as GCR events.

Fast neutrons (Category 2 events) are defined as neutrons with an energy $>$0.7 MeV [1] and can be detected by all four prisms. Category 1 events can be split into thermal and epithermal neutrons based on prism. Prism 1 faces the nadir direction and is covered with a 0.69 mm thick cadmium sheet, which absorbs neutrons below $\sim$0.4 eV. Therefore, Category 1 events from this prism are epithermal neutrons (0.4 eV - 0.7 MeV). As noted in [1, 6], due to the geometry of the prisms there are small gaps in the cadmium coverage allowing Prism 1 some thermal neutron sensitivity. Thermal neutrons are detected by exploiting the Doppler filter technique [7], which uses the fact that the spacecraft velocity (3.4 km/s) is faster than the velocity of thermal neutrons (most probable value of 1.9 km/s or 0.019 eV at Mars’ atmosphere temperature). Prism 2 is forward facing along the direction of spacecraft motion, and therefore detects both thermal and epithermal neutrons. Prism 4 faces backwards along the direction of spacecraft motion, and therefore only neutrons that have a velocity higher than the spacecraft velocity can be detected. This corresponds to neutrons in the epithermal range with an energy greater than 0.06 eV. Thermal neutron ($<$0.06 eV) rates are determined by subtracting the Prism 4 counting rate from the Prism 2 counting rate. Similarly, an alternate definition of epithermal neutrons (0.06 eV - 0.7 MeV) can be obtained from the Prism 4 counting rates. Finally, Prism 3 is shielded from Mars and therefore should be a good proxy for the spacecraft background. The sides of the prisms are also covered in cadmium.

The mapping phase of the MONS instrument began February 22, 2002 and has been operating nearly continuously since then, leading to 17 Earth-years of data. The most recent processing and analysis of MONS data [1] covered data through July 2009, corresponding to nearly 4 Mars-years of data ($L_{s}$ = 330 in Mars Year (MY) 25 to $L_{s}$ = 313 in MY 29, using the Mars calendar defined by [8]). While [1] showed some inter-annual comparisons of counting rates in the polar regions, their work focused primarily on creating an averaged CO2 frost-free map of two-layer water-equivalent hydrogen (WEH) based on the MONS data that subsequently was used in the most definitive MONS mapping of WEH and its depth distribution to date [5]. Another paper including inter-annual comparisons of the CO2 frost cap thickness for two Mars years towards the beginning of the MONS mission can be found in [6]. Here we present new data processing of the Category 1 MONS data that includes all data through the end of 2017 ($L_{s}$ = 108.3 in MY 34). This doubles the amount of data processed by [1] and quadruples the number of MY in a detailed comparison of inter-annual variability of the seasonal CO2 frost deposits in the polar regions.

3 New Data Processing

The data processing includes many steps to take the MONS data from raw binary data to prism counting rates registered with latitude and longitude. Much of the data processing follows and draws upon the work described in [1], but performed independently. Raw data were acquired from the Planetary Data System (PDS) Geosciences Node (pds-geosciences.pds.wustl.edu), which releases data quarterly for the Mars Odyssey mission and GRS instrument suite. The raw data or experimental data records (EDR) are organized into folders by calendar year and subsequently by day. Raw data for the MONS instrument are contained within the neutron_spectra files. Relevant engineering data for MONS are contained within the eng subdirectory. Information on the format of each binary EDR file is contained within the main label directory.

The MONS data are pre-packaged to contain ephemeris data in addition to the instrument data. Each data point is registered with a UTC time stamp and an “SLCK” clock value that is unique for each data point. The neutron data includes 64-channel histograms for Category 1 events and 32-channel histograms for the prompt (early) and delayed (late) Category 2 events. Counter data, which store the total number of counts over threshold in 19.75 second accumulation windows, include GCR, deadtime, and the number of, and which, prisms fired. There is additional information on the first 84 Category 2 events within each accumulation window, including time between the prompt and delayed pulses, and pulse heights. The data also contain sub-satellite latitude and longitude at the middle of each integration window and position and velocity of the spacecraft in different reference frames.

The raw data conversion was done using Python 3.5 and the unpacked data stored in a MySQL Database. Following conversion, data reduction takes place to remove bad data from the dataset, described in Section 3.1. After all bad data are removed, data corrections that result in the final dataset are applied, described in Section 3.2.

3.1 Data Reduction

There are several categories of “bad” data that must be removed before further processing can take place. These data cuts apply to both Category 1 and Category 2 MONS data, although only Category 1 data are processed here. The first and largest data cut is from solar energetic particle (SEP) events, that produce a large background in the prism counting rates. Stability cuts are also applied to the counter data, which remove outliers and transients in these datasets. Cuts on spacecraft orbit parameters are applied to also remove outliers or transients and remove data acquired during clock resets that corrupt our ability to normalize to counting rates. Finally, some additional data cuts related to various anomalous readings are applied. The final dataset contains only good data that passes all four of the following described cuts. A summary of how much data is removed by each cut is provided at the end of this sub-section.

3.1.1 SEP Event Cuts

Removal of SEP events is done manually by looking at the counting rate recorded by a dedicated GCR counter and removing periods of rate excursions. An example of the base procedure is described below for a SEP event in September 2004, shown in Fig. 2. The mean and standard deviation of the GCR counter (total counts in each 19.75 s accumulation window) is determined for 8 days before and 8 days after the event. Figure 2 shows black bands representing $\pm$3 standard deviations ($\sigma$) from the mean. The excursion is flagged as when the GCR counter extends beyond $\pm$3$\sigma$ from these means. To safely remove the full extent of each event, the SEP event cut range starts 4 hours before the start of the excursion and ends 4 hours after the end of the excursion. The final cut range is demonstrated as the gray shaded region.

There are several SEP events where the event cut method was adapted or event cut ranges were manually updated from the base method. These included SEP events that were low in strength or short in duration, but most frequently were when a decrease in the GCR counter was observed surrounding the peak of the SEP event (likely due to changes in the interplanetary magnetic field). This was observed in $\sim$20% of SEP events. In these cases if the dip was before the main SEP excursion, the start of the event cut was determined by eye. If the dip was after the main SEP excursion, the mean and standard deviation in counter from before the event was used to judge when the counter returned to nominal. An example of this type of event (July 2004) is shown in Fig. 3. This event exhibited the decrease in rates both before and after the event. The mean was determined from data between 6/26/2004 – 7/3/2004.

Ranges defining the removal of SEP events from the dataset are summarized in Table 1.

3.1.2 Stability Cuts

Stability cuts were applied to the GCR counters and the total counts in each of the four prism Category 1 histograms. The stability cuts are applied based on the deviation of each data point from a boxcar rolling median value. A rolling window is specified as the number of data entries to sum over, and the result is centered within the time range of the window. Based on the time scales over which observed rates can change, we chose to apply a “daily” rolling median window. Resampling data from 2003 through 2007 to a frequency of one day, the typical number of data entries in one day was 4185 entries.

An example of this technique is described using data from 2004. Figure 4 shows a histogram of the deviations from the rolling median value for the GCR counter (left) and the Category 1 Prism 4 total histogram counts (right). The line indicates the stability threshold, which was chosen as ten times the median of the deviations (sometimes called the “MAD”). Data greater than the stability threshold are cut, indicated by the + markers in Fig. 5. The stability cut most obviously affects what seem to be spurious readings in the GCR counter. For the Prism total histogram counts the stability cuts remove most spikes observed in the rates.

3.1.3 Spacecraft Cuts

Cuts were applied based on the spacecraft orbit data to flag when a variable goes out of a range or experiences an excursion from nominal. First, spacecraft orientation data were required to be available in the raw data (flags for pointing and intersecting equal to 1). The parameters subject to spacecraft cuts that must be between certain values are the latitude ($-90^{\circ}$ to $+90^{\circ}$), longitude ($0^{\circ}$ to $360^{\circ}$), and altitude (380 km to 460 km). In addition, issues related to the Northerly equatorial crossing is flagged under this cut. As described in [1] and illustrated in Fig. 6, when the spacecraft is moving Northward and crosses the equator the internal clock is reset and thus the measurement time interval is lost. These events cannot be properly normalized into counting rates and are removed. The most common bounding cuts come from pointing and intersecting data not being available.

Finally, some instances of transient deviations in the mars position and velocity as recorded in the instrument frame were observed and removed.

3.1.4 Other Cuts

There were a few data reduction cuts that do not fall under the above categories. This includes something identified in [1] as erroneous latitude registration errors. This can be seen in Fig. 7, which shows the latitude (for any longitude) registered near the equator. The black points indicate when the spacecraft is moving Northward (“up”) and the red points when it is moving Southward (“down”). Data in this plot have all other cuts already applied, so gaps in time mostly represent SEP event cuts and the effect of removing the Northward equatorial crossing can also be observed. The two regions of time cut under this error code are both in 2002, and are also shown in Fig. 7 with the cut region shaded gray.

There was also very rarely an issue with the UTC timestamp showing up at the seemingly incorrect time, affecting 63 data points overall. For example, several data entries with UTC timestamps indicating 2007 showed up at SLCK values near November 25, 2006. There are two other data points with a similar issue, one in 2009 and one in 2010. It is likely the data affected by this anomaly are fine and somehow only the UTC timestamp is erroneous. However, to be safe these data were removed.

Finally, cuts based on the MONS sensor temperature and high voltage power supply (HVPS) were applied. Note that there appears to be a mismatch in the mapping of engineering data to files starting April 1, 2004. Before this date, the engineering data are in the appropriately labeled files, which for the HVPS are the hvps_mntr_[1,2].dat files. From April 1, 2004 through the dataset that we considered, the data corresponding to the HVPS1 and HVPS2 were empirically found to be in the files plus_5v_crnt_dig.dat and plus_5v_anlg.dat, respectively. The other files in the eng directory are similarly mis-mapped. We did not determine the mapping for the other engineering data.

The temperature and HVPS voltage are shown for the beginning of data collection through the end of 2017 in Fig. 8. These data are registered at different SLCK clock values than the rest of the MONS data therefore no event cuts are applied in these plots. The temperature fluctuates seasonally with some excursions to lower temperature. The origin of the double-banded structure in the temperature data is not known at this time. The two HVPS channels track together with some excursions to high voltages. All of the HVPS excursions occur during SEP events and are therefore excluded based on that cut. Some of the temperature excursions partially overlap an SEP event, but those that are not are removed.

3.1.5 Summary of data cuts

From February 22, 2002 through the end of 2017 there are just over 23 million data points processed. The total percentage of data removed from all cuts is 14.1%. The majority of this, 12.1% of the data or 85.6% of the total cut, comes from the removal of SEP events. Stability cuts remove only an additional 0.2% of the data. Spacecraft cuts affect 1.8% of the data, but 12% of the data under this cut overlap with SEP event ranges. Within the spacecraft cuts, 0.56% of the data is cut due to bad pointing/intersecting flags, 0.4% of the data exhibits transients in the position or velocity parameters, and 0.28% of the data are removed from equatorial crossings. Temperature cuts remove 0.65% of the data and the erroneous latitude registration early in the dataset affects 0.26% of the data.

3.2 Data Corrections

Several data corrections are applied after bad data are removed. These corrections are necessary to extract the correct prism counting rates and appropriate latitude and longitude registration.

3.2.1 ADC Non-Linearity

The analog to digital conversion of the prism spectra introduces differential nonlinearities into the recorded histogram data. To see this best and to determine the correction, Prism 3 histograms were averaged poleward of 85*∘N during the Northern summer. Kilometer-thick perennial water-ice deposits cover most of the region poleward of 80∘*N in the Northern summer [9]. A cut selecting the period between solar longitude $L_{s}=110^{\circ}-140^{\circ}$ was used to be safely in Northern summer based on [10] observations of Northern seasonal cap growth and retreat. Prism 3 is oriented away from the planet and therefore is expected to show only a smooth, continuous background. Deviations from this smooth function allow the nonlinearity correction to be determined.

The dataset analyzed for determining the nonlinearity correction (2002–2007) covers three Northern summers within this $L_{s}$ range that span $\sim$63 days each: 12/18/2002 – 2/19/2003, 11/4/2004 – 1/6/2005, and 9/22/2006 – 11/23/2006. The raw Prism 3 spectra for these three summers are compared in Fig. 9.

There are some differences year-to-year, likely due to slight differences in sensor temperature and GCR flux. However, overall a repeating 16-channel pattern can be observed, a result of the ADC design used for the Category 1 histograms. The data were smoothed by applying a centered 7-channel boxcar filter. This means the first three channels and last three channels are not included in the smoothed data, which are shown in Fig. 10.

The correction factor was determined by calculating 1 + (Smoothed[i] - Raw[i])/Raw[i] for these 58 ADC channels. The extracted correction factors are shown in the left panel of Fig. 11. To determine the final correction factor, data from the appropriate channels above channel 16 are averaged (e.g. channels 16, 32, and 48). Channels below 16 were excluded to avoid any artifacts from the boxcar smoothing near the peak. The final correction factor for each set of 16 channels comes from averaging the results of the three summers and is shown in the right of Fig. 11 and given in Table 2.

Figure 12 shows an example of the ADC nonlinearity correction applied to Prism 1 counts. The data come from $\sim$12 days during the peak of Northern winter ($L_{s}=266^{\circ}-274^{\circ}$) and are averaged over this time period for latitudes poleward of 85*∘*N. The total number of counts is increased by 0.7%, consistent with the average correction factor. We expect other systematic errors to be much larger than this and therefore are not concerned.

3.2.2 Gain Correction

The gain of each prism drifts during the duration of the mission, due to degradation over time, high voltage, or temperature variations. To determine the gain correction the position of the 10B neutron capture peak must be identified. Unless there are abrupt changes in the prism high-voltage, the peak position is very stable and changes smoothly with time. To acquire good statistics on the determination of the peak location, and the shape of the background that is used later, entries within a 10*∘* $L_{s}$ bin were subdivided into 40 subbins. The peak location and background parameters were determined on the histograms summed within these subbins, which typically contained approximately 1500 entries.

The fitting procedure used in this analysis includes 5 steps and is illustrated using Prism 4. Steps 1 and 2 (Fig. 13) follow similarly to what was done in [1]; a linear fit in log-log space is performed and then a Gaussian is fit to the rotated data. In these steps only part of the spectrum is fit, from a minimum ADC channel as determined by the maximum value plus an offset (this safely places the minimum fit after the low-channel roll over) to a maximum ADC channel of channel 28. The offset was prism- and time-dependent, as the gain of the neutron capture peak varies, with typical values between 1–3 channels. The maximum ADC channel of 28 was chosen to be high enough above the peak to ensure a good fit but below a background feature observed at high ADC channels. The initial guess for the mean value of the Gaussian fit in Step 2 is determined by finding all the local maxima in the rotated array and requiring that the mean value is not immediately within a certain number of channels to the minimum fit location and of the remaining candidates is selected as the one with the highest counts value. This results in a single and most often correct guess for the mean value. In Step 3 the Gaussian fit mean and sigma are used to create an excluded region around the peak of $\pm$3$\sigma$ from the mean, and the linear fit is repeated over the same constrained minimum and maximum ADC fit values in an attempt to improve the linear fit parameters. Step 4 repeats the rotation based on the tuned linear fit and refits a Gaussian to extract tuned Gaussian parameters. The tuned Gaussian mean is the peak location and is used to apply the gain correction that lines up all data to have a peak in channel 10.

The peak positions for the four prisms from the start of the mission through the end of 2017 are shown in Fig. 14. There is a period of approximately 3.5 months during MY 31 (8/4/2012 – 11/16/2012) that the prism gains are too low and peak values cannot be reliably extracted; this data has been removed.

3.2.3 Peak Integration

Once the histograms have been gain corrected, background parameters can be fit and the histograms integrated to determine the signal counts. The background fit parameters were constrained by using the same summed histograms as used to determine the gain correction. The tuned log-log linear fit parameters and the gain-corrected Gaussian fit parameters from Steps 3 and 4 above are used to guide initial guesses of a Gaussian plus background fit of the spectrum in log-log space. The results are compared for a linear background and a quadratic background in Fig. 15. Both background fits yield essentially the same mean peak position, however, the quadratic background fits the curve better and yields a better background subtraction, in particular when the peak location is high within the ADC range. When fitting each individual data point, the Gaussian mean, Gaussian sigma, background slope, and 2nd-order background polynomial coefficient are constrained by knowledge of the higher-statistics fits, while the Gaussian height and background offset are allowed to float. For each data point, this new fit results in parameterizations for the signal and background functions, which are used to integrate and determine the number of signal counts. No deadtime correction is made, as [1] showed this is quite small, and therefore all counts are simply divided by 19.75 s to determine the count rate per second.

3.2.4 Altitude Correction

The orbit of MONS is slightly elliptical, with an altitude that changes from $\sim$380–460 km. Figure 16 shows a histogram of the spacecraft altitude for all of the processed data from the start of the mission in 2002 through 2017. The counting rates are normalized to an altitude of 400 km using the following equation which corrects for the solid angle observed at a given altitude $h$ [11]:

[TABLE]

where $R=3389.5$ km is the mean radius of Mars. The scale factor is calculated as $\Omega$(400 km)/$\Omega$(h) and varies from $\sim$0.99–1.05. We do not make any corrections for local elevation of the surface.

3.2.5 Ground-Track Correction

The MONS detector is positioned such that Prism 1 faces nadir and Prisms 2 and 4 are along and opposite to the direction of spacecraft motion, respectively. Because of the motion of the spacecraft and the fact that thermal and epithermal neutrons have similar velocity to the spacecraft, the latitude and longitude that each prism “sees” (i.e. ground track) is not the same as the sub-spacecraft location that is registered with each data point. To determine the angle offset between the registered location and the observed location, the counting rates summed over longitude as a function of latitude when the spacecraft is traveling Northwards and Southwards can be compared. An observable shift is seen in the Northward versus Southward data, as illustrated in the left panel of Fig. 17 using data from the first quarter of 2004. The angle that resolves this discrepancy is the angle offset, and this depends on the prism. As the spacecraft orbit has an inclination of 93.2*∘*, both the latitude and the longitude have to be corrected.

To determine the angle offset that best resolves the discrepancy for each prism, a $\chi^{2}$-minimization was performed which is shown in Fig. 18. The resulting angles are 0.4289*∘* for Prism 1, 2.5297*∘* for Prism 2, 0.8757*∘* for Prism 3, and -1.5667*∘* for Prism 4. A positive number indicates the prism is seeing ahead of the spacecraft motion, while a negative number is behind the spacecraft. These numbers are within 10-20% to those published in [1].

This correction is applied to each data point where we first calculate the bearing (direction) of the ground track by using the latitude and longitude of the current and subsequent data points. Given the latitudes $\theta_{1,2}$ and longitudes $\varphi_{1,2}$ of the initial point (1) and final point (2), the bearing $\lambda$ is given by:

[TABLE]

This is combined with the knowledge of the angle offset determined from the $\chi^{2}$ optimization ($\delta_{i}$ for Prism $i$) to calculate the new latitude and longitude to register with each prism:

[TABLE]

Applying these corrections to the first quarter of 2004 shows a resolved offset in the Northward and Southward counting rates versus latitude, as shown in the right panel of Fig. 17.

Rotations of the spacecraft about the velocity direction were performed during several periods in 2009 – 2011. This does not impact Prisms 2 and 4 which are oriented along the direction of motion, however, data from Prism 1 has been removed during these periods since it does not face nadir. The rotations occur between MY 29, $L_{s}\approx 265$ through MY 30, $L_{s}\approx 34$ and between MY 30 $L_{s}\approx 316$ through MY 31, $L_{s}\approx 61$.

3.2.6 GCR Correction

The GCR flux, which is the source for the measured neutron signals, varies over time mostly with solar cycle. To remove this effect, the data must be normalized to a particular date with a known GCR flux. To determine the GCR correction factor, the belly band procedure [1] is adopted, which assumes that near the equator the ground-surface processes are in equilibrium and therefore the neutron counting rates should be stable over time. In addition to the GCR flux changing over time, seasonal changes in the density of Mars’ atmosphere can lead to seasonal changes in the neutron counting rates. These must be accounted for through simulations before the GCR correction can be determined. A radiation transport tool to simulate the neutron leakage flux from Mars and a tool to transport this flux to the MONS spacecraft and predict the prism counting rates were developed. The details of these tools are discussed in the Appendix.

We divided the MONS data into $2^{\circ}\times 2^{\circ}$ latitude and longitude bins within $\pm$20*∘* of the equator. The Mars Climate Database (MCD) v5.3 global circulation model (GCM) [12] was used to determine the atmospheric density on the latitude and longitude grid as a function of seasonal $L_{s}$. The neutron counting rates were normalized to an average atmospheric density of 16 g/cm2, using the simulation described in the Appendix. With seasonal effects from the atmosphere removed from the data, the GCR correction was then determined.

We chose to normalize the GCR proxy to be unity in June 2008. This corresponds to a solar modulation of $\phi=463$ MV according to the latest Usoskin model [13] and was chosen because it is during the time period with the lowest uncertainty in the determination of the solar modulation. This is different than [1] where the data were normalized to the period October – November 2002. Based on the difference in solar modulation and therefore integrated GCR flux, we expect the counting rates presented in this work to be up to a factor of 2 higher than the counting rates determined in [1].

The multiplicative GCR correction is shown in Fig. 19. The GCR correction is prism- and time-dependent, as it depends on each prisms efficiency and response to the time-varying GCR flux, which depends on the prism gain. Since the chosen GCR normalization date is close to solar minimum, when the GCR flux is largest, the GCR correction factor is generally greater than 1. Between solar minimum and solar maximum, the GCR flux can change by a factor of 2.5.

4 Results

As an example of the temporal coverage available in the new dataset, the prism counting rates in the polar regions are shown in Fig. 20 for the entire dataset in 10*∘* $L_{s}$ bins. Data from the North pole (latitude $>$80*∘N) is overlayed with the data from the South pole (latitude $<-$80∘*N). These plots can be compared with the results from [1] (Figs. 14 and 23), which show similar plots extending part way through MY 29. As CO2 frost is deposited seasonally the neutron counting rates increase, due to CO2 having a low cross section for absorbing the GCR-induced neutrons, until the seasonal CO2-ice cap reaches peak mass. The counting rates then decreases as the CO2 frost is sublimed away.

Qualitatively the same trends are observed in this dataset when compared with [1]. The baseline count rate in the Northern summer is lower than the Southern baseline count rate, which relates to differences in the regolith and perennial ice caps at each pole. The counts at the peak of seasonal CO2-frost deposition in the South are higher than the peak counts in the North. There is also a slight reduction in the Northern peak counts in MY 28–29 relative to previous years observed in Prism 2 that is not seen in Prism 1, similar to observations in [1]. Quantitatively as mentioned previously, differences in the GCR normalization result in the present counting rates being larger than those in [1]; on average the Prism 1, 2, and 4 counting rates are a factor of 1.4, 2.6, and 1.7 higher than [1]. The larger difference in the Prism 2 counting rates is attributed to using a quadratic background fit instead of a linear fit in the determination of the counts, as Prism 2 was observed to have the largest 2nd-order polynomial coefficient. The ratios of peak counting rates in the South to North are generally consistent with [1].

The new dataset can be averaged over all years to produce averaged neutron counting rate maps for each of the prisms. These averaged maps are shown in Fig. 21 assuming 1*∘* binning in latitude and longitude. The data plotted are only CO2 free data, assuming a cutoff of 0.2 g/cm2, similar to [1]. Perennial and seasonal CO2 are not distinguished in this cut, therefore data from the permanent CO2 cap in the South are not included in these plots. The CO2 frost thickness was predicted for each latitude and longitude bin as a function of $L_{s}$ using the MCD GCM model [12]. These maps are very similar to those found in Fig. 12 in [1]. The counting rates in these frost-free maps are inversely proportional to water content in the near-surface.

Since Prism 1 and Prism 4 provide different measures of epithermal neutrons (and with Prism 1 having a small contamination from thermal neutrons), we plot the correlation of the frost-free counting rates in Prism 4 to the counting rates in Prism 1 in Fig. 22. There is a strong correlation that is well-fit by a straight line, with a slope of 0.54 in this dataset and a slope of 0.48 from the Maurice et al. dataset [1], which is shown for comparison. The effect of differences in the normalization of the data in this analysis is evident. Since Prism 1 has a larger dynamic range and smaller errors (see discussion of errors and Fig. 25 below), it is a better choice as the epithermal detector when frost free data are considered. However, as discussed in [6], when studying the polar regions the thermal neutron counting rate is extremely sensitive to changes in the atmospheric abundance of N2 and Ar, which can vary seasonally. Therefore, differences in the epithermal neutron behavior between Prism 1 (which has some thermal neutron sensitivity) and Prism 4 may provide some information on how the atmosphere changes seasonally, and this should be studied in future work.

Stereo-graphic polar projections of the averaged thermal neutron counting rate (Prism 2 - Prism 4) are shown in Figs. 23 and  24 for the South and North poles, respectively.

Each plot shows the thermal neutron counting rate during local summer (left) and local winter (right). The maps extend to a latitude of 45*∘* in the respective hemisphere and are smoothed with a Gaussian filter with 5*∘* FWHM to remove random variations between pixels. The South pole plots in Fig. 23 use cuts of $L_{s}=330^{\circ}-360^{\circ}$ to define Southern summer and $L_{s}=140^{\circ}-180^{\circ}$ to define Southern winter, similar to the plots in Fig. 2 of [4]. The North pole plots in Fig. 24 use cuts of $L_{s}=110^{\circ}-150^{\circ}$ to define Northern summer and $L_{s}=330^{\circ}-360^{\circ},0^{\circ}-20^{\circ},$ to define Northern winter. The summer plots for both poles can also be compared to Fig. 25 in [1], which shows frost free polar maps of thermal neutron counting rates down to $\pm 60^{\circ}$. In the winter, the seasonal caps are clearly identified; in the Southern hemisphere, the peak in neutron counting rate is slightly offset from center towards the Northwest in the plot. In the summer, the South pole exhibits the perennial CO2 cap that is offset to the West-Northwest in the plot. Outside of this cap, the counting rates poleward of $-60^{\circ}$ are generally much lower than the lower latitude terrain. This is similar in the summer at the North pole, and small enhancements in the counting rate above $60^{\circ}$N are similar to observations in [1].

This dataset is intended to be used for studying seasonal effects and comparison of inter-annual variability, therefore the typical data product will be a count rate binned not only in latitude and longitude, but also in year and seasonal $L_{s}$. Therefore, the typical uncertainty in each bin, defined as the standard deviation of the data divided by the square-root of the number of entries within each bin, will be larger than [1] which focused on removing frost effects to produce global time-averaged count rate maps. To limit the uncertainties in a given bin to less than 10%, a limit of 10*∘* binning in $L_{s}$ and 4*∘* binning in latitude and longitude is required. The uncertainty over all years for this binning is shown in Fig. 25. The average uncertainty (solid lines) vary slightly with latitude, due to there being more spread in the data due to frost effects, and turn up at the highest point because no data exist poleward of $\pm$87*∘*. The shaded region represents the extent of the inner 80% of the data.

Similar plots can be made for frost free data only (like Figs. 19 and 26 in [1]), and the errors in the polar regions will increase slightly due to a reduction in the number of data entries available, but the spread in the uncertainty will decrease. Assuming the data is evenly spread over the 8 Martian Years, the errors shown in Fig. 25 will increase by a factor of 2.8 ($\sqrt{8}$) when considering a single year, leading to errors of $\sim$4%–7% for Prism 1, $\sim$3.5%–5% for Prism 2, and $\sim$5.5%–10% for Prism 4.

We generated a map with 10*∘* binning in $L_{s}$ and 4*∘* binning in latitude and longitude and plotted counting rate trends as a function of $L_{s}$ to perform a preliminary comparison of inter-annual variability. The counting rate averaged over all Mars years versus $L_{s}$ for latitude bins in the polar regions are shown in Fig. 26, 27, and 28, for Prism 1 (epithermal neutrons), Prism 4 (alternate epithermal neutrons), and Prism 2 - Prism 4 (thermal neutrons), respectively. The counts are integrated over all longitude and the latitude bin noted in the legend is the center point within the 4*∘* bin. The counting rates have a rough summer-time background subtraction applied by using the average of counting rates during summer so that the counts above the summer-time baseline counting rate can be compared across latitudes. Both epithermal analogs have similar trends, with the Prism 4 counting rates about a factor of 2 lower than the Prism 1 counting rates. The thermal counting rates have a much different trend and much higher counting rate. At this stage, atmospheric effects have not been removed from the data.

The counting rates as a function of $L_{s}$ are separated by MY for the 86*∘* latitude band in Figs. 29, 30, and 31 for Prism 1 (epithermal neutrons), Prism 4 (alternate epithermal neutrons), and Prism 2 - Prism 4 (thermal neutrons), respectively.

These plots also show the counting rate averaged over all MY for reference as a solid line. In general, the counting rate trends from year to year are very reproducible and agree within the uncertainties of the data. The notable exception to the reproducibility, which was previously mentioned, is the thermal neutron counting rate from $L_{s}\sim 315^{\circ}$ in MY 28 through $L_{s}\sim 45^{\circ}$ in MY 29. At the peak of Northern winter, the MY 29 counting rates are about 14% lower than the average of the other Mars years. This drop in counting rate, which is suggestive of less CO2 deposition, occurs after the planet-wide global dust storm that emerged around $L_{s}\sim 265^{\circ}$ in MY 28 [14, 15]. Observations in the Mars Climate Sounder thermal infrared data [10] show the Northern seasonal cap is $\sim$10% smaller in spatial extent following the MY 28 global dust storm. In future work, it will be interesting to process and analyze the MONS data through 2018, as another global dust storm with an onset in $L_{s}=185^{\circ}$ occurred in MY 34.

5 Summary & Future Work

In summary, we have performed a full re-analysis of Mars Odyssey Neutron Spectrometer data, extending the coverage of MONS data through the end of 2017 to cover 8 Martian Years. This paper summarizes the data analysis procedure, including data reduction and data corrections, to document and provide the necessary understanding to process raw MONS data from the PDS and utilize this new dataset. Example results based on this new dataset include frost-free global counting rate maps and maps of the polar regions and were presented in Section 4. These maps were qualitatively compared to previous analyses of the MONS data performed in [4, 1] and found to show similar trends. Due to choices made in processing of this dataset, the overall normalization of the data is different than the previous analyses. We showed the averaged counting rates for different latitude bands in the polar regions, which show the typical latitude dependence of the counting rates as a function of $L_{s}$. Preliminary results on inter-annual variability of the seasons CO2 caps were also presented in Section 4, which show reproducibility in the Southern seasonal cap based on both thermal and epithermal neutron counting rates, and reproducibility in the Northern seasonal cap based on the epithermal neutron counting rates and in thermal counting rates in all years except MY 28 - MY 29, following a global dust event. The variation in the thermal neutron counting rates in MY 28 - MY 29 may be due to a combination of cap properties and atmospheric properties, and will be explored in future work that will include an inter-annual comparison across all the latitude bins.

Work utilizing this new dataset is ongoing by this team. We are currently performing the necessary simulation and modeling efforts to normalize the data and convert counting rates to CO2 frost thickness in both the North and South polar regions and better understand atmospheric effects. With these efforts, we will be better able to interpret the overall properties of the seasonal caps and how the MY 28 global dust storm impacted the overall mass of CO2 deposited and any changes in the extent in the Northern seasonal cap following this event. This work will be the subject of future papers.

6 Data Availability

The unbinned time-series dataset related to this article can be found at https://doi.org/10.5281/zenodo.3267970, hosted at Zenodo [Mesickdata]. The neutron counting rate maps binned by mars year, mars solar longitude, latitude, and longitude, will be published in the Planetary Data System Geosciences node in 2020 at the conclusion of funded NASA work. Questions related to the use of the data can be directed to the corresponding author.

Acknowledgments

Research presented in this paper was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20160672PRD3 and the NASA Mars Data Analysis Program under project number 80HQTR170004. This research used resources provided under the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 8923328CNA000001. The authors would like to acknowledge the reviewers for their helpful comments in improving the manuscript.

Appendix A: Simulation Tools

Neutron Flux Signal

The radiation transport simulation package Geant4 [16] was utilized to develop a tool to simulate the expected neutron leakage from Mars at the top of the atmosphere. Benchmarking of this simulation package [17] has shown that with the appropriate choice of physics model, Geant4 and MCNP6.2 [18] produce similar results but both slightly over-predict neutron density profiles measured in the Lunar Neutron Probe Experiment [19]. The shapes of the curves generally match the data well, and therefore with proper normalization these types of simulation tools are helpful in the interpretation of planetary neutron data.

For the simulations relevant to determining the change in neutron flux with atmospheric density, a single-layer model with varying amounts of water-equivalent hydrogen (WEH) and varying atmosphere thicknesses were performed. The soil composition was based on the S21 average composition from [20], with elemental compositions given in Table 3. A fixed Cl abundance of 0.517% was added to this composition, close to the average as measured by the Mars Odyssey gamma-ray instrument. The size of the atmosphere normalization correction to the counting rate depends on the average WEH of the soil, and published maps of average WEH from the Mars Odyssey gamma-ray detector that are available in the PDS were used to estimate this within the belly band region. After adding in Cl abundance and the appropriate amount of H2O, the S21 elemental abundances were scaled uniformly so that the total elemental abundance summed to unity. The density of the soil was assumed to be 1.8 g/cm3.

The composition of the atmosphere in the simulation was based primarily on values used in [4], which come from the Viking data [21] with minor modifications. The concentrations of CO2, H2O, N2, and Ar from [4] are 96.93%, 0.054%, 2.7%, and 1.6%, respectively. In addition, an O2 concentration of 0.13% was assumed, taken from [21]. The atmosphere was simulated out to 40 km above the surface, with 20 layers of exponentially decreasing density (exp($-h/H$), where $H$ is the scale height) and a scale height of 10.8 km, similar to the layering used in simulations by [22].

As described in [23], gravitational binding of neutrons can effect the measured flux spectra. On Mars, the gravitational binding energy is 0.132 eV. Neutrons below this energy can return to the surface and re-interact. Gravitational binding of neutrons was implemented in our simulation by a reflecting boundary at the top of the atmosphere that reflected neutrons below the binding energy back to the surface. Since the surface return time $\Delta$t (derived in [23]) can be on the order of the neutron decay lifetime, a weighting factor exp($-\Delta t/\tau$) was applied based on the probability of neutron decay assuming the most recent value of the neutron free lifetime, $\tau$ = 880.2 s [24]. The epithermal and fast neutron flux are not affected by gravitational binding, however, the thermal neutron flux is almost 50% higher at the top of the atmosphere when gravitational binding is included.

Simulations were run for WEH values of 1%, 3%, 4%, 5%, 6%, and 8% (spanning the range of measured values in the belly band region) with total atmospheric thicknesses ($\rho_{a}$) of 4, 8, 12, 16, 20, and 24 g/cm2 (spanning the range of values predicted by the GCM model [12] in the belly band region) at each WEH point. In these simulations, the GCR flux was modeled following the Castognali & Lal model [25] with a solar modulation of $\phi=900$ MV. More details about this model and other GCR models can be found in [17]. At this stage, the neutron energy and angle at the top of the atmosphere are recorded. For the atmospheric correction since only ratios are being compared, the absolute GCR flux is not important.

Detector Response

The MONS instrument response was modeled in a separate simulation utilizing Geant4; the geometry is shown in Fig. 32. The geometry of the four boron-loaded plastic prisms was modeled, including cadmium covering where appropriate. The gaps between the Prism 1 face and the cadmium sheet where thermal neutrons can leak in was also modeled. The external spacecraft was not included in the model based on results from [4], which describes that Prism 1 is well shielded from spacecraft background, and that Prism 2 and Prism 4 have the same response so that subtraction of Prism2 - Prism 4 to determine the thermal neutron counting rate will effectively remove the spacecraft background.

The efficiency for neutron capture was simulated with an $11\times 11$ cm2 planar source placed directly on top of a given prism face, which fully covers the active area including the gaps. Neutrons were uniformly generated from random points within the rectangular dimensions. Neutron energies of 1$\times$10*-10* MeV and increasing in increments of one decade to 1 MeV were simulated, and at each energy the neutrons originated at angles from 0 to 85 degrees in increments of 5*∘. A simulation with the source incident on Prism 1 provides the efficiency for epithermal neutron detection. A simulation with the source on Prism 3 provides the efficiency for thermal neutron detection, and serves as a proxy for Prism 2 and Prism 4 through an angular transformation. A simulation with the source on the side of the prisms was also performed. In addition to the efficiency being tabulated for events hitting the primary prism in each configuration, the efficiency of cross-talk events, e.g. events detected in Prism 2, 3, or 4 when the source was on Prism 1, was also tabulated. Figure 33 shows the effective area at 0∘* incident angle for Prism 1 (epithermal neutrons), Prism 3 (thermal neutrons), the side of Prism 1, and examples of cross talk events detected in Prism 2 when the source is on Prism 1 or Prism 3. The effective area is the simulated efficiency multiplied by the source area, and the total effective area is the summation of the primary + cross talk effective areas.

The simulated MONS counting rate was calculated by taking the simulated neutron current at the top of the Mars atmosphere (40 km) and combining it with the appropriate effective area tables from the Geant4 detector response simulations. However, several steps in between occur to account for the ballistic trajectory of the particles, the gravitational binding of neutrons, and the spacecraft velocity. The procedure to calculate the count rate for each prism includes equations derived in [23] and follows the steps below.

Transform the energy and angle to account for the ballistic trajectory due to the gravity of Mars, 2. 2.

Transform the energy and angle into the spacecraft motion frame, 3. 3.

Interpolate the effective area of the incident prism from the final energy and angle of the neutron, and 4. 4.

Account for the neutron lifetime.

The energy and angle of the simulated neutrons were transported to the spacecraft orbit assuming ballistic trajectories. The final energy of the neutron is [23] $K_{r}$ = $K-V\frac{(R-{R_{M}})}{R}$, where $K$ is the energy of the neutron leaving the surface of Mars, $V=GMm/R_{M}$ is the gravitational binding energy of Mars (0.132 eV), $G$ is the gravitational constant, $M$ is the mass of Mars, $m$ is the mass of the neutron, $R_{M}$ is the radius of Mars, and $R$ is the distance from the center of Mars to the spacecraft orbital altitude. The neutron incident angle was adjusted for the ballistic trajectory caused by the binding energy of Mars [23]: $\mu_{r}^{2}=1-\left(\frac{R_{M}}{R}\right)^{2}\frac{K}{K_{R}}(1-\mu^{2})$, where $\mu$ is the cosine of $\theta$, the angle the neutron leaves the surface, and $\mu_{r}$ is the cosine of $\theta_{R}$, the angle of the neutron at orbit.

The new energy was transformed into the reference frame of the spacecraft depending on which prism the neutron was determined to hit. Using Galilean velocity transformations the constant speed of the spacecraft (3380 m/s) is added (subtracted) to the component of the neutron velocity in the direction of the spacecraft motion if the neutron hit Prism 2 (4). For Prism 1 and the prism side, a random sampling of azimuthal $\phi$ values were assigned to each neutron angle and energy, and the incident velocity was calculated with the $\phi$ value included. This was to account for all possible points the neutrons may hit the detector from, and the results were averaged. This correction affects what energy each Prism “sees” based on the motion of the spacecraft, which impacts the effective area. Prism 2 is heading “toward” the neutrons and will crash into them, adding to their incident energy and Prism 4 is traveling “away” from the neutrons taking away some of their energy upon impact. In addition to this effect, the neutron flux itself must also be corrected for the shift in detected neutron energy. The original neutron flux and velocity at the top of the atmosphere were used to calculate the neutron density. This neutron density was then multiplied by the velocity of the neutron at the spacecraft to get the neutron flux in the reference frame of the moving spacecraft. Finally, given the definition of angles in the detector response simulations, an angle transformation of 90*∘* is applied to events hitting Prism 2, Prism 4, and the prism side.

The effective area tables were interpolated in energy and angle using the Python 2D SciPy interpolation package. For events hitting Prism 1, the Prism 1 effective area tables were used. For events hitting Prism 2 or Prism 4, the Prism 3 effective area tables were used. For events hitting the side of the prism cube, the Prism 1 side effective area tables were used. The primary signals (events hitting only the primary prism, and ignoring contributions from the side) account for 65% of Prism 1 and Prism 4 events and 85% of Prism 2 events.

Finally, a weighting factor is applied to account for the probability of neutron decay during the transport of neutrons from 40 km to the spacecraft orbit, using a neutron lifetime of $\tau$ = 880.2 s [24] and the neutron transit time [23] (note there is a sign error in this equation as found in [23], which has been corrected here):

[TABLE]

where $A=\left[4(K/V)(1-K/V)\mu^{2}\right]^{1/2}$, $B=2K/V-1$ for $K/V<1$, $m$ is the mass of the neutron, and the other constants have already been defined.

GCR Correction Details

The resulting ratio of count rates for Prism 1, Prism 2, Prism 3, and Prism 4 as a function of atmospheric density ($\rho_{a}$) for the different WEH values is shown in Fig. 34. For low WEH content, the correction factor is as much as 19% for the lowest atmospheric density, however, the typical correction is much smaller. These curves were fit to second-order polynomials, leading to the fit parameters given in Table 4.

To correct out the changes in neutron counting rates due to seasonal variations in the atmosphere, the data were binned for each year in $L_{s}$. At each $L_{s}$ point, the MCD GCM [12] was used to determine $\rho_{a}$ for each 2*∘$\times$2∘* latitude and longitude bin within the belly band region. The correction factor based on this $\rho_{a}$ was then calculated for each of the six simulated WEH values using the fitted parameters. The published map of WEH derived from the Mars Odyssey high-purity germanium gamma-ray spectrometer [Boynton2007, GRSmap] was then used to determine the WEH of the soil within each bin, and the final correction factor determined by an interpolation of the correction factors covering the range of WEH values. Once the neutron counting rates were normalized to 16 g/cm2, the GCR correction factor was determined.