Three Years of ACB Phase Function Observations from the Mars Science Laboratory: Interannual and Diurnal Variability and Constraints on Ice Crystal Habit

MSL's Navcams are two sets of cameras mounted on the rover's remote sensing mast (RSM). They operate in the 600–850 nm band with a field of view (FOV) of 45° × 45° and a signal-to-noise ratio of 200:1 (Maki et al. 2012), making them well suited to observing optically thin Martian clouds. Additionally, the RSM is able to rotate 360° in azimuth and 178° in elevation (Maki et al. 2012), allowing the Navcams to observe an extremely wide range of angles.

The Navcam PFSS was instituted in MY 34 to be used in the ACB season. It consists of nine three-frame movies each at different pointings forming a dome around the rover (Figures 1 and 2 in Cooper et al. 2019). This allows observation of a greater range of scattering angles in order to characterize the single scattering phase function of ACB season water-ice clouds at Gale Crater and constrain ice crystal geometries. A more detailed description of the observation is given in Cooper et al. (2019).

The PFSS has a preferred cadence of one observation per week during the ACB season (L_s ∼ 40° to L_s ∼ 150°), alternating between morning and afternoon observations. It has been taken at this approximate cadence through the ACB seasons of MYs 34, 35, and 36.

In order to not oversaturate the Navcams by imaging the Sun, there is a 2.5 hr solar keep-out around local noon, and the morning and afternoon observations include a gap in their coverage (as seen in Figure 1 in Cooper et al. 2019). However, some frames, which were taken too close to the Sun, still proved unusable, as it was difficult or impossible to determine if clouds were present. As a result, we imposed a lower limit of 45° on the scattering angles we could confidently observe, based on the FOV of the Navcams.

The scattering phase function, P(θ), describes the angular distribution of scattered radiation—in the case of this work, sunlight scattered by water-ice clouds—as a function of scattering angle (θ). The phase-function expression used in this work is derived by Cooper et al. (2019) and is based on the radiative transfer equation for plane parallel atmospheres (Equation (1.4.22) in Liou 2002). Since ACB clouds observed from Gale Crater are optically thin, we can use a single scattering approximation. Additionally, we need to account for light extincted by dust between the clouds and the rover, which is done by multiplying by a factor:

$$\begin{eqnarray}&&\exp \left(-\displaystyle \frac{{\tau }_{\mathrm{COL}}}{\mu }\right),\end{eqnarray} \tag{ 1 }$$

where τ_COL is an optical depth taken from observations by the MSL Mast Cameras (Lemmon et al. 2024).

We also need to assume a cloud optical depth, (Δτ_WIC), which Cooper et al. (2019) did by taking Mars Climate Sounder (MCS) water-ice column optical depths for the area of Gale Crater. These values were averaged over 10° of L_s and the last five MYs and normalized to 880 nm by multiplying a factor of 3.6 (Guzewich et al. 2017), which assumes an ice particle effective radius of 1.41 μm (Kleinböhl et al. 2011).

While we initially used the same method as Cooper et al. (2019), using MCS optical depths for Δτ_WIC, we also computed the phase function replacing these optical depths with those derived from MSL cloud observations (Hayes et al. 2023). This has the advantage that the MSL values are calculated for specific cloud observations, whereas the MCS tau values are average values. In addition, most MSL cloud opacity observations are taken very close temporally to the PFSS observations. More discussion of the strengths and weaknesses of each tau data set is given in Section 4.1. We present the results of both methods in Section 3 for completeness.

The expression for the phase function is given by (Equation (5) in Cooper et al. 2019):

$$\begin{eqnarray}&&P(\theta )=\displaystyle \frac{4\pi \mu {I}_{\lambda ,\mathrm{VAR}}{\rm{\Delta }}\lambda }{{\rm{\Delta }}{\tau }_{\mathrm{WIC}}{F}_{\lambda }\exp \left(-\tfrac{{\tau }_{\mathrm{COL}}}{\mu }\right)},\end{eqnarray} \tag{ 2 }$$

where μ is the cloud viewing angle, I_λ,VAR is the variation in spectral radiance, Δλ is the bandpass of the Navcams, and F_λ is the in-band solar flux.

The water-ice clouds seen during ACB season with MSL are often too optically thin to be seen in raw Navcam images, so each movie undergoes mean-frame subtraction (MFS), which removes the time-invariant component of each frame, leaving a perturbation movie of the time-varying component. Figure 1 shows an example of this. Similarly to Cooper et al. (2019), we assume that scattering from dust is uniform across the Navcam images, and thus, it is removed by the MFS process. However, some component of time-varying dust may still be present in the MFS frames.

Zoom In Zoom Out Reset image size

