Five Mars Years of Cloud Observations at Gale Crater: Opacities, Variability, and Ice Crystal Habits

All of the cloud observations to be discussed were taken using MSL's Navigation Cameras (Navcams), of which the rover has two pairs (one pair each for the primary and backup computer systems) mounted to the rover's remote sensing mast (RSM). Although they are classified as engineering cameras, as their primary function is to assist in rover navigation, the Navcams are radiometrically calibrated (Maki et al. 2012) and thus can be used to conduct scientific observations.

The Navcam sensors are flight-ready spare charge-coupled devices left over from the Mars Exploration Rover missions. They have a useful photosensitive area of 1024 by 1024 pixels and have a 45° field of view (FOV). The optical filter bandpass of each Navcam is centered around 650 nm, with sensitivity between 600 and 850 nm. Due to the expected low optical depths of the clouds, the Navcams were selected for the cloud observation campaign over the rover's other cameras due to their superior operational signal-to-noise ratio (S/N > 200:1; Maki et al. 2012), despite being sensitive to a wavelength range that is optimized for observing dust-covered surfaces rather than water-ice clouds. In addition to their excellent S/N, the Navcams provide the widest FOV of any of the rover's cameras (see Figure 1), allowing us to observe clouds over as much of the sky as is possible in a single image.

Zoom In Zoom Out Reset image size

The MSL cloud observation campaign has grown and evolved substantially over the course of the mission, as can be seen in Figure 2. Moores et al. (2015) and Kloos et al. (2016, 2018) all focused on the two longest-running observations: the Zenith and Suprahorizon movies (ZMs and SHMs, respectively). Since the publication of Kloos et al. (2018), a number of new activities have been implemented to study the properties of the clouds over Gale that the ZMs and SHMs are not optimized to observe, such as cloud altitudes and the scattering phase function. An additional two new observations—the Shunt Prevention ENV Navcam Drop-In (SPENDI) and the high data volume SuperSPENDI variant—have allowed us to supplement the nominal SHM cadence.

Zoom In Zoom Out Reset image size

In Figure 3, we present the temporal distribution of the ZMs and SHMs as a function of solar longitude and local true solar time (LTST). From the beginning of the mission until the beginning of MY 33, nearly all ZMs were taken near 16:00 LTST, with the majority of SHMs taken near 12:00 LTST. Both were driven by operational constraints; early-morning observations required additional power to warm the Navcams, while the near-vertical pointing of the ZMs and MSL's equatorial latitude prevent ZMs from being executed within 2.5 hr of local noon due to the risk of having the Sun present in the frame.

Zoom In Zoom Out Reset image size

Beginning in MY 33, dedicated once-weekly early-morning science time was allocated for environmental science observations, including a ZM and an SHM. This allowed for nearly full daytime diurnal coverage, with the exception of a gap near 10:00 LTST caused by the rover's uplink schedule. Beginning approximately halfway between the MY 35 and MY 36 ACB seasons, the timing of the AM and noon science blocks has become increasingly constrained, reopening gaps in the diurnal coverage. These gaps are a consequence of the fact that as the plutonium dioxide heat source that powers the rover's radioisotope thermoelectric generator (RTG) decays, simultaneous operation of the rover's survival heaters, actuator and camera warm-up heaters, and science instruments puts a greater strain on MSL's battery resources. Thus, science operations are preferentially limited to times when heating is less of a concern, with less flexibility to deviate from those preferred times than there was in the past.

Because environmental science time is now rarely allocated in the late afternoon, it has become challenging to ensure that our diurnal coverage does not become biased toward the AM and near-noon periods. Since MY 35, we have been aided by the implementation of two new observations—the Cloud Altitude Observation (CAO) and the SPENDI, discussed more in Sections 2.2.1 and 2.2.2, respectively—which contain ZM/SHM observations but are less subject to the power constraints faced by the standard ZM/SHM observations. The CAO is part of the standard ACB season observational cadence and can only be executed in the afternoon, while the SPENDI is specifically designed to drain excess charge from the rover's batteries. The additional flexibility afforded to us by these observations means that we have been able to maintain reasonably robust PM coverage despite the declining output from the RTG.

2.2.1. ZMs

2.2.2. SHMs

