Long-term variations of Venus' 365-nm albedo observed by Venus Express, Akatsuki, MESSENGER, and Hubble Space Telescope

An unknown absorber near the cloud top level of Venus generates a broad absorption feature from the ultraviolet (UV) to visible, peaking around 360 nm, and therefore plays a critical role in the solar energy absorption. We present a quantitative study on the variability of the cloud albedo at 365 nm and its impact on Venus' solar heating rates based on an analysis of Venus Express and Akatsuki's UV images, and Hubble Space Telescope and MESSENGER's UV spectral data; in this analysis the calibration correction factor of the UV images of Venus Express (VMC) is updated relative to the Hubble and MESSENGER albedo measurements. Our results indicate that the 365-nm albedo varied by a factor of 2 from 2006 to 2017 over the entire planet, producing a 25-40% change in the low latitude solar heating rate according to our radiative transfer calculations. Thus, the cloud top level atmosphere should have experienced considerable solar heating variations over this period. Our global circulation model calculations show that this variable solar heating rate may explain the observed variations of zonal wind from 2006 to 2017. Overlaps in the timescale of the long-term UV albedo and the solar activity variations make it plausible that solar extreme UV intensity and cosmic-ray variations influenced the observed albedo trends. The albedo variations might also be linked with temporal variations of the upper cloud SO2 gas abundance, which affects the H2SO4-H2O aerosol formation.


INTRODUCTION
The Solar radiance is the principal energy source for the atmosphere of terrestrial planets, such as Earth, Venus, and Mars. Inhomogeneous solar radiance absorption drives atmospheric motions, from small scale convection to large scale global circulation. These motions distribute excess energy, and transport mass and momentum in the atmosphere. Temporal variation of absorbed solar radiance is therefore an important indication of possible changes in the atmosphere. The long-term monitoring of solar energy absorption is particularly useful in radiative energy balance calculations as a major input energy into a planetary system.
On the atmosphere of Venus, the maximum solar energy deposition occurs in the upper cloud layer (60-70 km) rather than at the surface as in the case of the Earth (Crisp 1986;Titov et al. 2012). The maximum solar energy absorption in the clouds is due to an unidentified absorber, hereafter "unknown absorber", which has been an unsolved question in Venus research regarding its identity. Venus global scale clouds and upper haze are mainly composed of H 2 SO 4 -H 2 O (Titov et al. 2012;Allen 1964;Mills et al. 2007) which has a small imaginary refractive index (n i = [1 − 9] × 10 −8 ) in the UV to visible range (Palmer & Williams 1975;Hummel et al. 1988). As a result, the H 2 SO 4 -H 2 O clouds and haze absorb almost none of the solar radiance in this spectral range, but are effective scatterers, making a strong contribution to the total solar radiance scattered back to space (which is ∼75% of the incident flux) (Titov et al. 2012). UV images of Venus, however, show distinct patterns caused by the unknown absorber. The absorption spectrum produced by the unknown absorber is observed to reach maximum absorption levels between 340 and 380 nm, and then decrease smoothly with increasing wavelength from 380 nm through the visible range (Barker et al. 1975;Pérez-Hoyos et al. 2018). Some studies indicate that an absorption tail exists at wavelengths shortward of 340 nm (Pérez-Hoyos et al. 2018), and in the 170-320 nm range (Marcq et al. 2011). According to data from descent probes the unknown absorber may be located in the upper clouds (Tomasko et al. 1980;Esposito 1980), and absorbs about half of the solar radiance deposited at the cloud top level, accounting for ∼3 K/Earth day of the global mean solar heating around 65 km altitude, when the total global mean solar heating is ∼6 K/day (Crisp 1986). Many candidates have been proposed for the unknown absorber, including OSSO, S 2 O, S x , FeCl 3 , and iron-bearing microorganism (Mills et al. 2007;Frandsen et al. 2016;Krasnopolsky 2017;Pérez-Hoyos et al. 2018;Limaye et al. 2018). However, none of these species satisfy both the spectral features produced by the unknown absorber, and the lifetime and simulated vertical profile required to fit the observations (Krasnopolsky 2018;Pérez-Hoyos et al. 2018).
The UV patterns reveal clear temporal and spatial variations; for example, the well known strong zonal winds, or 'super-rotation', that rotate around the globe in 4-5 days (Barker et al. 1975;Rossow et al. 1990;Kouyama et al. 2013), and transport the unknown absorber horizontally. In addition to this background wind, Venus' infamous 'Yfeature' is explained with atmospheric waves, resulting in a short-term periodicity of 4 to 5 days (Boyer & Camichel 1961;Del Genio & Rossow 1982Peralta et al. 2015;Kouyama et al. 2015;Imai et al. 2016). In the mean time, short and long-term variations are also reported (Khatuntsev et al. 2013;Marcq et al. 2013;Kouyama et al. 2013;Lee et al. 2015a), which are tracked by the distribution or abundance of the unknown absorber. Temporal variations of latitudinal 365-nm contrasts are closely linked with the SO 2 gas abundance above the cloud top level (Lee et al. 2015a), suggesting influences of chemical processes on sulfuric acid cloud aerosol formations in the UV contrast (Esposito & Travis 1982;Parkinson et al. 2015). Since Marcq et al. (2013) reported a general decline in the SO 2 gas abundance in the same periods where a decline in the long-term 365 nm cloud top albedo was observed, Lee et al. (2015a) proposed that the rate of H 2 SO 4 production, dependent on SO 2 photolysis, may be the principal mechanism supporting the observed long-term 365 nm albedo variation trends. However, definitive claims regarding these relationships were not made by Lee et al. (2015a) due to uncertainties at that time regarding the impact of the instrument degradation on the retrieved 365 nm cloud top albedo (Shalygina et al. 2015).
In this study, we report that the long-term variations of the UV reflectivity are a real phenomenon through a comparison of four space-based instruments: imagers on board Venus Express and Akatsuki, and spectrometers on board MESSENGER and Hubble Space Telescope (Section 2). We carefully checked the same phase angle diskresolved data (Section 3.1.2), and update the calibration correction factor of the UV data of Venus Express through a cross comparison with the UV spectra of MESSENGER and Hubble Space Telescope (Section 3.1.4). We evaluate the updated UV image data of Venus Express by a comparison with UV images of Akatsuki, showing a successful performance. We find common long-term variations in both of the disk-resolved albedo and the whole-disk albedo (Sections 4.1 and 4.2). The results are employed in our radiative transfer model calculations to understand possible solar heating variations (Sections 3.3 and 4.3). We discuss the significance of these results on the relationship with atmospheric winds at the cloud top level, including Venus global circulation model calculations, and possible reasons for the observed 365 nm albedo variations in Section 5.