Following MFS, each perturbation movie is contrast stretched and manually examined to determine if clouds are present. If they are, a radiance map of frame 2 is created from the original images to determine the variation in spectral radiance between a region of assumed cloud and a region of assumed empty sky (Figure 2). We select points of high spectral radiance (presumed cloud) and low spectral radiance (presumed empty sky) and average over a circle of 5 pixel radius to determine the value of I_λ,VAR to be used in Equation (2). The value of I_λ,VAR should be taken as a lower limit due to uncertainty in whether the low-radiance region is truly empty sky or merely a thinner cloud. The selected location of the cloud within the frame is also used to calculate the scattering angle, θ.

Zoom In Zoom Out Reset image size

While it is difficult to account quantitatively for all errors associated with Equation (2), we can approximate the uncertainty from several sources. The MSL opacities used for this study have a mean uncertainty of around 20%, while the uncertainty on the I_λ,VAR term is likewise ∼20%. The contribution to uncertainty from τ_COL is much lower, averaging ∼6% for the ACB seasons of MYs 34, 35, and 36. Taken together, these give an uncertainty of around 30% for our calculated phase function, when we use the MSL-derived cloud opacities.

Cooper et al. (2019) take an upper limit on MCS optical depths of 10%. Combined with the uncertainties on I_λ,VAR and τ_COL, this gives an uncertainty of around 23% for the phase function. However, this does not take into account any uncertainties that may be associated with the fact that we are averaging the MCS opacities over 10° of L_s and five MYs.

In determining an ice crystal habit, and in comparing to previously derived water-ice cloud phase functions, the overall shape of the phase function is of more interest than its absolute magnitude. Typically, a phase function is normalized over scattering angles 0°–180° such that the area under the curve is equal to 4π (Hansen & Travis 1974). For instance, HG functions (discussed in more detail in Section 4.2) are normalized in such a way (Hapke 2012a). However, since more scattering from particles of comparable radius to Martian aerosols occurs in the forward scattering direction (scattering angles less than 40°), which the PFSS is not able to observe, it is not possible to determine what normalization values are appropriate for our derived phase function. As a result, our computed phase functions are normalized around the mean values of the phase functions to which we compare them.

Figure 3 shows the temporal distribution of the last three MYs of PFSS observations. The MY 34 data are described in detail in Cooper et al. (2019). In MY 35, there were a total of 26 PFSS observations, with 13 morning and 13 afternoon observations, spanning L_s = 55.0°–159.4°= (sols 2471 to 2691). In MY 36, there were a total of 30 PFSS observations, with 13 morning and 17 afternoon observations, spanning L_s = 43.7°–164.0° (sols 3115–3368). Of these, the two penultimate observations (sols 3353 and 3359) were excluded as they were taken during a regional dust storm, and it was not possible to determine if the features seen in the frames were water-ice clouds or dust features. Aside from the above exceptions, every observation contained at least one movie in which clouds could be positively identified, for a total of 167 movies in MY 35 and 191 movies in MY 36 showing cloud features.

Zoom In Zoom Out Reset image size

The radiance points chosen from the movies and their corresponding calculated phase functions spanned scattering angles 46.4°–154.1° in MY 35 and 45.0°–139.9° in MY 36. These angular ranges are much narrower than that of MY 34, due mainly to our change in methodology that rejects scattering angles less than 45°.

We produced a mean phase-function curve for both MYs using the MCS cloud opacities, using a moving average with a window size of 15° of scattering angle. This window size remains fixed regardless of the number of points within the 15° bin. Figure 4 shows these curves with their 95% confidence intervals. Both phase functions have a greater magnitude at scattering angles >60°, which decreases with increasing scattering angle until the curve becomes approximately flat after about 100°.

Zoom In Zoom Out Reset image size

When we compare our two derived curves to the Cooper et al. (2019) MY 34 curve (Figure 5), we see that all three phase functions have a similar shape, with small differences in normalized magnitude generally falling within their respective 95% confidence intervals.

Zoom In Zoom Out Reset image size

We also derived mean phase-function curves using MSL water-ice tau values in the place of MCS column tau values. We rederived individual phase-function values and applied the same moving average to produce curves for all three MYs, as shown in Figure 6. As in the case of MY 35 above, we removed any outliers that had values greater than 3 standard deviations from the mean.

Zoom In Zoom Out Reset image size

Again, all three phase functions agree closely with each other (Figure 7), and we see the increase in magnitude in the forward scattering direction and the flattening of the curve as scattering angle increases. The MY 36 phase function shows some increase in magnitude toward 140°; however, this is likely due to the relative paucity of data points at that range of scattering angles.

Zoom In Zoom Out Reset image size