Due to the low optical depths of the clouds observed over Gale, they are very rarely visible to the naked eye in the raw movie frames. To enhance the contrast between the clouds and the empty sky for ease of analysis, we use a technique known as mean-frame subtraction (MFS). We first construct a mean frame as a pixel-by-pixel average of all eight frames. This mean frame is then subtracted from each individual frame, leaving behind only the component of the image that varies with time (i.e., the clouds). Although all analysis is performed on the MFS data, we take the additional step of renormalizing each frame by dropping the highest and lowest 2% of the values and stretching the remaining values to assist in the visual identification of cloud features and the selection of high- and low-radiance points for our opacity retrieval models. An example ZM following each of these three steps is shown in Figure 4.

Zoom In Zoom Out Reset image size

Once each movie has been processed, it is assigned a quality value between 0 (featureless) and 10 (clouds visible in the raw frames). A value of −1 is also used to indicate movies that are unusable due to either the presence of the Sun in the frame or an RSM failure that left the Navcams pointing in a direction where the sky could not be observed. Movies with quality values of −1 or 0 are excluded from further analysis. In Figure 5, we show the quality values assigned to each movie, averaged across 50 bins of solar longitude and LTST. The ACB season is clearly visible as the band of higher-quality movies between L_s ∼ 50° and 170°. A second, smaller peak centered around L_s ∼ 300° represents the peak of atmospheric dust loading during the perihelion dusty season.

Zoom In Zoom Out Reset image size

To determine the cloud opacities, we use two models that make different assumptions about the nature of the clouds being observed. The first, reproduced here as Equation (1), is known as the High Clouds (HC) model (see Equations (2)–(9) of Kloos et al. 2016 for a full derivation). The HC model assumes that the clouds are at high altitudes above the bulk of the dust-scattering centers and are composed of water-ice crystals:

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

Stepping through this equation term-by-term, μ is the cosine of the zenith viewing angle and I_λ,VAR is the difference in spectral radiance between a high- and low-radiance point in an MFS frame. These two points, corresponding respectively to a cloudy and cloudless portion of the image, are chosen through visual inspection. Note that the low-radiance point could be a region of the sky with lower-opacity clouds rather than open, cloudless sky, so our modeled opacities should be viewed as lower limits. Δλ is the spectral bandpass of the Navcams (250 nm), P(Θ) is the scattering phase function of the clouds, F₀ is the in-band solar flux at the top of the atmosphere, and τ_COL is the atmospheric column density as measured by Mastcam imaging of the solar disk (Lemmon et al. 2024).

For an in-depth discussion of the uncertainties involved in the HC model, refer to Kloos et al. (2018). Here, we will note that the largest uncertainty comes from P(Θ), as the scattering phase function of Martian clouds is poorly constrained. Previously, Moores et al. (2015) and Kloos et al. (2016, 2018) assumed that this term was fixed at a constant value of 1/15. This assumption was justified as being representative of the flat phase function of terrestrial cirrus clouds at large scattering angles (Wang et al. 2014) and was verified by comparing the opacities derived from ZMs and SHMs taken less than 10 minutes apart. These movie pairs observe the same cloud at scattering angles up to 30° apart, and no significant difference in opacity was observed in all but three pairs (Kloos et al. 2018).

Despite the apparent validity of the flat phase function assumption in the scattering angle regimes explored in the literature to date, we will not be adopting this assumption because the MY 34–36 data set covers a wide range of scattering angles (16°–142°; see Figure 6). Previous studies using MSL data (Cooper et al. 2019; Innanen et al. 2024) have shown that the scattering phase function of ACB clouds over Gale deviates significantly from a flat curve at low scattering angles near the Sun, indicating that continued use of this assumption is not appropriate at low scattering angles. Kloos et al. (2016) determined a lower limit for the scattering phase function using the ZMs and SHMs but was limited to scattering angles between 70° and 115°, with the bulk of the observations between 80° and 90°. These values are outside the range where the value of the phase function is expected to deviate from the flat assumption. Their analysis was also limited by the availability of Mars Climate Sounder (MCS) cloud opacity retrievals to use as a reference, which cut the number of useful data points down to 12 ZMs and 8 SHMs.

Zoom In Zoom Out Reset image size

The second model, reproduced here as Equation (2), is known as the Whole Atmosphere (WA) model (see Equation (3) of Moores et al. 2010 for a full derivation). It assumes that the clouds are at lower altitudes and have similar scattering properties as the bulk atmosphere (i.e., are dust-based rather than ice-based):