DATA
Our 365-nm reflectivity analysis covers a decade of period, from 2006 to 2017, with only a one-year gap in 2015. The data were acquired from four instruments. The longest period of monitoring, 2006-2014, was covered with the Venus Monitoring Camera (VMC) on board Venus Express (Markiewicz et al. 2007). Two sets of 365-nm images were taken with the UV Imager (UVI) on board Akatsuki; one set of images was acquired from far distance in 2011 after the first Venus orbit insertion failure of Akatsuki (Nakamura et al. 2014), and the other set of images was taken from 2015 December to 2017 May, after the successful orbit insertion (Nakamura et al. 2016). Regular star observations of UVI have been conducted from 2010 for the calibration purpose. This revealed steady sensitivity over time (Yamazaki et al. 2018), and the mean error range in 2010-2017 is 18%. Near-UV spectra were taken with the Mercury Atmospheric and Surface Composition Spectrometer (MASCS) on board MErcury Surface, Space ENvironment, GEochemistry, and Ranging (MESSENGER) during its Venus flyby on 2007 June 5, using the VIS-VIRS channel (320-950 nm, spectral resolution 4.7 nm) (Pérez-Hoyos et al. 2018). An error of 5-10% is estimated at 340-390 nm from star observations during MESSENGER's cruise phase (Holsclaw et al. 2010). Other near-UV Venus spectra were taken with the Space Telescope Imaging Spectrograph (STIS) of Hubble Space Telescope with the G430L grating mode (290-570 nm, spectral resolution 0.54 nm) on 2011 January 2 (Jessup et al. 2015). 5% error is estimated for STIS measurements from based on regular standard star observations (Jessup et al. in press).
MASCS, STIS, and UVI data overlap observations with VMC data in 2007 and 2011. This is important for radiometric comparisons, as all of three performed star observations, and retrieved radiometric uncertainties independently. VMC data is highly uncertain in terms of radiometric calibration, as its star observations revealed 82% error (Titov et al. 2012;Shalygina et al. 2015). Cross calibration between VMC and the Visible and Infrared Thermal Imaging Spectrometer (VIRTIS) on board Venus Express was conducted, using simultaneous overlapped spectral range observations between VIRTIS and VMC (Titov et al. 2012;Shalygina et al. 2015). However, the absolute calibration of VIRTIS was not done at the time of these publications, resulting in questions about this cross calibration (Lee et al. 2015a). Comparisons in our study show that this VMC and VIRTIS cross calibration factor and the retrieved sensitivity degradation ratio of the 365-nm filter results in too large a difference in VMC data with respect to MASCS, STIS and UVI (Section 3.1.4). That will be discussed in detail (Section 3).
VMC images were manually selected in order to filter out those exhibiting artifacts, which could not be successfully corrected with additional flat field corrections (Titov et al. 2012). In addition, we selected VMC and UVI images having the dayside of Venus completely within the field of view. There are two UVI flat fields of the CCD matrix, and this study used the one measured on the ground before the launch of Akatsuki (Yamazaki et al. 2018). The same data were used for the star flux calibration (Yamazaki et al. 2018), and the mean calibration correction factor (β) is 1.525 in 2010-2017. We multiplied this calibration correction factor β by the measured radiance of UVI. The other UVI flat field data were made with the on-board diffuser, and the public data in DARTS and PDS are generated using the on-board diffuser flat field.
Except VMC, which mostly observed southern hemisphere (Titov et al. 2012), the other three instruments observed both northern and southern hemispheres. For disk-resolved data, our comparison is done only for southern hemisphere, keeping the consistency with VMC data (see Section 3.1.3). Since 365-nm images show persistent dark low latitudes and bright high latitudes (figure 2), we compare low (0-30 • ) and high (50-70 • ) latitudes separately for all images. Spectral data were taken using a narrow slit, and observed Venus across from local noon towards the terminator of a half illuminating phase, as shown in figure 1. All spectral data sets fall into low and middle (30-50 • ) latitudinal bins, and we use only the southern low latitude bin data in this study. For disk-integration (Section 3.2), we did not distinguish hemispheres, and used all valid pixels on Venus disk. Table 1 shows the configurations of four data sets used in this study.

