Thermal Forcing of the Nocturnal Near Surface Environment by Martian Water Ice Clouds

We explore the potential role of clouds in moderating the nighttime temperature within Gale crater, as observed by the Rover Environmental Monitoring Station (REMS) instrument suite aboard the Curiosity rover. Just prior to aphelion, the decreasing trend in minimum daily temperature within Gale slows down. We investigate if this is due to increased formation of twilight and nighttime clouds, re‐radiating heat and reducing atmospheric and surface cooling. While diurnal analysis of REMS temperatures shows brief atmospheric warmings of 3–5 K post‐sunset during the occurrence of these clouds, an absence of similar warmings in the ground temperature measurements make it unlikely that clouds are the primary source. Seasonally, however, clouds that persist overnight could cause a warming of daily minimum surface temperatures in the Ls∼20°–50° season. This period can serve as a baseline and allow the potential effects of clouds to be more clearly discerned in the REMS temperature measurements. For this season, and in the peak of the aphelion cloud belt season, our modeled atmospheric energy budget shows a nocturnal decay signature of downward IR reflected and re‐emitted flux consistent with the presence and impact of clouds. The expected approximate exponential decay of this flux post‐sunset is damped more heavily in cloudier seasons than less cloudy or dusty seasons, suggesting formation and thickening of ice clouds as atmospheric nighttime temperatures cool.

clouds (likely water ice) have an in situ measurable effect on the surface temperature and near-surface air temperature within Gale crater (latitude ∼5°S), particularly in the nighttime.
Ground temperature sensor (GTS) and air temperature sensor (ATS) data from REMS show that the recorded minimum diurnal temperature over the course of a martian year does not always evolve in what might be considered a typical seasonal or diurnal pattern. To first order, nighttime temperatures are warmer around perihelion (occurring at solar longitude, L s = 251°), and cooler around aphelion (L s = 71°) but in Mars Years (MYs) 32 and 33, during the onset of the aphelion season (∼L s = 20°-50°) when the martian atmosphere is relatively clear Smith et al., 2001), a perturbation develops in the diurnal minimum temperature curve for both ground temperature (GT) and air temperature (AT). While MY 34 is also contained within the available data set, there are large periods with missing data making it less useful (as seen in Figure 1).
Here, we aim to test two separate hypotheses about these clouds on the local meteorology: (a) the presence of aphelion cloud belt (ACB) perturbations deviating from the first-order expected seasonal temperature trends (Figures 1b and 1c) are consistent with the effect of reflected and re-radiated IR flux from the increased presence of clouds, and (b) clouds contribute to shorter-timescale increases in twilight AT via reflection of outgoing solar flux.
We correlate REMS AT and GT measurements with the presence of twilight clouds observed by the REMS Ultraviolet Sensor (UVS) and augmented by nighttime Mars Climate Sounder (MCS) water ice opacity profiles. These all support the hypotheses that twilight and nighttime clouds are consistent with these seasonal nighttime temperature perturbations seen in the REMS data. Additionally, we probe the relative magnitude and seasonal evolution of nighttime atmospheric reflected and re-emitted IR flux, by using an energy balance model, and compare our results to the observed seasonal temperature and opacity trends. Section 2 provides an overview of the data utilized in this study, while Section 3 discusses the numerical methods we have employed for processing of the REMS data. Section 4 discusses our results, as constrained by the limits of the available data, and we present our conclusions in Section 5.

Observed Temperature Perturbations
Figures 1b and 1c shows the Ls∼20°-50° perturbations highlighted with blue arrows, and a sum of four-slowest intrinsic mode functions from an empirical mode decomposition (Huang & Wu, 2008), plotted atop the data, all of which is discussed further in Section 2.2. These perturbations temporarily reverse the otherwise secular, seasonal cooling trend as Mars approaches aphelion. They coincide with another slight change in the trend of MCAM opacities as seen in Figure 1a, the onset of the seasonal peak in water ice cloud opacity ( Figure 1a shows total opacity of all aerosols observed by MSL Mastcam, while water ice opacity derived from MSL Navcams and ChemCam is shown in Kloos et al., 2018 andMcConnochie et al., 2018 for MYs 31-33, respectively), including the equatorial ACB that typically extends from L s ∼50° to 150° (e.g., Clancy et al., 1996).
Specifically at twilight, high altitude clouds can reflect solar radiation, while further into the night, clouds can both reflect and absorb outgoing infrared radiation from the surface, re-radiating heat back towards the surface and reducing nighttime cooling (as is observed on Earth with high altitude cirrus; e.g., Yu et al., 2018). Around this seasonal perturbation, the REMS ATS detects increases in AT on minute-to-hour-long timescales during twilight, when the diurnal temperature is otherwise trending downwards. These measured AT perturbations occur on evenings when noctilucent clouds are observed at twilight by MSL cameras, and detected by the REMS UVS, and higher values in high-altitude water ice opacities are measured by MCS (Figure 1a for total opacity observed by MSL Mastcam; Kloos et al., 2018 for water ice opacity derived from MSL Navcams MYs 31-33; Giuranna et al., 2021).