$$\begin{eqnarray}&&{\rm{\Delta }}{\tau }_{\mathrm{WA}}=-\mathrm{ln}\left[1+a-\displaystyle \frac{a}{\exp (-{\tau }_{\mathrm{COL}})}\right].\end{eqnarray} \tag{ 2 }$$

Here, a is I_λ,VAR (as in the HC model) divided by the mean radiance of the frame.

Given the assumptions made by these models, it is most appropriate to use the HC model during the ACB season, when we are most likely to be observing high-altitude water-ice clouds, and the WA model at all other times of year, when we are most likely to be observing lower-altitude dust clouds. To determine the temporal extent of the ACB, we adopt the same approach as Moores et al. (2015) and Kloos et al. (2016, 2018). The Mars Color Imager (MARCI) instrument on board the Mars Reconnaissance Orbiter (MRO) produces weekly weather reports, which we monitor for the presence of ACB clouds. During the study period covered here, there was one extended gap in MARCI coverage, lasting from L_s ∼ 75° of MY 35 to L_s ∼ 125° of MY 36. This covers the end of the MY 35 ACB season through the beginning of the MY 36 ACB season. As we are thus unable to determine when the MY 35 ACB season ended and the MY 36 ACB season began using orbital data, we instead adopt the range of solar longitude values where ACB clouds are most typically observed: L_s = 40°–150° (Kloos et al. 2018). We note that in MYs 32, 34, and 35, ACB clouds were observed by MARCI as soon as L_s = 10°–30°, so setting the beginning of the MY 36 ACB season at L_s = 40° may be later than is appropriate.

Because we are most interested in ACB clouds, the remainder of this paper will be focused on the values derived using the HC model. Comparisons between ACB HC values and non-ACB WA values are restricted to Figure 12 and Table 1.

We can invert Equation (1) to solve for the scattering phase function rather than the opacity:

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

Rather than assuming a value for the scattering phase function, we now have to assume a value for the cloud opacity. Kloos et al. (2016) used this form of the HC equation along with water-ice optical depth retrievals from MCS and the ZMs and SHMs taken during the first MY of the MSL mission to determine a lower limit for the scattering phase function across a narrow range of scattering angles (70°–115°). MCS water-ice optical depths were also been used by Cooper et al. (2019) and Innanen et al. (2024) in combination with a dedicated observation known as the Phase Function Sky Survey (PFSS) to determine the value of the scattering phase function across a much wider range of scattering angles (20°–150°), though the PFSS produces increasingly unreliable results for scattering angles below approximately 45°.

Using optical depths from MCS is potentially problematic for several reasons. First, MCS has difficulty retrieving values near the surface due to increased dust opacity in the lower atmosphere, with the water-ice column depths often limited to altitudes above 30 km (Guzewich et al. 2017). The clouds we observe over Gale Crater are frequently lower than that (Campbell et al. 2020, 2021a), so MSL and MCS may not be observing the same clouds. Indeed, we know that MCS and MSL are not observing exactly the same clouds, as MCS measurements are not taken simultaneously with the ZMs/SHMs/PFSS. Second, MRO, to which MCS is mounted, is in a Sun-synchronous orbit, meaning that it observes each part of the surface at the same time each day. Kloos et al. (2018) suggested that cloud opacities over Gale vary diurnally, so the MCS data would not capture this detail. Furthermore, the footprint of MCS measurements is quite large, with typical nadir observations having a horizontal resolution of 1–6 km (Hayne et al. 2012) and typical limb observations having a much worse horizontal resolution of 200 km (McCleese et al. 2007).

Ideally, we would like the opacity input to Equation (3) to be derived from observations taken of the same clouds that we see from MSL. If it is indeed the case that the flat, constant phase function assumption holds at large scattering angles, then we can assume that the opacities calculated using Equation (1) in that regime are representative of opacities at all scattering angles. To account for potential diurnal or interannual variability in opacities, we can split the observations taken in the flat phase function regime by MY and time of day, then use an average of each as an input to Equation (3) to determine the phase function at all observed scattering angles.