METHODS
For disk-resolved images, we calculated a radiance factor (Hapke 2012) (Section 3.1.1), and then albedo, applying photometric correction (Section 3.1.2). Spectral data were convolved using the filter transmittance function of the 365-  nm channel of VMC, and we applied the same photometric correction as for the images. In order to take into account UVI's 2011 data, we calculated whole-disk albedo (Sromovsky et al. 2001;García Muñoz et al. 2014, 2017, without photometric correction due to the small apparent size of Venus (Section 3.2). Radiative transfer model calculations were performed using the model and gaseous database in Lee et al. (2015bLee et al. ( , 2016 (Section 3.3), to estimate the abundance of the unknown absorber that explains the observed 365-nm albedo, and to calculate solar heating rates.

Reflectivity
We converted the observed radiance to a radiance factor, r F , (Hapke 2012), the ratio of bidirectional reflectivity of a surface to the perfectly diffuse Lambertian surface illuminated normally. We calculate the average solar flux at 1 AU (Chance & Kurucz 2010) at each of UV filters of VMC and UVI, S (W m −2 µm −1 ), and the radiance factor as where R obs is the observed radiance (W m −2 sr −1 µm −1 ), β is a calibration correction factor of VMC or UVI, and d V is the distance of Venus to the Sun (in AU).