We were interested in exploring if there was any diurnal variability in ice crystal habit. Previous studies have shown the potential for diurnal variability at Gale Crater in both opacity (Kloos et al. 2018) and altitude (Campbell et al. 2020). These differences could be due to differences in temperature and formation mechanism, with morning clouds being a result of overnight clouds formation at lower temperatures and afternoon clouds formation due to enhanced convection (Tamppari et al. 2003). As a result, we wished to see if this would impact the cloud ice crystal habit. We derived mean morning and afternoon phase functions by taking morning and afternoon observations separately and applying the same moving average. Figure 8 shows a morning versus afternoon comparison for each MY using the MSL-derived tau phase function.

Zoom In Zoom Out Reset image size

For all three MYs, the shapes of the mean curves are very similar to each other, which does not indicate a difference in phase function between the morning and afternoon clouds. It is not possible to definitively say whether the ice crystal habit is constant over the course of a sol, but such behavior would be consistent with our data.

Like Cooper et al. (2019), we initially used MCS retrievals of column water-ice optical depth for the tau value of water-ice clouds to derive the phase function. However, our use of these MCS values does present some potential problems in our calculation of the phase function. Most obviously, the MCS tau values we use do not correspond one-to-one to the PFSS observations, and indeed the values we use are averaged over 10° of L_s and five MYs. As a result, we are not using tau values for the individual clouds we see in the PFSS observations but instead taking an average tau value for the apporximate time of year. Additionally, at low altitudes, the atmosphere becomes opaque and can become obscured by topography, which restricts the instrument's ability to detect water ice (Kleinböhl et al. 2009).

The use of tau values derived from MSL observations has several advantages over the use of MCS tau values. The MSL cloud observations from which the tau values are derived (zenith and suprahorizon movies, ZMs and SHMs) are often taken with the PFSS observation, and in most cases we were able to use tau values for clouds that were observed on the same sol and nearly at the same time as the PFSS observations. Only in a few cases were tau values taken from observations several sols (and thus only a few degrees of L_s ) before or after the PFSS. Additionally, the ZM and SHM also use the NavCams, meaning that the observations are in the same wavelength and have the same viewing geometry as the PFSS, namely looking at the clouds from below. This means in most cases we can be sure that we are observing the same clouds and deriving a tau value that is representative of those specific clouds.

That being said, neither approach is without weakness. The MSL-derived tau values are calculated using an assumed phase function, the weaknesses of which are discussed in Hayes et al. (2023). However, given the advantages in seeing either the same clouds, or at least temporally nearby clouds, we have made the choice to continue using the phase-function values calculated with MSL-derived tau values.

Cooper et al. (2019) attempted to constrain a dominant ice crystal habit through comparing their derived phase function with modeled phase functions of seven randomly oriented water-ice crystals (Figure 7 in Cooper et al. 2019). They found that the phase function did not fit any one of the curves for the ice crystal geometries particularly well. As our derived phase functions are very similar in shape to that of Cooper et al. (2019; Figure 9), we may expect to see similar results. However, it is likely that this method of comparison is not optimal for Martian water-ice clouds. The water-ice crystal data set used by Cooper et al. (2019) assumed very large water-ice particle sizes on the order of tens of microns (Yang & Liou 1996; Yang et al. 2003), while Martian water-ice crystals tend to be much smaller (Clancy et al. 2003). In addition, the formation of Martian water-ice clouds, with much less water vapor aggregating on a much greater number density of dust particles that serve as cloud condensation nuclei, may create some unusual or irregular geometries for which there is no good terrestrial analog. While the halo seen by Perseverance does strongly imply a crystal size and shape (Lemmon et al. 2022), it is a single occurrence that has not been seen in over five MYs of observations at Gale Crater.

Zoom In Zoom Out Reset image size

We instead use a different method of determining likely ice crystal habit, by fitting a two-term HG function to our calculated phase-function values. HG functions can be used to approximate single scattering phase functions, with a two-term HG function representing bidirectional scattering (Hapke 2012a). HG functions have a history of use in planetary applications though mostly for modeling the phase functions of dust and regolith (Hapke 2012b). Despite this, we believe that it is a good fit for our application as it is, at its core, used to represent single scattering phase functions, which we have derived for the PFSS results. Additionally, it allows us to easily compare our results to various airborne dust particles, which have been modeled in similar ways.

We use a two-term HG function similar to that used by Moores et al. (2015), who examined dust cover on the Phoenix Lander telltale mirror:

$$\begin{eqnarray}\begin{array}{rcl}P{(\theta )}_{\mathrm{HG}} & = & \displaystyle \frac{(1-c)(1-{b}^{2})}{{\left(1+2b\cos (180-\theta )+{b}^{2}\right)}^{3/2}}\\ & & +\,\displaystyle \frac{c(1-{b}^{2})}{{\left(1-2b\cos (180-\theta )+{b}^{2}\right)}^{3/2}},\end{array}\end{eqnarray} \tag{ 3 }$$

where c is the backscatter fraction and b is the asymmetry parameter. Small values of b give low, wide forward- and backscattering lobes, while larger values give higher, narrower lobes. Higher values of c indicate smaller amounts of forward scattering (Hapke 2012b).

Different authors have defined c in different ways. Hapke (2012a) uses a value of c, which gives the two-term HG function as

$$\begin{eqnarray}\begin{array}{rcl}P{(g)}_{\mathrm{HG}} & = & \displaystyle \frac{(1+{c}_{{\rm{H}}})(1-{b}^{2})}{2{\left(1-2b\cos (g)+{b}^{2}\right)}^{3/2}}\\ & & +\,\displaystyle \frac{(1-{c}_{{\rm{H}}})(1-{b}^{2})}{2{\left(1+2b\cos (g)+{b}^{2}\right)}^{3/2}},\end{array}\end{eqnarray} \tag{ 4 }$$

where g is the phase angle. Thus, our value of c can be related to the Hapke (2012a) value (c_H) by

$$\begin{eqnarray}&&{c}_{{\rm{H}}}=2c-1.\end{eqnarray} \tag{ 5 }$$

Hapke (2012b) created a plot of b versus c values for 495 regolith and dust particles, which lie within a region of physically allowable values resembling a hockey-stick shape (Figure 2 of Hapke 2012b). The centroid of this region is given by a "hockey-stick relation" (Equation (8) of Hapke 2012b).

We fit the HG function (Equation (3)) to the full range of data for each MY. Figure 10 shows these b and c values, with error bars representing the 95% confidence bounds in b and c. While the asymmetry parameter, b, appears fairly well constrained especially compared to c, this is only the case within the bounds of our data, between scattering angles of 45° and 155°. These limitations on our data set mean that we are not seeing the areas where we would expect to see most scattered light, and therefore the fit HG functions for each MY are only our best guesses with the data available.

Zoom In Zoom Out Reset image size

Our b and c values fall near the hockey-stick curve and within the allowable bounds established by Hapke (2012b; Equation (6) in that work). In addition, our values fall near those reported by Souchon et al. (2011), which looked at volcanic samples, particularly irregular, rough, opaque glass and round, smooth, translucent glasses. While the cloud particles we observe are not necessarily exactly these materials, we may see similar irregular, rough shapes. This could also indicate we are seeing igneous dust particles with ice accreted onto it.

As a point of comparison, we fit the same HG function to the phase functions of six modeled water-ice crystals with particle sizes around 3 μm (Bi & Yang 2014). These are different phase functions than were used by Cooper et al. (2019), which had much larger particle sizes. The ones we use are on the order of likely Martian water-ice crystal sizes (Clancy et al. 2003). We plot the determined b and c values alongside those we found for the PFSS results.

We also compare our results with those for Martian dust particles. HG functions have been used to fit Martian dust like that seen at the Spirit rover site (Johnson et al. 2006), the Phoenix lander site (Moores et al. 2015), and, most relevantly to this study, Gale Crater (Chen-Chen et al. 2019). Johnson et al. (2006) found an HG function for the Martian simulant HWMK919 at various wavelengths. We use their results at 750 nm, which falls within the bandpass of the MSL Navcams. The Chen-Chen et al. (2019) results also use observations from the Navcams but use a three-term HG function. We fit our two-term expression to their three-term curve in order to extract b and c values. The Surface Stereo Imager on board Phoenix used in the Moores et al. (2015) study is centered at 450 nm, outside the bandpass of the Navcams. We also include the b and c values for a theoretical perfect, clear sphere as determined by McGuire (1995) for comparison.

We see (Figure 11) that our b and c values do fall relatively near several of the other phase functions, especially the Navcam-observed dust. We also note that our phase function does not fit particularly well to that of the theoretical sphere. This is of note given that many climate models assume spherical water-ice crystal geometry. The closest of the modeled water-ice crystals to our own phase function is droxtals, which is of interest as that geometry is used in the single scattering phase function of Wolff et al. (2019).

Zoom In Zoom Out Reset image size

Due to our narrow range of scattering angles, there is some inherent uncertainty associated with our fit HG functions as we do not see either the forward- or backscattering lobes. Additionally, the method of fitting an HG function has been used primarily in planetary applications to look at the phase function of regoliths and dusts, and as such, there are not many good points of comparison for water-ice particles.