To ensure a consistent methodology across all observations taken to date, we will need to apply our new phase function to the previously published data sets. Unfortunately, neither Moores et al. (2015) nor Kloos et al. (2016, 2018) documented the specific points selected for analysis in each movie. However, the instrument pointing is recorded in the Planetary Data System (PDS) label for each observation, so we can determine the scattering angle at the center of the frame. The wide FOV of the Navcams means that the actual scattering angle observed can differ by up to 30° from this central value. If we assume that the highest-contrast points in the MY 32 and 33 data occurred near the Sun and close to the edge of the Navcam frame, as was frequently the case in the MY 34–36 data, we can apply a ∼20° correction to the central scattering angle as a reasonable guess for what the true observed scattering angle was. Doing so, we find that the dependence of the opacity on scattering angle for MYs 32 and 33 (assuming a flat phase function) is nearly identical to that in MYs 33–36 (see Figure 7). In MY 33, nearly all of the high-opacity, low scattering angle observations were taken in the morning, suggesting that if a proper phase function correction was to be applied to these data, the Kloos et al. (2018) diurnal variability would likely be much smaller than what was reported, if it is present at all.

Zoom In Zoom Out Reset image size

In Figure 8, the problems resulting from the flat phase function assumption are immediately obvious. Typical water-ice phase functions peak strongly in the forward-scattering direction, meaning that they would have values higher than the previously assumed 1/15. Because the phase function appears in the denominator of Equation (1), having a smaller value for the phase function will result in a higher value for the opacity, leading to a nonphysical increase in opacities at scattering angles below approximately 60°.

Zoom In Zoom Out Reset image size

Based on the results seen in Figure 8, we decided to set the minimum cutoff for determining the "true" opacity at a scattering angle of 70°. The average of the values further from the Sun than this cutoff (after splitting by MY and time of day) were used as the opacity input to Equation (3) when deriving the phase function.

Our six derived phase functions are presented in Figures 9 and 10, divided by MY and time of day, respectively. The phase function curves and associated 95% confidence intervals were found using a sliding average with a window size of 5°. In general, the phase functions behave as we expect: flat at high scattering angles and increasing as the measurements approach the Sun. This is qualitatively similar to the phase functions derived by Cooper et al. (2019) and Innanen et al. (2024), as well as to observed phase functions of terrestrial cirrus clouds (Wang et al. 2014). The MY 34 AM and PM phase functions cover a similar range of scattering angles, but this is not the case for MYs 35 and 36, whose AM phase functions extend approximately 40° and 20° further from the Sun, respectively, than their PM counterparts. These unusually high scattering angle observations are a combination of SHMs that were taken in the very early morning (before 06:30 LTST) and those that had their azimuthal pointing adjusted to avoid looking at high-elevation local terrain features, consequently observing an area of the sky further from the Sun than usual.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There does not appear to be any significant interannual or diurnal variability in the phase functions, as the curves for all 3 yr lie within the 95% confidence interval of the other 2 yr. This is consistent with other work undertaken on MSL that found no significant variability in the phase function of water-ice clouds over Gale (Innanen et al. 2024).

Because we do not see any significant interannual or diurnal variability in the shape of the scattering phase function, we assume that all clouds observed during the MY 34–36 ACB seasons have the same phase function and thus adopt a single phase function for all years at all times of day.

Before our phase function can be used in the HC model, it must be normalized. Traditionally, phase functions are normalized to a value of either 4π or 1. However, because our observations do not cover the full range of scattering angles, it is not possible for us to normalize our phase function on its own. Instead, it is normalized to a set of analytically derived phase functions. For a full discussion of the normalization procedure used, see Section 4.4.

Cloud opacities for the MY 34–36 ACB seasons are presented in Figure 11 as a function of scattering angle, solar longitude, and LTST, along with the values derived using the flat phase function assumption. As expected, the magnitude of the difference is largest at small scattering angles and for ZMs taken close to the ZM keepout window.

Zoom In Zoom Out Reset image size

In Figure 12, we present all measured cloud opacities from MYs 31–36. In Table 1, we present the mean preferred cloud opacities for the ZMs and SHMs, as well as the L_s ranges over which each method was applied. The HC opacities during the ACB season are much higher than the WA opacities during the dusty season, which is consistent with the lower-quality values seen in Figure 5. In reality, the transition between the two regimes is unlikely to be as abrupt as it is portrayed here. However, because we do not have the ability to distinguish between water-ice and dust clouds, we cannot model a transitional period when both dust and water-ice clouds might be observed.

Zoom In Zoom Out Reset image size

The early-morning science blocks that were first implemented during MY 33 have continued through MY 36, so we now have four full MYs of good diurnal coverage from 06:00 to 18:00 LTST. As was the case in Kloos et al. (2018), the gap in diurnal coverage near 10:00 LTST due to rover uplink activities persists. In MY 36, an additional coverage gap emerged near 14:00 LTST, with many fewer observations taken near that time than in previous MYs due to operational constraints imposed by the rover's age.