Specifications for the REMS GT and AT
MSL REMS data are used for this analysis, covering the start of the mission through Sol 2,224 (MY 31 at solar longitude ∼155° (L s ) to MY 34 at L s ∼284°). REMS data records have been downloaded from the Planetary Data System (PDS) through the atmospheres node (data set ID: MSL-M-REMS-4-ENVRDR-V1.0; 3 of 13 Gómez-Elvira, 2013); the AT as measured by the REMS tip thermistor from sensor 1 is used. The REMS UVS data are the revision described by de Vicente-Retortillo et al., 2017 For a typical sol, five-minute REMS observations containing GT and AT measurements at a frequency of 1 Hz are made every hour with respect to local mean solar time (LMST), in addition to a varying number of extended-duration observations of at least one hour of continuous 1 Hz observations. Thus, the timing of REMS GT and AT measurements with respect to local true solar time (LTST) changes with season, following the equation of time over the Mars year. Ideally, this coverage would lead to a comprehensive data set over the mission; however, there are several identified limitations in the quality of the data. The REMS GT data, for example, is especially noisy during aphelion (electronic noise during the cold nighttime around aphelion dramatically decreases the precision of the GT data) and the first few minutes of each hour during this period present large uncertainties. To remove these noisy data for the present study, only GT values given a REMS GTS quality flag in the PDS data set (signaling an adequate power signal and sufficiently warm temperatures in the GTS electronics) are used. This quality control generally leads to additional temporal gaps in the nighttime data set, beyond the ∼55-min spacing between hourly observations. To reduce high-frequency electronic noise in the temperature data (as done in, e.g., Miller et al., 2018), a 45-s running average is used, providing a 1 K RMS error. Assuming that neither air nor surface temperatures change significantly from one second to the next, the difference between adjacent temperature measurements should be dominated by electronic noise. Thus, this difference is mostly a measure of electronic noise, and, from Miller  (2018), no significant improvement in reducing the noise is found beyond using a 45-s running average window. This running mean is applied to all continuous regions of the 5-min and hour(s)-long blocks of AT and GT data, buffered by 22 s on each end of the block to account for the halfwidth of the 45-s window. The running averages are then used in place of the 1-Hz measurements to investigate both the effect of nighttime and twilight clouds on AT over a minutes-to-hours timescale within a single sol, and for performing the seasonal and nocturnal energy balance analysis with AT and GT discussed in Section 3.
To assess the robustness of the aphelion temperature perturbations in the minimum daily REMS GT, we use a combination of the previously mentioned multisol running mean, a Singular Spectral Analysis (e.g., Ghil, 2002), and a Hilbert-Huang empirical mode decomposition (Huang & Wu, 2008). All of these show the pre-aphelion onset perturbation (L s ∼20°-50°) in the REMS minimum GT and AT, when the modes varying on the timescale of one to a few sols or faster are removed. This means retaining only the first three principal components of the Singular Spectrum Analysis or the first four intrinsic mode functions of the Hilbert-Huang. For this work, we show the filtering obtained with the Hilbert Huang modes in Figure 1, since they provide a lower standard deviation from the data than the Singular Spectrum Analysis. The Hilbert-Huang transform is a filter roughly based on the number of maxima, minima, and "mean crossings" or crossings through the mean value of a signal, as the signal goes from a minimum to a maximum (Huang & Wu, 2008). A decomposition into Intrinsic Mode Functions (IMFs) returns a slowest (last) mode that has one mean crossing, typically associated with, at most, one oscillation for the full temporal range of the analyzed data. Increasing (faster) Intrinsic Mode Functions (IMF) orders have increasing mean crossings with a decreasing time difference between them and, therefore, increasing oscillation frequencies.
The third-slowest mode decomposition presents a time difference between each two mean crossings similar to the multi-sol difference between annual temperature maxima and minima in the observed slower seasonal scale of the signal. The fourth-slowest IMF mode adds some subseasonal time scale oscillations while still on a multisol, almost seasonal, scale, and subsequently the sum of only the slowest four modes is shown atop the data in Figure 1.