Photometric corrections
365-nm images show a combination of a smooth gradient from the sub-solar point to the terminator, and dark features owing to the presence of the unknown absorber. The smooth gradient depends on the incidence (i), emergence (e), and phase (α) angles, which can be described by a photometric law (disk function), D(µ, µ 0 , α), where µ = cos(e) and µ 0 = cos(i). We can separate albedo, A(α), and disk function, D(µ, µ 0 , α), from the radiance factor (r F ) that is derived from the measured radiance (Shkuratov et al. 2011), This albedo A, the equigonal albedo, depends on α (Shkuratov et al. 2011).
Previous studies showed that the Lambert and Lommel-Seeliger law (LLS) performs better in describing the gradient depending on geometric angles, compared to the Lambert law and the Minnaert laws (Lee et al. 2015a(Lee et al. , 2017. Therefore, in this study we adopted the LLS (Buratti & Veverka 1983;McEwen 1986), where k is a coefficient depending on α, and α is in [degree]. Eq. 4 is updated from Lee et al. (2017), using more images to find the mean condition along phase angle. Figure 2 shows albedo A at three 5-degree phase angle bins: 75-80, 80-85, and 85-90 • , from top to bottom. Left four columns are VMC images, having middle-latitudinal views, close to UVI's equatorial views on the right two columns. While UVI data have quality navigation using limb-fitting (Ogohara et al. 2017), VMC data do not. So we restricted VMC data satisfying i < 84 • , e < 81 • , and r F > 0.05 (Lee et al. 2015a). The last condition causes non-smooth terminator for VMC images in figure 2. As shown in these example images, we find no systematic tendency of albedo along local time, but temporal variations in the brightness and morphology. We also attempted to search for systematic variations in the albedo along longitude, as the surface topography may affect the 365-nm albedo, particularly over Aphrodite Terra (Bertaux et al. 2016), but the longitudinal coverage of the VMC data is not evenly distributed over time, and additionally this depends on phase angle selections. A detailed analysis along longitude, latitude and time, requires a different approach from the broad range average utilized in this study. Here we focus on temporal variations using data obtained over a broad range of longitudes. We derive the mean latitudinal albedo from disk-resolved VMC and UVI images, divided into two broad latitude bins: low (0-30 • ) and high (50-70 • ) latitudes. We derive a low latitudinal albedo also from the MASCS and STIS spectral data to complete the cross comparison with VMC data (Section 3.1.4).

Areas of disk-resolved albedo for a comparison
From the equatorial orbit, UVI images show that cloud top albedo and contrast patterns primarily displayed northsouth symmetry. However, as the example in figure 3 shows, we also observed cloud top albedo patterns that were asymmetric across the equator. This asymmetry was observed frequently in January-February 2017. Wind fields retrieved from cloud motions also detected the similar asymmetry between northern and southern hemispheres over the same period (Horinouchi et al. 2018). Thus, we restrict our disk-resolved data comparison only for southern hemisphere, which was observed by all of four instruments.

Cross-comparison of disk-resolved VMC, MASCS, STIS, and UVI data
The observed long-term decrease of albedo had been previously attributed to the sensitivity degradation of VMC by Shalygina et al. (2015). Following these authors, we used the 2.32 calibration correction factor (their Fig.12), which was retrieved from the comparison with VIRTIS-IR, and their value for the degradation ratio, k d (= −16.2 × 10 −5 Figure 2. Example 365-nm images of Venus. These present the cloud top albedo as observed by VMC and UVI between 2006 and 2017 when the observational phase angle was comparable to either STIS or MASCS Venus observations. In these example images, only Venus dayside a.m. or p.m. quadrant was visible within the camera field of view. First row shows images obtained in the 75-80 • phase angle range, second row images were obtained at 80-85 • , and third images were obtained at 85-90 • . Left four columns are images taken by VMC covering latitudes extending from 90 • S to ∼10 • S, and the right two columns are pole-to-pole dayside images obtained by UVI. Though the higher latitudes are brighter in the UVI observations than the equator, the intense polar hood brightening detected at 40-70 • S by VMC in 2006 has not been observed by UVI. All images are photometrically corrected (see Section 3.1.2), and share the same color bar on the left. orbit −1 , and 1 orbit of Venus Express equals 1 Earth day). We correct all VMC data to the value at an initial sensitivity condition (0 orbit number of Venus Express) using the below equations given in Shalygina et al. (2015). where t is time (orbit), B and B 0 are the corrected and observed radiance respectively (W m −2 sr −1 µm −1 ), and β is the 2.32 calibration correction factor. Temporal sensitivity degradation correction is given as where k d is the sensor degradation factor (orbit −1 ). We can get B(t = 0) as So the final form is (8) Figure 4 shows a comparison of low and high latitudinal albedo at the same phase angle bins, using the corrected VMC data with Eq. 8 (Shalygina et al. 2015). These corrected VMC data are significantly brighter than any of the independently calibrated MASCS, STIS, and UVI observations. The large difference between 2006 VMC and 2016 UVI is especially noticeable, while data in 2006 supposed to have least sensitivity degradation. We therefore discard this correction on VMC data due to the inconsistency with other calibrated MASCS, STIS, and UVI data.
Instead, we use the star calibration correction factor, 2.0±0.822, for the initial calibration correction factor (β) of the 365 nm channel of VMC (Titov et al. 2012;Shalygina et al. 2015). The large error of 82% results in the ambiguous definition of the absolute radiance. In this study, we improve the calibration of the data using the data points of MASCS in 2007 and STIS in 2011 as a reference to fit VMC values (Table 1). To limit uncertainties that may arise from the influence of the aerosol scattering along phase angle (α), we utilize only VMC data obtained at the same phase angle as either the MASCS (α = 85 − 90 • ) or STIS (α = 75 − 90 • ) observations. There are differences of 13 and 31 days from the closest VMC image at the time of the MASCS and STIS observations, respectively, at corresponding phase angles. To compensate for possible short-term fluctuations that may have occured during those time periods, we derive the 80-day mean albedo observed by VMC in the low-latitude bin at the two phase angle bins of MASCS and STIS. We then use the 80-day mean albedo to calculate ratios of A(MASCS)/A(VMC) at α = 85 − 90 • , and A(STIS)/A(VMC) at α = 75 − 80 • , where A is mean low-latitude albedo. Using this process, the ratios of 0.74 for A(MASCS)/A(VMC) and 0.84 for A(STIS)/A(VMC) were inferred. Differences in the ratios may result from differences in the latitudinal coverage of the MASCS and STIS observations (figure 1), and may also incorporate possible temporal variation of sensitivity of VMC. Even though the latter is possible, the decreasing trend of low latitude albedo between 2007 and 2011 changes from 29% to 24% in the low latitudinal polynomial fit (figure 5), which is yet a minor effect on the results of this study. Using the mean value of the ratios, 0.79, over all the VMC data, the VMC and UVI albedo retrievals become reasonably comparable (figure 5). Thus, we adopt this value and apply a new modified VMC calibration correction factor of 1.58 (β = 2.0 × 0.79) to the VMC data used in our study. As figure 5 shows, when the modified calibration correction factor is applied, the 365 nm albedo observed by VMC and UVI are reasonably aligned; those in 2008 are overlapped data in 2016 (UVI).

Whole-disk albedo
In order to evaluate robustness of the 1.58 modified VMC calibration correction factor, we employ 82 images taken with UVI in 2011 to compare with VMC. These UVI images were obtained after the first failure of the planned Venus orbit insertion of Akatsuki (Nakamura et al. 2014). In those images, the apparent size of Venus is a few pixels across, but sufficient signal-to-noise ratios were achieved. We calculated the whole-disk albedo, A whole−disk , which is a function of phase angle (α), following the equation below (Sromovsky et al. 2001), where d Venus is the Venus distance to the Sun (AU), F Venus is the measured disk-integrated Venus flux (W m −2 µm −1 ), Ω Venus is the solid angle of Venus as seen from spacecraft (sr), and S is the solar irradiance at 1 AU (W m −2 µm −1 ), which is calculated for either UVI or VMC, using each of transmittance functions. . The observed solid angle of Venus, Ω Venus , is calculated as where r Venus is the cloud top level radius of Venus (6052+70 km), and d V−sc is the distance of spacecraft from Venus (km). Observed Venus flux, F Venus , is calculated as where (x, y) is a location of pixel in image, Ω pix is a solid angle of one pixel of either VMC or UVI, R obs is radiance (W m −2 sr −1 µm −1 ) in the target area (r < r o ), and r is the distance of (x, y) from the Venus disk center (emission angle 0 • ). r o is the radius range in which the measured radiance are summed considering the point spread function of the instrument (5 pixels) that is the required radius of aperture photometry of UVI star flux analysis. The whole-disk albedo can be expressed as A g Φ(α), where A g is geometric albedo, and Φ(α) is the phase law of Venus, describing the disk-integrated scattering efficiency as a function of phase angle (García Muñoz et al. 2014, 2017.

Comparison of whole-disk albedo of VMC and UVI
We calculated the whole-disk albedo, A whole−disk (Eq. 9), of the 82 UVI images in 2011, and all of the disk-resolved VMC and UVI images. The latter is possible because our analysis is restricted to the images for which the observable dayside, defined by the observation elongation angle, is fully captured within the field of view of cameras (see figure 2). The 1.58 modified calibration correction factor is applied to the VMC data. Figure 6 shows the results of the whole-disk albedo versus phase angle. The 82 UVI images taken in February-March (circles) and in May 11-20 (triangles), are compared to the VMC images obtained contemporaneously. As a reference, the whole-disk albedo results derived from 2015-2017 UVI data and 2006-2007 VMC data are also included in the plot. The vertical bar of UVI data indicates the 18% error in the absolute radiance (Section 2). Fractions of disk illuminated by the Sun changes with the solar phase angle, from 100% at 0 • phase angle to 0% at 180 • phase angle, so there is a dominant decreasing trend as phase angle increases. At small phase angles, a local minimum related to glory is apparent (García Muñoz et al. 2014;Lee et al. 2017). The polynomial fit of UVI's 2015-2017 data shows the empirical phase function Φ(α) in the 40-110 • phase angle range (brown solid line, Eq. is shown in Table 2). We shift this phase function vertically to fit the maximum and minimum VMC whole-disk albedo at 75-80 • phase angle (cyan lines), assuming a change of geometric albedo A g over time whereas the scattering properties Φ are the same. The lower cyan line that fits the VMC whole-disk albedo in 2011 February-March, encompasses UVI data at the same period. This means that the 1.58 modified calibration correction factor for VMC data works well.
In addition, our analysis of the 2011 and 2015-2017 UVI whole-disk albedo at 0 • phase angle, which is a geometric albedo A g , increases from ∼0.33 in 2011 to ∼0.40 in 2015-2017. This result confirms the increasing albedo trends from 2011 to 2015-2017, observed at high and low latitudes based on the recalibrated disk-resolved 2011 VMC and diskresolved 2015-2017 UVI data (figure 5). The observation of an increasing albedo over these time periods completely opposes the behavior that would be produced by progressive long-term UV sensor degradation (Shalygina et al. 2015). This cannot be used to investigate the influence of surface topography on the 365-nm albedo because for each date a broad range of longitudes is included in the derivation of the albedo at small phase angles; this includes the 150-315 • E range in 2011 March, the 180-330 • E in 2016 May, and the 165-(360)-45 • E range in 2016 December-2017 January. Additionally, the 2011 UV data have insufficient spatial resolution to segregate specific longitude and latitude topographic regions.

Radiative transfer model
We use a one-dimensional line-by-line radiative transfer model (SHDOM, Evans 1998) in the 0-100 km altitude range and in the 2000-50000 cm −1 (= 0.2 − 5 µm) range to estimate the abundance of the unknown absorber that fits the observed 365-nm albedo, and to calculate the solar heating rate at low latitudes at the local noon time. The configurations are the same as used in Lee et al. (2015bLee et al. ( , 2016. CO 2 line parameters were taken from a combined HITEMP2010 (Rothman et al. 2010), and one developed by Wattson & Rothman (1992) and Pollack et al. (1993), as described in Lee et al. (2016). We included collision-induced CO 2 absorption in near-infrared, and that of H 2 O (Lee et al. 2016). Line parameters of other gases, N 2 , SO 2 , OCS, HCl, CO, HF, and H 2 S, were taken from HITRAN2012 (Rothman et al. 2013), and vertical profiles of gaseous abundances were taken from Titov et al. (2007). We included Rayleigh scattering (Pollack 1967;Hansen & Travis 1974), and UV range absorption cross sections of SO 2 (Wu et al. 2000). Microphysical properties of the cloud aerosols (mode 1, 2, 2, and 3) were taken from Zasova et al. (2007). We took the vertical structures of clouds' extinction coefficient from Crisp (1986), as shown in figure 7. For the unknown absorber, we assumed Crisp (1986)'s absorption coefficient (Q abs ) of the mode 1 particle (0.15-µm mean radius cloud particles) (Knollenberg & Hunten 1980;Kawabata et al. 1980;Wilquet et al. 2009;Luginin et al. 2016) in the 0.3-0.8 µm spectral range in the upper cloud layer (57-71 km). This assumed vertical location of the unknown absorber has been widely adopted in previous solar heating calculations (Crisp 1986;Lee et al. 2015b;Haus et al. 2016).

Temporal variations of low and high latitudinal albedo
The albedo A (Eq. 9) is phase angle α dependent due to strong backward and forward aerosol scatterings (Lee et al. 2015a;Shalygina et al. 2015;García Muñoz et al. 2014). Therefore, we restrict a comparison of A(α) to data obtained While direct comparison of the mean A should be done using data at a specific phase angle bin, often there are missing data over time due to regular changes in phase angles along the orbit of the spacecraft. In order to have a better temporal coverage of the mean A variations, we calculate the percent deviation of A from the 2016-2017 mean phase curve,Ā UVI (α), which is derived as a polynomial fit of UVI data in 2016-2017 in the 50-100 • phase angle range (figure 8). Table 2 shows the coefficients of the polynomial fit to the mean phase curves at low and high latitudes. Deviations from the UVI phase curve are defined as:   beginning of 2011 are confirmed by the overlap with the January 2011 STIS data. At high latitudes, periods of albedo decrease are less pronounced and appear to be shorter lived than those at low latitudes. This may be an indication of the combined influence of the unknown absorber abundance and the meridional circulation (Hadley circulation).
In particular, the latter would remove older aerosols downward below the cloud top level, following the descending branch at high latitudes (Imamura & Hashimoto 2001), and leaving behind only the bright newly formed aerosols that support the existence of Venus' bright high latitudes.  We apply Eq. 12 also for the whole-disk albedo A whole−disk and use the mean phase curve (figure 6 and Table 2) to get the percent deviations of A whole−disk from the 2016-2017 mean phase curve. Figure 10 shows the result as a function of time. The same albedo decreases inferred from the disk-resolved data from 2006 to 2011 are shown in this whole-disk albedo analysis. The independent UVI data obtained in early 2011 (dots with error bar) are well overlapped with those of VMC, including the same short-term albedo variations between March and May 2011. A darker 365-nm albedo in 2011 than 2006-2007 or 2016-2017 is a common feature in both the disk-resolved latitudinal and disk-integrated data, implying that the 365-nm darkening that occurred between those dates was a global phenomenon.

Solar heating variations
Because of the influence of the unknown absorber on solar heating rate near the cloud top level (Crisp 1986;Titov et al. 2007;Lee et al. 2015b) owing to the broad absorption spectrum from UV to visible (Crisp 1986;Pérez-Hoyos et al. 2018), the long-term 365 nm albedo variation we present here should have had a significant effect on Venus' solar heating rate. We use the radiative transfer model (Section 3.3) to determine a required abundance of the unknown absorber, which is assumed to be fixed to the mode 1 particles in the 57-71 km altitude range (Section 3.3), and to estimate the resulting solar heating rate changes at low latitudes. We incorporate the observed low latitudinal 365-nm albedo variations into our model using a scaling factor f that is multiplied by the assumed initial absorption coefficient Q abs taken from Crisp (1986) for the unknown absorber. f = 1.0 means the initial value. f > 1.0 means more abundant unknown absorber, which results in darker albedo, and vice versa. So we control only the single scattering albedo of the mode 1 particles, but not the size of the particles. In order to fit the observed A, we calculate 5 sets of emission (e) and incident (i) angles that satisfy α = 88 • : (e, i) = (28 • , 60 • ), (38 • , 50 • ), (48 • , 40 • ), (58 • , 30 • ), and (68 • , 20 • ), to simulate different combinations of e and i. The resulting radiance factors are corrected using the same photometric law which was applied to the calculation of the albedo (Section 3.1.2), and the mean value is used to find f values that match the observed albedo. Figure 11 shows the relationship of f and the mean value of calculated 365-nm albedo for the 88 • phase angle conditions ( Figure 5). The mean albedo observed at low latitudes is 0.33, requires f = 1.18. This is close to the initial value, f = 1.0, that results in a calculated albedo (A) = 0.35. The approximate maximum albedo observed at low latitudes, A = 0.40, requires f = 0.65, and the minimum albedo, A = 0.25, requires f = 2.51. We employ these f values in net solar flux profile calculations at 15 • latitude, which is the middle of the low latitude bin, at local noon time. Figure 12 shows the net solar flux divergence spectrum as a function of altitude z, −∇ · F net = −(dF net /dz), where F net is net solar flux. The strongest influence of the unknown absorber appears around 400 nm, where decreasing absorption and increasing solar irradiance overlap. Figure 13 shows the calculated local noon time solar heating rate at 15 • latitude, derived from Figure 12. This represents that the solar heating rate varied from −25% to +40% from the mean (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017), which is a significantly large range. The peak of solar heating rate is 36 K/(Earth)day for the mean albedo condition, 49 K/day for the minimum albedo, and 27 K/day for the maximum albedo at this local noon time.
We note that the vertical structure of the upper clouds may affect the solar heating rate (Lee et al. 2015b), but the structure of the upper clouds at low latitudes from near infrared observations is shown to be rather stable during Figure 11. Expected albedo as a function of f . Vertical error bars correspond to uncertainties depending on incidence and emission angles (see text for details).
the time of Venus Express (Ignatiev et al. 2009;Cottini et al. 2015;Fedorova et al. 2016), validating the use of a fixed cloud structure. Other analyses suggest that the vertical distribution of the unknown absorber may sometimes extend vertically above the cloud top level (Molaverdikhani et al. 2012;Lee et al. 2015a). However, the upper haze and vertical locations of the unknown absorber are not changed in our calculations as these are beyond the scope of this study. Further sensitivity studies on solar heating rate would explore in depth the effects of variable vertical and latitudinal distributions of the unknown absorber on Venus' global solar heating rate.

Relationship between solar heating and zonal winds at the cloud top level
The observed 365-nm albedo variations should result in solar heating variations near the cloud top level as shown in Section 4.3. We note that two long-term trends occurred in parallel; in the period when the 365 nm albedo had declined, leading to increases in solar heating, the long-term cloud top zonal wind was observed to increase from 80-90 m/s to ∼110 m/s between 2007 and 2012 around local noon at 20 • S (Khatuntsev et al. 2013;Kouyama et al. 2013;Hueso et al. 2015). Likewise, between 2014 and 2016, when the 365 nm albedo had increased, leading to decreases in solar heating, the zonal wind speed was observed to slow down from 110 to 100 m/s (Horinouchi et al. 2018). This wind speed variation is qualitatively consistent with the expected change in vertical shear associated with cyclostrophic balance. As the strongest solar heating occurs at low latitudes, the low latitudinal heating can alter the pole-to-equator gradient of temperature. The increased low latitudinal solar heating inferred from the 2011-2013 period should have increased the meridional temperature gradient, which then increases the equilibrium vertical shear, leading to increase of wind speed. Additionally there would be contributions of change of momentum flux associated with the thermal tide, which has been expected to participate in the angular momentum budget. To test these ideas, a simulation was  run with the latest version of the IPSL Venus GCM (Garate-Lopez & Lebonnois 2018). Starting from the reference simulation (see Appendix A for details), the solar heating rate was modified as a function of time, slowly reducing it (the whole profile at the same time, a simple test adjustment) over 20 Venusian Solar days, i.e. 2340 Earth days or 6.4 Earth years. The rate of solar heating change is −40% in 6.4 years (from 7.5 K/day to 4.5 K/day in zonal mean solar heating rate), consistent with the evolution proposed in this study. The time evolution of the zonally-averaged solar heating rate in the latitudinal band between 10 • S and 20 • S and at 30 mbar (∼70 km, near the cloud top level) is plotted in figure 14, together with the simultaneous evolution of the temperature and zonal wind in the same region. There is a clear correlation between these variables with an amplitude only slightly higher than the observed zonal wind variations. This correlation is caused by a reduced meridional circulation that directly affects the transport of angular momentum upward and poleward, resulting in a reduction of the cloud-top zonal wind peak. Detailed analysis is described in Appendix A.

Possible causes of the observed albedo variations
There are many intervening agents which may act in combination to produce the observed albedo variations, for example, the chemical composition and reaction rate of the unknown absorber, its interaction with or dependency on the chemical state of other atmospheric constituents, and the variability of the cloud and haze structure as a function of time. We note that the SO 2 abundance above the cloud top level was observed to decline from 2007 to 2012 (Marcq et al. 2013), and then subsequently increase around 2016 (Encrenaz et al. 2019). Since the 365 nm albedo will become relatively brighter with more abundant pure sulfuric acid aerosols which are formed through photolysis of SO 2 , it is plausible that the observed albedo trend is linked to the long-term trend of SO 2 abundance.
The period of low 365 nm albedo in this study overlaps the known maximum time of solar activity cycle (Jiang et al. 2018) as shown in figure 15. This resembles the correlation of Neptune's reflectivity and the solar activity cycle (Aplin & Harrison 2016). Since solar EUV radiation might affect photochemical reactions involving SO 2 that are necessary for aerosol formation on Venus (Mills et al. 2007;Parkinson et al. 2015), further study is required to explore influences of solar activity cycle on the Venusian atmosphere. It is also possible that production rate of sulfuric acid aerosols is altered by galactic cosmic rays via ion-induced nucleation. Electrostatic interactions between ionized acid molecules can enhance new aerosol formations by reducing critical size and by increasing collision possibility (Lovejoy et al. 2004;Kirkby 2007), and there are observations on H 2 SO 4 -H 2 O ultrafine aerosols of less than 9 nm in diameter in Earths upper troposphere and stratosphere that were explained by the ion-induced nucleation (Lee et al. 2003). The peak of ion production rate in the Venusian atmosphere due to galactic cosmic rays is predicted at 62.5 km (Nordheim et al. 2015) with the 46-58 ion pairs cm −3 s −1 range of variations between solar minimum to maximum. So the upper haze aerosol formation might be effectively triggered by such ion-induced nucleation, and vary following the solar activities. Figure 15 shows comparisons of the 365-nm albedo at low latitudes, neutron cosmic ray detected from the Figure 14. Temporal variations of zonal mean solar heating rate (orange), wind speed (red), and temperature (blue) at 10 • S-20 • S at 30 mbar (∼70 km). The solar heating rate is controlled to decrease smoothly along time by 40% from the reference condition in the IPSL-Venus GCM (Appendix A). Simultaneous variations in temperature and zonal wind speed are shown together.
Oulu station 1 , and Lyman-α flux at the Earth location 2 . A 5-day mean is applied to compare the 365-nm albedo and cosmic ray to compensate the atmospheric rotation rate on Venus. The solar rotation rate (25-day) mean is applied for the comparison of the 365-nm albedo and Lyman-α flux, to take into account the different locations of Venus and of the Earth respect to the Sun. The results shows that the 365-nm albedo A has a negative correlation with Lyman-α flux, and a positive correlation with neutron cosmic rays, but these may act together with the mesospheric SO 2 gas influences.

Comparison with other planets
Regardless of the cause of the observed albedo changes, the range of albedo variation on Venus is surprising. On the Earth, clouds play a considerable role as a buffer of possible climate variations and are also a regulator of the solar energy distribution (Stephens et al. 2015). However, the clouds on Venus are different; rather than supporting a stable solar heating rate, drastic variations of solar heating seem to occur as inferred from the 365-nm albedo. The astounding nature of the albedo variation results we present here is further emphasized by results derived from other planetary albedo studies in the Solar System where weaker long-term albedo variations were observed. For example, at Neptune the observed magnitude varied by ±0.02 (corresponding to ±2% changes in flux) at the blue (472±10 nm) and red (551±10 nm) filters over 1972-2014 (Aplin & Harrison 2016), and at Mars where the surface albedo varied by 10% at the red filters (575-675 and 550-700 nm) from 1976-1980to 1999-2003(Geissler 2005.

Further studies
In addition to the impact of the solar heating rates at the cloud top level, the vertical profile of solar flux on Venus, down to the surface, should also be altered by the observed cloud top albedo changes. Such solar flux variations Figure 15. Comparison of the 365-nm albedo A, neutron cosmic rays, and Lyman-α flux density evolution with time. On each date the 5-day mean of low latitudinal A at 85-90 • phase angle bin (a); the 5-day mean of cosmic ray (neutron) measured at the Oulu station (b); the 25-day mean of Lyman-α flux (c) is shown. A comparison between (a) and (b) is shown in (d). A comparison between (a) and (c) is shown in (e), but using the 25-day mean for consistency. See text for details. may explain unbalanced net radiative energy below the clouds (Lee et al. 2017), and therefore further investigation is needed to understand the true impact of the albedo changes on the entire lower atmosphere of Venus.
Our study focuses on the observed 365-nm albedo and its direct impacts on solar heating at the equator. This does not cover detailed modeling of net radiative forcing, such as cooling rate changes (Haus et al. 2017), that would impact on microphysical and photochemical processes. Such studies must be completed to accurately infer the true impact of the solar heating on cloud formation and climate. The work we present here provides a foundation for future in-depth studies of links between Venus 365 nm albedo and processes that directly impact Venus climate. . This trend is consistent among four independent UV instruments; VMC and MASCS in 2007, and VMC, STIS, and UVI in 2011, either using disk-resolved or disk-integrated data. We discard the previously suggested sensitivity degradation of VMC (Shalygina et al. 2015), and propose a new calibration correction factor for VMC in this study. The ranges of albedo variations are ∼0.2-0.4 at low latitudes, and ∼0.3-0.6 at high latitudes in 2006-2017, so albedo had been varied by a factor of two over the last decade. The whole-disk albedo also shows a similar trend, changed from −30% to +20% compared to the mean value in 2016-2017, meaning that the albedo variation occurred on a global scale.
Our one dimensional line-by-line radiative transfer model calculations reveal this level of albedo variation can alter solar heating rate from −25% to +40% due to the broad absorption spectrum of the unknown absorber from UV to visible range. We suggest that this solar heating rate variation can be a cause of the observed long-term zonal wind speed variation at low latitudes that increased from 80-90 m/s in 2007 to ∼110 m/s in 2012, and then decreased to 100 m/s in 2016-2017. Wind speed increased during the low albedo time, when solar heating was stronger than average, implying that increased solar heating may play a role in wind speed changes through cyclostrophic balance, enhanced thermal tide, or vertical momentum transport. We show the results of Venus GCM simulations, which support the linear relationship between solar heating rate and zonal wind speed.
The observed 365-nm albedo variations might be caused by variations of SO 2 gas abundance above the clouds. We also suggest links between the 365 nm albedo, and the Solar Cycle and consequent galactic cosmic ray density variations. Continuous 365 nm observations are necessary to clarify the mechanism of the 365 nm albedo variations.  (Haus et al. 2014). The authors also used solar heating look-up table (Haus et al. 2015) according to the latitudinal cloud structures. Thermal cooling is based on the net-exchange rate (NER) formalism (Eymet et al. 2009) with additional continua. This last version of the IPSL Venus GCM was able to simulate prominent cold band surrounding the poles of Venus close to observations (Garate-Lopez & Lebonnois 2018).
In this latest IPSL Venus GCM, only solar heating has been reduced, mimicking the expected solar heating rate variations shown in this manuscript. Figure 14 presents the linear correlation among solar heating, temperature, and zonal wind speed. To interpret these correlations, the time variations of the latitudinal profiles of zonally and temporally (over 2 Venusian solar days) averaged heating rates (solar, infrared and dynamical terms of the energy budget), temperature, and vertical and zonal winds are plotted in figure 16 at 30 mbar.
Looking at the different heating rates ( figure 16a), it appears that in average, the decrease of the solar heating is mostly compensated by a decrease of the infrared cooling, corresponding to the decrease in temperature (figure 16b). At mid-to high-latitude, though, the dynamical term associated with averaged meridional and vertical motions is not negligible, and therefore the decrease of the solar heating is compensated by an impact on the averaged meridional and vertical winds. A reduction of the amplitude of the vertical wind is seen in figure 16c, except between 30 and 50 • of latitude, which may indicate some impact here of changes in the meridional energy budget. This reduction of the mean meridional circulation has a direct impact on the transport of angular momentum upward and poleward, inducing a reduction of the cloud-top zonal wind peak (figure 16d). Regardless of the simplicity of the simulation set-up, the results explain how solar heating variations can affect zonal winds. It will be important to compare with long-term temperature trend analysis that may be available in near future using night side temperature field retrieved from VIRTIS-H/Venus Express (Migliorini et al. 2012) or radio occultation temperature profiles from VeRa/Venus Express (Tellmann et al. 2012).