In the left panels of Figure 13, we present the opacities as a function of LTST, averaged across half-hour bins. There does not appear to be any obvious diurnal variability. MYs 33 and 34 do have a larger number of high-opacity movies in the morning, but these observations are outliers; the bulk of the clouds observed in the morning during these 2 yr are similar to those observed in the afternoon. Consequently, we conclude that while high-opacity clouds may be more common in the morning, there is no systematic difference in the opacity of morning and afternoon clouds over Gale.

Zoom In Zoom Out Reset image size

These results suggest that the variability reported by Kloos et al. (2018) is likely the result of the scattering angle effect demonstrated in this paper rather than a true difference in opacity. Unusually, the highest-opacity outliers in MY 33 were almost all SHMs. Because the SHMs are pointed at a more oblique angle than the ZMs, there is a longer physical distance between the rover and the clouds being observed. Consequently, we are looking through more atmospheric dust, which tends to lower the sensitivity of the observation. This effect can be seen most acutely during the non-ACB season, when the ZMs regularly record higher opacities than the SHMs taken at the same time (note in Figure 12 how the ZM observations clearly lie above the SHM observations).

There is a slight drop in opacity near both sunrise and sunset, but this is an illumination effect rather than a true change in opacity. When the Sun is near the horizon and the sky is darker than it is with the Sun at high elevation angles, the contrast between the clouds and the background sky decreases, lowering the value of I_λ,VAR for a given opacity. Properly addressing cloud opacities at these times of day will likely require a 3D radiative transfer model, which is beyond the scope of this paper.

Following the analysis of Kloos et al. (2018), we averaged our opacities into bins of width 5° L_s to examine variability across the temporal extent of each ACB season (the right panel of Figure 13). There appears to be no meaningful variability during the MY 34 ACB season, while both the MY 35 and MY 36 seasons have a peak near L_s = 140°–160°. However, the MY 36 peak is likely due to there being a small number of data points in that bin, one of which is a high-opacity outlier. We do not observe the L_s = 50°–90° peak seen by Kloos et al. (2018) in MY 33, which was attributed to the rapid transport of water vapor to cloud-forming altitudes by the Hadley cell winds, whose velocities reach their maximum near this time (Wolff et al. 1999).

We assumed that the "true" cloud opacity clusters around the average value of the measurements made at scattering angles >70°, so there is some circular reasoning in then concluding that there is no variability. We did attempt to account for potential diurnal and interannual variability by deriving a separate phase function for each MY and time of day, but there is more that we can do to justify this assumption.

To do so, we repeated our analysis using just those points that lie in the regime where the flat phase function assumption is valid, without applying our derived phase function correction. The results of this analysis are shown in Figure 14. Due to the smaller number of data points (many bins only have one or two observations), it is harder to draw robust conclusions, but there is no obvious variability, either diurnally, intraseasonally, or interannually. As such, we conclude that the observed low variability is a real feature and not a product of our assumptions.

Zoom In Zoom Out Reset image size

The apparent lack of diurnal variability in the opacity of ACB clouds over Gale is unusual in the global context. There is substantial evidence suggesting that Martian cloud formation is driven at least in part by thermal tides (e.g., Hinson & Wilson 2004; Szantai et al. 2021), implying that we would expect to observe diurnal changes in opacity. Indeed, diurnal variability in Martian aerosols has been previously documented, e.g., over Tharsis (Madeleine et al. 2012) and Jezero (Patel et al. 2023; Smith et al. 2023). The most frequently observed pattern is an enhancement in the early morning, attributed to fogs that form overnight, and a smaller enhancement in the early afternoon, potentially due to increased dust lifting activity delivering more cloud condensation nuclei to cloud-forming elevations (Tamppari et al. 2003).

However, these patterns are not necessarily consistent across the entire planet. Recent analysis of ACB optical depths retrieved by the Emirates Mars Infrared Spectrometer (EMIRS) on board the Emirates Mars Mission (EMM) orbiter has shown that while the double-peaked pattern with maxima in the early morning and late afternoon with a midsol low are observed over much of the spatial extent of the ACB, the magnitude of the differences can vary significantly from location to location, with some areas showing little to no diurnal variability (Atwood et al. 2022). The largest differences (Δτ ∼ 0.1–0.2) are seen over Tharsis and the northern midlatitudes, while the southern edge of the ACB (i.e., where MSL is located) experiences less variability between morning and afternoon clouds, similar to what we have observed. It should be noted that the EMIRS data do still show decreased cloud opacities near local noon in the vicinity of Gale (albeit a smaller difference than elsewhere on the planet), which we do not observe.