Thermal Forcing in Twilight GT and AT Data
The effects of clouds should be seen in both AT and GT, as the atmosphere is warmed by water ice clouds via thermal re-emission from the ground. Higher thermal inertia of the martian surface compared to the atmosphere makes GT on Mars less easily perturbed than AT. As such, GT should more clearly reveal longer-duration temperature fluctuations caused by continuous nighttime cloud layers, compared to AT, which would capture more rapid fluctuations caused by winds and the ready exchange of air masses in the martian atmosphere.
The running averages of REMS AT data from 17:30 to 20:30 LTST are compared for sols when twilight clouds were, and were not, observed by the REMS UVS and Mastcam imagery in the pre-aphelion season. This threehour timespan was chosen so as to begin slightly before sunset, and to continue to 2.5 hr past sunset, when the Sun is well below the horizon, thus capturing potential effects of any noctilucent clouds (e.g., Donahue et al., 1972 for Earth;Lemmon et al., 2019or Lorenz et al., 2020 for Mars) present through the various stages of twilight (ending with the Sun 18° below the horizon), and into a totally dark sky.
To understand the amount of shortwave radiation that could be reflected from such a twilight cloud while the Sun is below the horizon, we use the Mars Weather Research and Forecasting (MarsWRF) general circulation model (Richardson, 2007) in its 1-D mode to simulate scattering from water ice clouds of various opacities, heights, thicknesses and particle sizes. The version of MarsWRF we use here, and below in Section 3.2, is nearly the same as that described in Toigo et al. (2012, although in its single-column implementation, and includes a two-stream, correlated-k radiative transfer scheme (Mischna et al., 2012), use of the Yonsei University boundary layer scheme for vertical mixing in the lower atmosphere (Hong et al., 2006), and a simple microphysics scheme that allows for condensation/sublimation of atmospheric CO 2 , and formation of water ice/CO 2 ice clouds when the atmosphere becomes saturated. The model also permits the user to prescribe cloud parameters, which is the approach we implement here.
Because MarsWRF always considers the Sun to be "above" the atmosphere, it is not possible to directly simulate incident solar radiation when the zenith angle >90°. To get around this limitation, and to understand the 10.1029/2020JE006737 5 of 13 maximum shortwave radiation that could be reflected by a cloud illuminated, and observed, from "below", we instead invert the model orientation, and calculate the shortwave flux reflected by a similar cloud illuminated, and observed, from above (this requires removing Rayleigh scattering and zeroing out the surface albedo to ensure our calculated flux derives only from the cloud itself, and not the atmosphere or surface). In the backwards hemisphere, the peak in scattering, and thus peak in signal, will occur around a phase angle of 180°. In this way, the model can inform us about the amount of reflected shortwave flux that would otherwise be attributable to these twilight clouds. The impact of this reflected flux on GT and AT can then be evaluated separately, as detailed in Section 4.1.

Thermal Forcing in Nighttime GT and AT Data
Given the trends observed in the seasonal REMS AT and GT measurements shown in Figures 1b and 1c, we want to assess if the nighttime (sunset to sunrise) flux due to clouds, (Wm −2 ), varies with respect to season in a manner consistent with the presence of nighttime clouds. The REMS UVS provides evidence of the presence of clouds after sunset, but not of their persistence through the night. To identify whether the presence of such clouds can be confirmed in the nighttime atmosphere, we calculate an energy balance at the martian surface, incorporating information both from numerical climate models and available REMS data. Because REMS does not contain instrumentation sufficient to solve for the full surface energy balance, we apply output from the MarsWRF GCM to represent the missing energy terms.
We assume nighttime conditions for our energy balance, which means we equate outgoing longwave radiation from the surface with all other terms contributing energy into the surface. These other terms include downwelling IR radiation from the overlying atmosphere, eddy conduction and sensible heating, subsurface conduction, and, potentially, contributions from atmospheric clouds. The energy balance, then, is shown in Equation 1: The term on the left hand side represents the surface IR emission. The time-varying values of are acquired from the running mean of the REMS GTS data. Surface emissivity is commonly assumed to range from between 0.9 and 1.0 (e.g., Christensen et al., 2004); here, we have selected a value of 0.9. The first term on the right-handside is eddy conduction between the surface and atmosphere. Both calculations from REMS data (e.g., Martínez et al., 2014) and numerical modeling show the eddy flux to be almost vanishingly small at night, generally <1 Wm −2 , and so largely negligible as compared to the other terms in Equation 1. For our energy balance, we assume a constant value of 2 Wm −2 , although, as we shall see, the choice of this value is not particularly important.
The second term on the right-hand side is subsurface conduction. This term will be negative during the daytime as the Sun warms the surface and the surface temperature is higher than in the first few cm of the subsurface (which corresponds with the diurnal penetration depth). At night, however, the surface cools, and falls below the subsurface temperature soon after sunset, yielding a positive value for this term. We use the running mean REMS GTS temperatures, again, for T sfc , but apply MarsWRF subsurface temperatures calculated for a depth at which the temperature remains constant on the diurnal timescale, which represents the temperature to which the surface value will relax. The depth of the thermal wave is equal to (KP/π) 1/2 , where K = k/(ρc) is the subsurface thermal diffusivity, and P is the period of thermal perturbation (86,400 s for the diurnal cycle). For a wide range of surface properties, the product of ρc ≈ 10 6 . Furthermore, k = I 2 /(ρc), where I is the thermal inertia of the surface material. With this information, we can both calculate the depth of the diurnal thermal wave at Gale crater (∼3 cm), and the overall subsurface term. Temperature output from MarsWRF is obtained at 5 cm depth every 10 sols for the full martian year. The subsurface thermal model of MarsWRF that calculates these temperatures is the same as described in Mischna and Piqueux (2020), and has been shown to match closely to more detailed thermal models (e.g., the KRC subsurface model, Kieffer., 2013) over the full diurnal and seasonal cycle, with departures of only a few K throughout the year, similar to that of other published models (e.g., Haberle et al., 1993;Mellon et al., 2000).
The third term on the right-hand side captures the downwelling IR from the atmosphere excluding the contribution from water ice clouds and dust. We note that this term is not measured by REMS. Here, too, we obtain the necessary downwelling IR flux from the MarsWRF GCM (i.e., the "sky temperature"), output every three hours for the full martian year.
The final term on the right-hand-side, , represents those fluxes contributed by the presence of aerosols such as water ice clouds or dust, in the atmosphere. There are two components to this contribution: first and most significant, reflected IR flux from outside of the 15 μm CO 2 band, and second, IR emission from the water ice clouds or dust, also from outside of the 15 μm CO 2 band. Clouds, when present, are sufficiently high and cold that this second component is significantly smaller than the first, and also relies on unknown values, such as the cloud emissivity (a function of cloud thickness and spatial coverage, which are, themselves, mostly unknown).
The fourth term, however, also serves as something of a residual "catch-all" for uncertainties in the components that make up the prior four flux terms in the energy balance. Uncertainty in surface emissivity, for example, will lead to some small fraction of the "true" flux being captured by this fourth term, while the majority is captured in the surface emission term. Model uncertainties are folded into residual fluxes in , too, along with the effects of atmospheric dust. For our purposes, however, we are more interested in the relative flux differences between those periods when clouds are present, and when they are not, and so absolute accuracy is not particularly significant; rather, it is the relative seasonal, nocturnal magnitude of in which we are interested. We can therefore proceed ahead, examining the mean value of across the seasons within the MY 31-34 data set. This is done by inputting our calculated or observed terms in Equation 1, and solving for . Recall that REMS GTS measurement blocks occur hourly for different extended durations, and are smoothed over a 45-s window, while the model-derived terms are calculated every three hours for atmospheric temperatures or every 10 sols (for subsurface temperatures). Values for are calculated for each running mean REMS GTS data point, using Equation 1 rearranged to solve for , with the most appropriate values used for the other terms.
The effect of clouds can be potentially determined in two ways. First, we can compare the shape of the mean value of throughout the seasons with the seasonal REMS minimum GTS, and Mastcam atmospheric opacity curves displayed in Figure 1. We can look at seasons where the presence of noctilucent twilight clouds was first observed before the ACB reaches its peak, the peak ACB, and where water ice clouds are less common, to see if the shape and relative magnitudes of are consistent with clouds impacting the nighttime radiation balance.
We can also use a second approach, which is to examine nighttime diurnal trends for cloudy, dusty, and relatively clear seasons. Relative differences in the magnitude and shape of between these seasons may reflect the influence of nighttime clouds (either from reflected IR or cloud emission), and separate it from the impact of dust aerosols, if we can ensure that all other residual terms are not appreciably variable between the two periods. Here, the eddy flux term is negligible at night, irrespective of season, based upon the analysis of measurements and numerical modeling Savijärvi et al., 2020).
Having separated the effect of aerosols from the clear atmosphere, atmospheric emissivity due to CO 2 will not change across seasons significantly, meaning that uncertainty introduced in the term will cancel out between seasons. Likewise, the chosen value for (0.9) is fixed across seasons, and there is no reason to believe it will change at Gale crater during the year.
The only term that will potentially appreciably change in a way unrelated to season is the subsurface term. Changes in surface properties (thermal inertia, albedo) along the rover traverse will change the local subsurface temperature, T sub . Our selection of a model value of T sub , at 5 cm depth eliminates most of the variability due to surface composition for most surfaces traversed by MSL during the aphelion seasons in MY 32-34, and should be fairly representative of subsurface temperatures regardless of surface property at any local point. In other words, our time-varying choice of subsurface temperature is an approximation of the temperature to which the surface will relax, regardless of location. Hence, even if a significant portion of is made up of measurement and model flux uncertainties, they should be mostly constant between the clear and cloudy seasons, with the only variable flux term being from the clouds or dust, themselves.

Effects of Twilight Clouds on AT
In Figure 2, rows a and b show REMS GT and AT measurements between 17:00 and 21:00 LTST for six sols in the aphelion onset season, when the signature of twilight clouds is observed with the REMS UVS (Gómez-Elvira, 2013; Lemmon et al., 2019), and, for comparison, six adjacent sols where it is not (rows c and d). MSL Navcam and Mastcam images show twilight clouds on 7 sols out of 9 attempts within the 0° < L s < 40° range 10.1029/2020JE006737 7 of 13 for MY 34 during sols 2,400-2,434, when the presence of these clouds was first directly observed by MSL (Curiosity's Three-Frame Mosaic of Clouds, 2019; Lemmon et al., 2019). Excess brightness in UVS data with the Sun 4-8° below the horizon is consistent with twilight cloud Mastcam and Navcam imagery, and with UVS data from that L s range were then analyzed for similar signatures in MY 32 and 33. A slightly extended time range Figure 2. Rover Environmental Monitoring Station (REMS) ground temperature (GT) and air temperature (AT) data and running means between 17:00 and 21:00 LTST for six sols where the signature of twilight clouds was detected by the REMS UV sensor (rows a and b) and for six adjacent sols where they were not confirmed (rows c and d), for comparison. A short-lived AT warming is present in all sols where the UVS had the signature of twilight clouds, but is less common on sols where it had no REMS signature; the effect on GT on all sols is negligible around twilight.
(starting at 17:00 instead of 18:00 LTST) is used in the plots in Figure 2, simply to show the diurnal trends in AT and GT just preceding the sampled time period (18:00-06:00 LTST). Note that the amount and timing of the data within this four-hour LTST window varies from sol to sol, as rover observations are planned with respect to local mean solar time, and the timing of the extended, hour-long observations changes for each sol due to rover power resources and REMS science team requests. From Figure 2, prominent increases, or jumps, of 3-5 K in the AT running mean temperatures (red curve) are observed for around 15 min after sunset (∼18:00 LTST) for each sol in rows a and b. For the GT running mean temperatures (green curve), a far more subtle (or absence of a) signal is observed over this same short timescale. The perturbations and spikes in AT around sunset (in rows a and b) are also observed on some sols adjacent to the ones where the REMS UV sensor suggests the presence of twilight clouds (rows c and d), but not all of them.
At later hours, nighttime (∼02:30 LTST) MCS ice opacity profiles taken near Gale crater show multiple cases, on multiple evenings, where peaks in ice opacity indicate clouds at high levels (∼30 km) for solar longitudes spanning 33. . This solar longitude range includes MSL sol 2,425, when the presence of noctilucent clouds was confirmed and observed by MSL Mastcam, and the altitudes of the ice opacity peaks are consistent with those where noctilucent clouds would be expected to form. These observations are consistent with the persistence of clouds through much of the night.
The coincidental increase in AT is interesting to note, but the almost negligible effect on the GT is inconsistent with twilight clouds producing this warming effect, as the surface would need to first absorb the reflected solar radiation, and then re-emit it in the infrared to be absorbed by the atmosphere. In that process, an increase in GT, similar to AT, should be present.
Additionally, from our 1-D model discussed in Section 3.1, we find that reflection from water ice clouds can only contribute a small fraction (∼21%-30%) of the solar flux necessary (∼14 Wm −2 ) to explain the full ∼5 K temperature difference at the surface seen in Figures 2a and 2b. Thin clouds, such as those observed during the aphelion onset in Kloos et al. (2018) with visible opacity 0.05, are capable of reflecting ∼3-5 Wm −2 across a range of particle sizes from 2 to 10 m (within the range of particle sizes found in Clancy et al., 2003 using TES;Guzewich & Smith, 2019using CRISM). More opaque clouds (visible opacity 0.15, suggested by Kloos et al., 2018 for Gale crater in this season) can reflect as much as ∼8-12 Wm −2 . This suggests that other effects beyond these tenuous clouds must contribute a significant fraction of the energy needed to explain the observed atmospheric temperature perturbations.
At sunset, GT starts to exponentially decay towards an asymptotic temperature, where it remains for the night. Downwelling IR heating remains slowly varying throughout the night. There is a possibility that advecting air parcels or katabatic winds flowing over the rover deck and warm RTG could lead to these sudden warmings; however, we observe no correlation between rover azimuthal orientation with respect to the upslope direction and warming within the sols shown in Figure 2. Additionally, winds after sunset are not driven by convection, and usually have a timescale of seconds, not minutes like the warmings in Figure 2. An exploration of the pressure signal shows the absence of eddy-like signatures similar to those associated with dust devils (Ellehøj et al., 2010).
Another potential source of the twilight warmings are orographic winds, which would be expected to repeat their pattern daily, but we do not observe this AT perturbation every sol (Figures 2c and 2d), nor is the same pattern of temporary heating observed over sequential sols. Additionally, downslope winds capable of causing sudden short-lived increases in AT do so via oscillations during the evening transitions (e.g., Lehner et al., 2015 and references therein), but the signal within Figures 2a and 2b does not always show oscillations, and occurs at a time where on Earth (Lehner et al., 2015) and on Mars (Miller et al., 2018) near surface temperatures typically show stable behavior.

Nighttime Energy Balance Results
To assess if the nighttime atmospheric reflected and re-emitted flux, , varies diurnally and seasonally in a manner consistent with the signatures expected from nocturnal clouds (e.g., Hinson & Wilson, 2004), we have generated two-dimensional histograms of the energy balance results with respect to LTST and L s , in Figure 3 (for the whole MY) and 4 (for particular seasons). . Data corresponding to all solar longitudes for all analyzed Mars years are plotted in panels a and b, against L s and LTST, respectively. In panel a, the Savitzky-Golay smoothed mean flux (using a window of 7 bins) for the corresponding L s bin, is plotted atop the histogram.
From Figure 3a, the seasonal atmospheric reflected and re-emitted flux evolves seasonally. The minimum flux of the smoothed mean (20 ± 20 Wm −2 ) occurs during a season that is not particularly cloudy nor dusty (cf. Figure 1a), centered on L s ∼250°.
values around this season may extend below 0 Wm −2 , but they are within two standard deviations of the mean, which is reasonable given the uncertainties around discussed in Section 3.2. The maximum of the smoothed mean flux (∼50 ± 20 Wm −2 ) occurs around the L s ∼20°-50° perturbation present in Figures 1a-1c second within L s ∼120°-150° during the ACB, coinciding with increased Mastcam daytime opacities ( Figure 1a). As expected, a perturbation in is also present in the peak of the dusty season, and overlaps the peak of the MY 34 dust storm around L s ∼200°. is shown with respect to L s in panel a with a Savitzky-Golay filtered mean plotted overtop in black (using a window of 7 bins). Panel b shows the flux over all L s with respect to LTST.
Since MSL is not stationary, we note that the rover's activities and locations during the L s ∼20°-50° season in MYs 32 and 33 properties could be responsible for part, or all of these observed temperature and flux perturbations. In MY 32 during the pre-aphelion season, the rover was quickly driving toward Mount Sharp, traversing many units with variable thermal properties. In MY 33, the perturbation corresponds to the duration over which MSL was on the "Stimson formation," a hard sandstone surface with high thermal inertia (Fraeman et al., 2016). This interpretation is supported by depressed daytime temperature ranges over the exact same interval and contributes to the effect observed in GT. By contrast, the slight flattening of the decay rate in daytime Mastcam opacity near the UVS observation of twilight clouds in this season is consistent with the presence and intervention of nighttime clouds. Figure 3b, compares the reflected and re-emitted flux over all L s , as a function of LTST. The reflected and re-emitted flux follows a decaying exponential trend just after sunset. As time progresses, the concentration of data points in the histogram shows a damping and deformation of the normal exponential decay that would occur in the absence of clouds and dust. To try and isolate the superimposed effects of clouds and dust in Figure 3b, similar plots were created for four 30° seasons of interest: L s = 20°-50° covering the pre-aphelion GT and AT perturbation (labeled the "ACB onset" season), L s = 120°-150° as the "ACB Peak" season (Figure 1a Mastcam opacities;McConnochie et al., 2018), L s = 190°-220° as the typical "peak dusty" season (Battalio & Wang, 2021;Guzewich et al., 2018), and L s = 230°-260° being the "decaying dust" season (Battalio & Wang, 2021;Giuranna et al., 2021), in Figures 4a-4d.
The range of 30° was selected to find a compromise in isolating seasonal aerosols, while reducing the effects of variability in REMS observation timing and duration, for the generation of mean value curves. The average nocturnal curve for each season contained data from MYs 31 to 34, was smoothed with a 9-bin Savitzky-Golay filter, and plotted in panel e with its respective standard deviation as a continuous shaded errorbar.
As the surface cools following sunset, the thermal flux to be reflected and re-emitted by clouds declines; however, the steady development and growth in both opacity and size of clouds can oppose this effect (Wilson et al., 2007), especially within the ACB season. The nocturnal cooling of the atmosphere creates conditions for ice saturation over a greater range of altitudes, thickening the clouds as the night progresses (Wilson et al., 2007), which dampens the exponentially decreasing trend in . Of course, it is not possible to entirely isolate cloud and dust activity, as sustained dust suspension at high altitudes can occur with dust storms. These suspended particles can act as cloud condensation nuclei, inciting cloud formation once the atmosphere is cooled to saturation, something that is easier to achieve at higher altitudes in warmer seasons (Guzewich et al., 2013). Furthermore, the "decaying dust" season we have selected occurs after the dustiest season, and it is possible that there are still some lasting effects from large-scale storms; however, this is our best possible approximation based upon the timing of the ACB (Clancy et al., 1996), dust activity (Battalio & Wang, 2021), and seasonal averages of aerosols (Giuranna et al., 2021) for a season with less impact from both types of aerosols.
The large range of fluxes in panels a-d can be attributed to the fact that these seasonal data contain multiple MYs, and that is a sort of "catch-all" for energy balance model errors. Despite this, there are interesting trends observed within the smoothed nocturnal means and their standard deviations in panel e.
The "peak dust" and "decaying dust" mean curves have similar shapes and nocturnal trends, with the "decaying dust" curve having a mean magnitude approximately 15-20 Wm −2 lower than the "peak dust" curve throughout the night. Both seasons follow an approximate exponential decay in to first order, which is consistent with an absence of significant dust injection to the atmosphere after sunset.
In contrast, the "ACB onset" and "ACB peak" smoothed average curves are very similar in trend throughout the night, with little deviation in shape or magnitude. The point density plot (Figure 4d) for the "ACB peak" is less noisy than the "ACB onset" season ( Figure 4a), which could suggest less influence from dust.
Comparing the shape of the smoothed mean "peak dust" and "decaying dust" nocturnal curves with the two ACB curves, we notice the two former mean curves decline approximately exponentially, rounding off slightly around ∼23:00-24:00 LTST, slowing the decline toward sunrise. In contrast, the two ACB mean curves approach a more pronounced rounding-off point closer to ∼19:30-20:00 LTST, with much greater damping in the decline of for the remainder of the night. The ACB curves end the night with mean values greater than the "peak dust" mean by ∼10 Wm −2 around sunrise. These results are consistent with relatively rapid cloud formation and growth at various altitudes following sunset in the colder aphelion season (Wilson et al., 2007), but minimal, or delayed, formation of clouds that may be limited to higher altitudes during the warmer dusty seasons (Guzewich et al., 2013).
Within the previously outlined limitations associated with our energy balance analysis and the REMS data, we are able to show perturbations in our resultant nocturnal reflected and re-emitted atmospheric flux that are consistent with the hypothesis that nighttime thermal perturbations detected in the REMS GT and AT are signatures of nocturnal clouds, though there remain other potential contributors to these thermal perturbations that prevent a definitive confirmation of the source of these features.

Conclusions
In this study, we show that while it is less likely that twilight clouds observed with MSL's Mastcam, REMS UVS, and MCS are driving short-lived temperature perturbations in REMS AT post-sunset in the aphelion onset season, there is support for the hypothesis that nocturnal clouds present in the leadup to and duration of the ACB season leave a measurable thermal signature on the surface and the near surface atmosphere of Gale crater (Figures 1, 3  and 4). Figure 2 shows sudden short-lived AT warmings during sols when twilight clouds leave a signature in the UV signal; similar warmings are absent in most sols without the cloud UV signal. These warmings occur when the atmosphere is typically in its calm diurnal period; however, they are not present to the same degree in GT, and so it is unlikely that clouds are directly responsible for these warming events. There are several other processes ranging from topographic downslope flows to horizontal transport of air and wind gusts that can otherwise cause this warming effect. On the sols we analyzed, the apparent thermal signature of wind gusts are observed on all sols, even those without a perturbation, but they tend to be long lasting warmings. Thermal oscillations are also reported, but they are repetitive instead of single peaks and observed at a different time of sol. It is interesting to note that these warmings occur more frequently on sols with twilight clouds. Figure 3 shows seasonal peaks in IR reflected and re-emitted flux , in the same seasons that GT and AT daily minimum temperatures show a thermal perturbation in Figure 1, in addition to the peak ACB season. A modeled energy balance at the martian surface shows a stronger and earlier damping of the first order exponential decay of in the aphelion onset and peak ACB seasons, compared to peak and decaying dusty seasons. Additionally, the cloudier season mean values at the end of the night are greater in magnitude than the "peak dust" mean by ∼10 Wm −2 . Results are consistent with ice clouds that form and thicken as the upper atmosphere cools throughout the night, slowing the escape of upwelling IR to space, especially during the coldest seasons of the MY on Gale. Clouds are seen in the daytime during the entire range of L s between the end of the seasonal dust storms at L s ∼0° and the beginning of storm activity around L s ∼150°; however, this analysis focused on the increase in cloud formation in the onset of the aphelion season, using it as a baseline. Future work can build upon what we have presented here to further assess the potential effects of thermal forcing of the surface and atmosphere by twilight and nighttime martian clouds.

Data Availability Statement
The REMS data used in this work is publicly available from the Planetary Data System Atmospheres Node (REMS: Mars Science Laboratory Rover Environmental Monitoring Station RDR Data V1.0, MSL-M-REMS-4-ENVEDR-V1.0, https://atmos.nmsu.edu/PDS/data/mslrem_1001/). The data shown in Figures 3 and 4 are available in an OSF Repository via DOI 10.17605/OSF.IO/WBM32 (Cooper et al., 2021). The MarsWRF model used has temperature and atmospheric parameters detailed in Appendix A of Mischna & Piqueux, 2020.