To validate our opacity retrievals, we can compare them with those from orbital spacecraft. For this purpose, we use retrievals from MARCI on board MRO and from the Emirates Exploration Imager (EXI) on board the EMM Hope probe. The retrieval pipelines for these data sets are described in Wolff et al. (2019, 2022), respectively.

Because MRO is in a Sun-synchronous orbit, MARCI observes all points on the surface at approximately the same local time (∼15:00 LTST, though the exact time can vary significantly across a single MARCI observation due to the instrument's large FOV). As such, the MARCI retrievals cannot capture any diurnal variability, but they can still be used to assess intraseasonal and interannual variability.

Unlike MRO, EMM is in an elliptical supersynchronous orbit, which allows for observations of the surface to be carried out at different local times, albeit with varying resolutions as the spacecraft's altitude changes. This means that EXI optical depths can be used to validate our claim of low diurnal variability over Gale. However, we can only directly compare with MY 36, as EMM arrived at Mars shortly after the beginning of MY 36.

Both data sets are of sufficiently high resolution that it would be inappropriate to extract the pixel value immediately over the rover, particularly for the SHMs. Depending on the altitude of the clouds (15–30 km during the ACB season; Campbell et al. 2020), the clouds seen in the SHMs may be as far as ∼290 km away from the rover, well outside of Gale (see Figure 15 for an example). As such, selecting the cloud opacity retrieved from immediately over the rover may not be representative of the cloud opacities retrieved from our movies, as the two retrievals would come from clouds with significant spatial separation. To account for this, we adopt the average of all MARCI/EXI pixels within the Navcam FOV as our comparison data set.

Zoom In Zoom Out Reset image size

The results of these comparisons can be seen in Figures 16 (MARCI) and 17 (EXI). Our opacity retrievals closely track the orbital measurements; we particularly note the increase in opacities near the end of the ACB season (L_s = 150°) in MY 35 and the slight dip near L_s = 60° in MY 36.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our ground-based retrievals occasionally experience upward deviations from the orbital retrievals near the beginning and end of the ACB season. This may indicate increased atmospheric dust loading at these times of year, as we are not able to distinguish between water-ice and dust clouds, while the MARCI and EXI retrievals can. This is particularly noticeable at the end of the MY 36 ACB season, when an early-forming regional dust storm swept over Gale.

In MY 32 and MY 36, it appears that our derived opacities are slightly lower than those measured from orbit, while still tracking the general evolution of opacities over the ACB season. There are three possible explanations. First, our opacity model is sensitive to the phase function normalization. Because we do not have the ability to properly normalize our phase functions, adjusting the normalization term will increase or decrease the derived opacities. Second, it may indicate that the points we assumed to be cloudless when selecting the high- and low-radiance points were not actually cloudless. Rather, they may still contain an exceptionally thin cloud layer. This would artificially decrease the value of I_λ,VAR in our model and thus the derived opacity. However, it would be surprising for this effect to occur consistently over the entire duration of the ACB season, so the phase function normalization is a more plausible explanation. Finally, the orbital measurements themselves may be too high (see discussion below).

We should also consider the fact that EXI and MARCI optical depths are an imperfect data set with which to compare. First, both retrievals were performed in the UV at 320 nm, well outside of the spectral range of MSL's Navcams. Second, much like our retrievals, both had to make assumptions about the nature of the phase function, choosing to adopt the phase function of a 3 μm droxtal ice crystal. For a more in-depth comparison of our adopted phase function to those of various ice crystal geometries, see Section 4.4.

It is also possible that neither the MARCI nor the EXI retrievals are truly representative of the water-ice optical depths, particularly over Gale itself. The surface reflectance model of the retrieval pipeline includes a spatially dependent treatment of the Hapke w parameter (the surface's single-scattering albedo) to account for the bright and dark regions present across the planet. However, the w parameter data are averaged into 1° bins, with an additional 5 by 5pixel median filter applied to remove artifacts. Gale is just over 25 wide and thus will be covered by ∼9 pixels. Gale is partially filled with dark basaltic sand, particularly south of Aeolis Mons, and it is unlikely that the 1° bins are sufficient to capture the substantial surface albedo differences at smaller scales caused by these sand deposits.

In the MARCI/EXI data, this problem manifests itself as unusually high water-ice optical depths over Gale. The unexpectedly high values are concentrated primarily south of Aeolis Mons (where the majority of the sand is located) but not over the mountain itself (see Figure 18). It appears as if Gale is nearly constantly filled with a low-altitude fog, but this has never been observed either from orbit or by Curiosity (Gale is notably one of the few large topographic depressions on Mars that does not fill with fog in the early morning; Möhlmann et al. 2009). Consequently, the higher water-ice optical depths measured by MARCI and EXI over Gale may be driven in part by surface reflectance effects rather than an actual increased amount of water ice in the atmosphere.

Zoom In Zoom Out Reset image size

We can attempt to quantify the magnitude of this effect by comparing the opacities in the north-facing SHM FOV with those in the south-facing SHM FOV. If the dune fields (which lie to the south of MSL) are artificially inflating the orbital water-ice optical depths, we would expect that the south FOV average would generally be higher than the north FOV average.

This comparison over the full year can be seen in Figure 19. Indeed, it does appear that the south-facing average is frequently higher than the north-facing average, sometimes by as much as ∼50% of the total reported optical depth. However, it appears that this problem is much more pronounced outside of the ACB season, where nearly every measurement reported to contain water ice is consistent with this surface reflectance effect (notably, the difference during this period appears to be worst near L_s = 300°, when dusty season cloud opacities are highest). During the ACB season, the story is more complicated. Although many MARCI measurements do have a higher optical depth to the south of MSL, a large number have a higher optical depth to the north.

Zoom In Zoom Out Reset image size

Because our cloud movies point southward during the ACB season, it is important to keep in mind that this surface reflectance effect may be leading to an imperfect comparison between the ground-based and orbital optical depth retrievals. However, if the orbital measurements are indeed too high during the ACB season, we can account for this by increasing the value of our phase function normalization, which will lower the optical depth values to match without affecting the overall fit.

Determining the shape of the scattering phase function is useful not just for our opacity calculations. Because the shape is a consequence of the microphysical properties of the cloud scattering centers, it can reveal details about the physical nature of the ice crystals that make up the cloud. This is of critical importance to the development of Martian global climate models, which generally assume the presence of spherical ice crystals.

Cooper et al. (2019) compared the results of the MSL PFSS observation to analytically derived phase functions for seven randomly oriented ice crystal geometries frequently found in terrestrial cirrus clouds (Yang & Liou 1996; Yang et al. 2010). Based on the results of their analysis, they concluded that ACB clouds are unlikely to consist solely of any of the seven geometries investigated. Instead, their derived phase functions could be the result of a mix of ice crystal geometries or an exotic geometry not seen on Earth.

However, these terrestrial data sets do have one serious problem that may limit the validity of their comparison—they were computed for particle sizes of ∼50 μm. Although this particle size is appropriate for terrestrial applications, it is approximately 1 order of magnitude larger than the particle sizes that have been previously determined for the ACB (3–4 μm; Clancy et al. 2003). This difference is problematic because at the wavelengths to which MSL's Navcams are sensitive, it takes us from a regime where geometric refraction of light through the ice crystals is appropriate to one where Mie scattering must be used. The switch from geometric refraction to Mie scattering will have major consequences for the shape of the phase function, most notably suppressing the peaks at 22° and 46° that are responsible for the wide variety of optical scattering phenomena (e.g., halos and parhelia) that are observed frequently on Earth but have only been conclusively seen once on Mars (Lemmon et al. 2022).

We made use of a database of analytically derived phase functions covering a spectral range of 0.2–100 μm for 189 particle sizes between 2 and 10⁴ μm (Ping et al. 2013; Yang et al. 2013; Bi & Yang 2017). This database includes nine different ice crystal geometries (droxtals, hollow and solid bullet rosettes, plates, single columns, and aggregates consisting of five plates, ten plates, and eight columns) with three different surface roughness conditions (smooth, moderately roughened, and severely roughened). This database allowed us to more realistically approximate what MSL's Navcams would observe by averaging across multiple wavelengths and particle sizes.

To construct our phase functions for each geometry/roughness combination, we averaged together the relevant phase functions for all wavelengths within the Navcam bandpass and particle sizes of 1, 3, and 6 μm. The contribution of each phase function to the final curve is weighted by the relative responsivity of the Navcams at that wavelength. We found that after averaging over the Navcam bandpass, the surface roughness has little to no effect on the shape of the phase function at these particle sizes. Consequently, all figures here were generated using just the smooth geometries (no surface roughness).

For each geometry, the phase functions must be normalized for the results to be meaningful. Ordinarily, phase functions are normalized to either 1 or 4π over scattering angles between 0° and 180°. Since our derived phase functions do not cover this entire range, we instead normalized our phase functions by the median value of the comparison phase function. This normalization procedure was applied to all comparisons.

The ice crystal geometry comparisons are presented in Figure 20. At these particle sizes, the averaging process removes many of the identifiable features in the individual phase functions, leaving us with smooth, generally featureless curves. Because our derived phase functions are also fairly smooth and featureless, it is challenging to say whether or not we should prefer any one geometry over any other. Additionally, none of the comparison phase functions show the "flattening out" behavior that we observe for scattering angles greater than 60°. This difference may suggest that ACB clouds are not composed of nearly pure water-ice crystals but rather something with radically different light-scattering properties, like ice-rimed dust particles.

Zoom In Zoom Out Reset image size

We also made comparisons with other phase functions that have been derived for Martian clouds. These included water-ice cloud phase functions from Viking (Clancy & Lee 1991) and the Thermal Emission Spectrometer (TES) on board MGS (Clancy et al. 2003), as well as Viking dust phase functions (Clancy & Lee 1991). The six comparison phase functions are plotted against our MSL-derived phase functions in Figure 21, with the same normalization procedure applied as in Figure 20.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because we are observing aphelion clouds, we would expect to see the closest agreement with the TES type 2(A) and 2(B) aphelion clouds. Indeed, we can see reasonably good agreement, though the limited range of the TES phase functions makes a robust assessment of the fit challenging. It should also be noted that TES phase functions were derived from observations made in the thermal infrared, while MSL's Navcams operate in the visible. It is known that the optical properties of terrestrial cirrus clouds are quite different in the infrared versus the visible—for example, they are optically thin in the visible but comparatively optically thick in the infrared, leading to their being the only cloud type that creates a net warming effect on Earth's surface (Stephens & Webster 1981). As such, a comparison between the MSL and TES phase functions may not be one-to-one.

The midlatitude and polar cloud phase functions do also appear to follow our derived phase functions relatively closely, though, as with the ice crystal phase functions, they do not level off at large scattering angles in the same way that ours do. Once again, the fact that all of the phase functions are fairly smooth makes a meaningful comparison difficult as they lack specific, distinctive features whose presence or absence could be used as a point of reference.

Additionally, we compared with the phase functions from the PFSS that is also executed as part of the MSL cloud observation campaign. The PFSS consists of nine three-frame movies that form a dome around the rover, covering nearly the entire sky, with the exception of the region surrounding the Sun. A phase function is derived from this observation using Equation (3), so it faces the same question we did in determining which cloud opacities to use as an input. Ultimately, Innanen et al. (2024) opted to derive six phase functions for MYs 34–36: three using the opacities presented herein and an additional three using MCS opacities. The comparison between the PFSS and our derived phase functions for each year is presented in Figure 22.

We would expect to see a close agreement between our phase functions and those derived from the PFSS, as both were constructed from observations of the same clouds at the same location using the same instrument on board the same spacecraft, unlike the Viking and TES phase functions. Indeed, we do see generally good agreement between our phase functions and those derived from the PFSS using MSL opacities. The MCS/PFSS phase functions do not fit as well, as the forward-scattering peak begins about 20° further from the Sun. We previously noted the temporal and spatial resolution challenges of MCS data, which may be responsible for the difference.

In both cases, we cannot evaluate the fit over the full range of scattering angles. Due to its design, the PFSS can observe closer to the backward-scattering direction than the ZMs and SHMs do, but it also cannot observe the forward-scattering peak. Although the PFSS can, in principle, look as close as ∼20° from the Sun, the fact that each pointing uses only three frames (as opposed to the eight used by the ZMs and SHMs) means that in the 20°–40° range, the background signal from the sky becomes increasingly dominant over that from the clouds. This makes it more challenging to robustly extract the clouds from the sky using MFS, adding enough noise to the derived values that they are excluded from the final phase function. By stacking more images to create the mean frame in the ZMs and SHMs, we are able to improve the S/N of the observation, resulting in a more physically plausible phase function at low scattering angles.