---
author:
- Yeon Joo Lee
- 'Kandis-Lea Jessup'
- 'Santiago Perez-Hoyos'
- 'Dmitrij V. Titov'
- Sebastien Lebonnois
- Javier Peralta
- Takeshi Horinouchi
- Takeshi Imamura
- Sanjay Limaye
- Emmanuel Marcq
- Masahiro Takagi
- Atsushi Yamazaki
- Manabu Yamada
- Shigeto Watanabe
- 'Shin-ya Murakami'
- Kazunori Ogohara
- 'William M. McClintock'
- Gregory Holsclaw
- Anthony Roman
bibliography:
- 'paper.bib'
title: 'Long-term variations of Venus’ 365-nm albedo observed by Venus Express, Akatsuki, MESSENGER, and Hubble Space Telescope'
---

Introduction {#sec:intro}
============

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 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 [@Crisp86; @Titov12]. 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 [@Titov12; @Allen64; @Mills07] which has a small imaginary refractive index ($n_i=[1-9]\times10^{-8}$) in the UV to visible range [@Palmer75; @Hummel88]. 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 $\sim$75% of the incident flux) [@Titov12]. 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 [@Barker75; @Perezhoyos18]. Some studies indicate that an absorption tail exists at wavelengths shortward of 340 nm [@Perezhoyos18], and in the 170–320 nm range [@Marcq11]. According to data from descent probes the unknown absorber may be located in the upper clouds [@Tomasko80a; @Esposito80], and absorbs about half of the solar radiance deposited at the cloud top level, accounting for $\sim$3 K/Earth day of the global mean solar heating around 65 km altitude, when the total global mean solar heating is $\sim$6 K/day [@Crisp86]. Many candidates have been proposed including OSSO, S$_2$O, S$_x$, FeCl$_3$, and iron-bearing microorganism [@Mills07; @Frandsen16; @Krasnopolsky17; @Perezhoyos18; @Limaye18]. 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 [@Krasnopolsky18; @Perezhoyos18].

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 [@Barker75; @Rossow90; @Kouyama13], and transport the unknown absorber horizontally. In addition to this background wind, Venus’ infamous ‘Y-feature’ is explained with atmospheric waves, resulting in a short-term periodicity of 4 to 5 days [@Boyer61; @DelGenio82; @DelGenio90; @Peralta15; @Kouyama15; @Imai16]. In the mean time, short and long-term variations are also reported [@Khatuntsev13; @Marcq13; @Kouyama13; @Lee15a], 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 [@Lee15a], suggesting influences of chemical processes on sulfuric acid cloud aerosol formations in the UV contrast [@Esposito82; @Parkinson15]. Since @Marcq13 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, @Lee15a 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 @Lee15a due to uncertainties at that time regarding the impact of the instrument degradation on the retrieved 365 nm cloud top albedo [@Shalygina15].

In this study, we report that the long-term variations of the UV reflectivity are a real through a comparison of four space-based instruments. We carefully checked the same phase angle disk-resolved data . The results are employed in our radiative transfer model calculations to understand possible solar heating variations (Sections \[subsec:RTM\] and \[subsec:RTMresults\]). We discuss the significance of these results on the relationship with atmospheric winds at the cloud top level, and possible reasons for the observed 365 nm albedo variations in Section \[sec:discussion\].

Data {#sec: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 with Venus Monitoring Camera (VMC) on board Venus Express [@Markiewicz07a]. Two sets of 365-nm images were taken with 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 [@Nakamura14], and the other set of images was taken from 2015 December to 2017 May, after the successful orbit insertion [@Nakamura16]. Regular star observations of UVI have been conducted from 2010 for the calibration purpose. This revealed steady sensitivity over time [@Yamazaki18], and the mean error range in 2010–2017 is 18%. Near-UV spectra were taken with 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) [@Perezhoyos18]. An error of 5-10% is estimated at 340–390 nm from star observations during MESSENGER’s cruise phase [@Holsclaw10]. near-UV Venus spectra were taken with 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 [@Jessup15]. 5% error is estimated for STIS measurements from based on regular standard star observations 2[[@Jessup19]]{}.

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 [@Titov12; @Shalygina15]. Cross calibration between VMC and Visible and Infrared Thermal Imaging Spectrometer (VIRTIS) on board Venus Express was conducted, using simultaneous overlapped spectral range observations between VIRTIS and VMC [@Titov12; @Shalygina15]. However, the absolute calibration of VIRTIS was not done at the time of these publications, resulting this cross calibration [@Lee15a]. 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 in VMC data with respect to MASCS, STIS and UVI (Section \[subsubsec:CrossComp\_Disk-resA\]). That will be discussed in detail (Section \[sec:methods\]).

VMC images were manually selected in order to filter out those exhibiting artifacts, which could not be successfully corrected with additional flat field corrections [@Titov12]. 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 [@Yamazaki18]. The same data were used for the star flux calibration [@Yamazaki18], and the mean calibration correction factor ($\beta$) is 1.525 in 2010-2017. We multiplied this calibration correction factor $\beta$ 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 [@Titov12], 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 \[subsubsec:disk-res\]). Since 365-nm images show persistent dark low latitudes and bright high latitudes (figure \[fig:img\_comp\]), we compare low (0–30$^{\circ}$) and high (50–70$^{\circ}$) 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 \[fig:globe\_spec\]. All spectral data sets fall into low and middle (30–50$^{\circ}$) latitudinal bins, and we use only the southern low latitude bin data in this study. For disk-integration (Section \[subsec:whole-diskA\]), we did not distinguish hemispheres, and used all valid pixels on Venus disk. Table \[tab:data\] shows the configurations of four data sets used in this study.

  --------------------------------- --------------------- -- ------------- ------------- ----------
             Instrument                      VMC                 MASCS         STIS      
          Filter or channel                365-nm              VIS-VIRS        G430L     
                                          2006 May             2015 Dec.       2007         2011
                                         $-$2014 May          $-$2017 May     Jun. 5      Jan. 27
      Number of images/spectra              17742                1619           152          45
     Phase angle \[$^{\circ}$\]         52.3$-$138.6          0.9$-$144.0   89.7$-$90.0     79.3
       Latitude \[$^{\circ}$\]       Southern hemisphere                      0$-$30S     0$-$30S
   Longitude \[2[$^{\circ}$E ]{}\]         0$-$360              0$-$360      75$-$148     91$-$144
  --------------------------------- --------------------- -- ------------- ------------- ----------

  : Summary of 365-nm observations used in this study. \[tab:data\]

Methods {#sec:methods}
=======

For disk-resolved images, we Spectral data were convolved using the filter transmittance function of the 365-nm channel of VMC, and 2[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 [@Sromovsky01; @Garciamunoz14; @Garciamunoz17], without photometric correction . Radiative transfer model calculations were performed using the model and gaseous database in @Lee15b [@Lee16] .

Disk-resolved albedo {#subsec:disk-res_albedo}
--------------------

### Reflectivity {#subsubsec:reflectivity}

### Photometric corrections {#subsec:photo-corr}

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 ($\alpha$) angles, which can be described by a photometric law (disk function), $D(\mu, \mu_0, \alpha)$, where $\mu=\cos(e)$ and $\mu_0=\cos(i)$. We can separate albedo, $A(\alpha)$, and disk function, $D(\mu, \mu_0, \alpha)$, from the radiance factor ($r_{\rm F}$) that is derived from the measured radiance [@Shkuratov11], $$\centering
r_{\rm F}=A(\alpha)D(\mu, \mu_0, \alpha).
\label{eq:RF_AD}$$ This albedo $A$, the equigonal albedo, depends on $\alpha$ [@Shkuratov11].

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 [@Lee15a; @Lee17]. Therefore, in this study we adopted the LLS [@Buratti83; @McEwen86], $$\centering
D_{\rm LLS}=k(\alpha)\frac{2\mu_0}{\mu_0+\mu}+(1 - k(\alpha))\mu_0,
\label{eq:D_LLS}$$ where $k$ is a coefficient depending on $\alpha$, $$k(\alpha)=0.216004+0.00194196\times\alpha-(2.11589\times10^{-5})\times\alpha^{2},
\label{eq:k_LLS}$$ and $\alpha$ is in \[degree\]. Eq. \[eq:k\_LLS\] is updated from @Lee17, using more images to find the mean condition along phase angle.

### Areas of disk-resolved albedo for a comparison {#subsubsec:disk-res}

Figure \[fig:img\_comp\] shows albedo $A$ at three 5-degree phase angle bins: 75–80, 80–85, and 85–90$^{\circ}$, 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 [@Ogohara17], VMC data do not. So we restricted VMC data satisfying $i<84^{\circ}$, $e<81^{\circ}$, and $r_{\rm F}>0.05$ [@Lee15a]. The last condition causes non-smooth terminator for VMC images in figure \[fig:img\_comp\]. 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 [@Bertaux16], 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 images, two broad latitude bins: low (0-30$^{\circ}$) and high (50-70$^{\circ}$) latitudes. We derive a low latitudinal albedo also from the spectral data to complete the cross comparison (Section \[subsubsec:CrossComp\_Disk-resA\]).

From the equatorial orbit, UVI images show that cloud top albedo and contrast patterns primarily displayed north-south symmetry. However, as the example in figure \[fig:img\_hasym\] 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 [@Horinouchi18]. 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 {#subsubsec:CrossComp_Disk-resA}

The observed long-term decrease of albedo had been previously attributed to the sensitivity degradation 2[of VMC ]{}by @Shalygina15. 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\times10^{-5}$ orbit$^{-1}$, and 2[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 @Shalygina15. $$\centering
B(t)=B_0(t)\beta,
\label{B_t}$$ where $t$ is time (orbit), $B$ and $B_0$ are the corrected and observed radiance respectively (W m$^{-2}$ sr$^{-1}$ $\mu$m$^{-1}$), and $\beta$ is the 2.32 calibration correction factor. Temporal sensitivity degradation correction is given as $$\centering
B(t_2)=B(t_1)\left(\frac{1+k_d t_2}{1+k_d t_1}\right),
\label{eq:B_corr}$$ where $k_d$ is the sensor degradation factor (orbit$^{-1}$). We can get $B(t=0)$ as $$\centering
B(t=0)=B(t)\left(\frac{1}{1+k_d t}\right).
\label{eq:B_0}$$ So the final form is $$\centering
B(t=0)=B_0(t)\beta\left(\frac{1}{1+k_d t}\right).
\label{eq:B_0_end}$$

Figure \[fig:pha78N88\_discarded\] shows a comparison of low and high latitudinal albedo at the same phase angle bins, using the corrected VMC data with Eq. \[eq:B\_0\_end\] [@Shalygina15]. 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$\pm$0.822, for the initial calibration correction factor ($\beta$) of the 365 nm channel of VMC [@Titov12; @Shalygina15]. 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 \[tab:data\]). To limit uncertainties that may arise from the influence of the aerosol scattering along phase angle ($\alpha$), we utilize only VMC data obtained at the same phase angle as either the MASCS ($\alpha=85-90^{\circ}$) or STIS ($\alpha=75-90^{\circ}$) 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 $\alpha=85-90^{\circ}$, and $A$(STIS)/$A$(VMC) at $\alpha=75-80^{\circ}$, 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 \[fig:globe\_spec\]), 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 \[fig:pha78N88\]), 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 \[fig:pha78N88\]). Thus, we adopt this value and apply a new modified VMC calibration correction factor of 1.58 ($\beta=2.0\times0.79$) to the VMC data used in our study. As figure \[fig:pha78N88\] shows, when the modified calibration correction factor is applied, the 365 nm albedo observed by VMC and UVI are reasonably aligned; those in 2008-2009 (VMC) are overlapped data in 2016 (UVI).

Whole-disk albedo {#subsec:whole-diskA}
-----------------

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 [@Nakamura14]. 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_{\rm whole-disk}$, which is a function of phase angle ($\alpha$), following the equation below [@Sromovsky01], $$\centering
A_{\rm whole-disk}(\alpha)=\frac{\pi}{\Omega_{\rm Venus}}\frac{{d_{\rm Venus}}^2F_{\rm Venus}}{S_{\odot}},
\label{eq:APhi}$$ where $d_{\rm Venus}$ is the Venus distance to the Sun (AU), $F_{\rm Venus}$ is the measured disk-integrated Venus flux (W m$^{-2} \mu$m$^{-1}$), $\Omega_{\rm Venus}$ is the solid angle of Venus as seen from spacecraft (sr), and is the solar irradiance at 1 AU (W m$^{-2} \mu$m$^{-1}$)2[, which is calculated for either UVI or VMC, using each of transmittance functions. ]{}.

The observed solid angle of Venus, $\Omega_{\rm Venus}$, is calculated as $$\centering
\Omega_{\rm Venus}=\pi \left(\sin^{-1}\left(\frac{r_{\rm Venus}}{d_{\rm V-sc}}\right)\right)^2,
\label{eq:Omega_venus}$$ where $r_{\rm Venus}$ is the cloud top level radius of Venus (6052+70 km), and $d_{\rm V-sc}$ is the distance of spacecraft from Venus (km). Observed Venus flux, $F_{\rm Venus}$, is calculated as $$\centering
F_{\rm Venus}=\sum_{r<r_o}R_{\rm obs}(x,y)\times\Omega_{\rm pix},
\label{eq:F_venus}$$ where $(x,y)$ is a location of , $\Omega_{\rm pix}$ is a solid angle of one pixel of either VMC or UVI, is radiance (W m$^{-2}$ sr$^{-1}$ $\mu$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$^{\circ}$). $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\Phi(\alpha)$, where $A_g$ is geometric albedo, and $\Phi(\alpha)$ is the phase law of Venus [@Garciamunoz14; @Garciamunoz17].

### Comparison of whole-disk albedo of VMC and UVI {#subsubsec:comp_whole-diskA}

We calculated the whole-disk albedo, $A_{\rm whole-disk}$ (Eq. \[eq:APhi\]), 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 the field of view of cameras (see figure \[fig:img\_comp\]). The 1.58 modified calibration correction factor is applied to the VMC data.

Figure \[fig:diskint\_phase\] shows the results of the whole-disk albedo 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 \[sec:data\]). Fractions of disk illuminated by the Sun changes the solar phase angle, from 100% at 0$^{\circ}$ phase angle to 0% at 180$^{\circ}$ 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 [@Garciamunoz14; @Lee17]. The polynomial fit of UVI’s 2015–2017 data shows the empirical phase function $\Phi(\alpha)$ in the 40–110$^{\circ}$ phase angle range (brown solid line, Eq. is shown in Table \[tab:poly\_ci\]). We shift this phase function vertically to fit the maximum and minimum VMC whole-disk albedo at 75–80$^{\circ}$ phase angle (cyan lines), assuming a change of geometric albedo $A_g$ over time whereas the scattering properties $\Phi$ 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$^{\circ}$ phase angle, which is a geometric albedo $A_g$, increases from $\sim$0.33 in 2011 to $\sim$0.40 in 2015–2017. This result confirms the increasing albedo trends 2[from 2011 to 2015–2017, ]{}observed at high and low latitudes based on the 2[recalibrated disk-resolved 2011 VMC and disk-resolved 2015–2017 UVI data]{} (figure \[fig:pha78N88\]). 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 [@Shalygina15]. 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 of the albedo at small phase angles150–3152[$^{\circ}$E]{} range in 2011 March, 180–3302[$^{\circ}$E]{} in 2016 May, 165–(360)–452[$^{\circ}$E]{} 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 {#subsec:RTM}
------------------------

We use a one-dimensional line-by-line radiative transfer model (SHDOM, @Evans98 [-@Evans98]) in the 0–100 km altitude range and in the 2000–50000 cm$^{-1}$ ($=0.2-5~\mu$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 @Lee15b [@Lee16]. CO$_2$ line parameters were taken from a combined HITEMP2010 [@Rothman10], and one developed by @Wattson92 and @Pollack93, as described in @Lee16. We included collision-induced CO$_2$ absorption in near-infrared, and that of H$_2$O [@Lee16]. Line parameters of other gases, N$_2$, SO$_2$, OCS, HCl, CO, HF, and H$_2$S, were taken from HITRAN2012 [@Rothman13], and vertical profiles of gaseous abundances were taken from @Titov07. We included Rayleigh scattering [@Pollack67; @Hansen74], and UV range absorption cross sections of SO$_2$ [@Wu2000]. Microphysical properties of the cloud aerosols (mode 1, 2, 2’, and 3) were taken from . We took the vertical structures of clouds extinction coefficient from @Crisp86. For the unknown absorber, we assumed @Crisp86’s absorption coefficient ($Q_{\rm abs}$) of the mode 1 particle (0.15-$\mu$m mean radius cloud particles) in the 0.3–0.8 $\mu$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 [@Crisp86; @Lee15b; @Haus16].

Results {#sec:results}
=======

Temporal variations of low and high latitudinal albedo {#subsec:disk-resolved-results}
------------------------------------------------------

The albedo $A$ (Eq. \[eq:APhi\]) is phase angle $\alpha$ dependent due to strong backward and forward aerosol scatterings [@Lee15a; @Shalygina15; @Garciamunoz14]. Therefore, we restrict a comparison of $A(\alpha)$ to data obtained at near equivalent phase angle bins. Figure \[fig:pha78N88\] shows the temporal evolution of Venus’ low and high latitudinal mean albedo obtained at phase angles of 75–80$^{\circ}$ and 85–90$^{\circ}$. In this figure a strong and steady decline in the albedo occurs between 2006 and 2011, from 0.4 to 0.25 at low latitudes and from 0.6 to 0.3 at high latitudes. Albedos at low and high latitudes in 2015 December–2017 May are restored to the 2008-2009 values.

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, $\bar{A}_{\rm UVI}(\alpha)$, which is derived as a polynomial fit of UVI data in 2016–2017 in the 50–100$^{\circ}$ phase angle range (figure \[fig:EQ\_CC\_phasefn\]). Table \[tab:poly\_ci\] 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: $$\centering
{\rm Deviation~[\%]}=\frac{A(\alpha)-\bar{A}_{\rm UVI}(\alpha)}{\bar{A}_{\rm UVI}(\alpha)}\times100.
\label{eq:dev_A}$$

           $\bar{A}_{\rm UVI, low-lat}(\alpha)$   $\bar{A}_{\rm UVI, high-lat}(\alpha)$   $\bar{A}_{\rm UVI, whole-disk}(\alpha)$
  ------- -------------------------------------- --------------------------------------- -----------------------------------------
   $c_0$                0.18708465                                                                      0.62539023
   $c_1$               0.0086404233                                                                   $-$0.022833883
   $c_2$             $-$0.00017491717                                                                  0.00063482472
   $c_3$         1.1079428$\times10^{-6}$                                                       $-$1.2777032$\times10^{-5}$
   $c_4$                                                                                         1.7206082$\times10^{-7}$
   $c_5$                                                                                        $-$1.4916523$\times10^{-9}$
   $c_6$                                                                                         7.9181817$\times10^{-12}$
   $c_7$                                                                                       $-$2.3230007$\times10^{-14}$
   $c_8$                                                                                         2.8698482$\times10^{-17}$

  : Coefficients of a polynomial fit, 2[$y(\alpha)=\sum_{i=0}^{n}{c_i}{\alpha}^i$, where $n\leq8$,]{} in Figures \[fig:diskint\_phase\] and \[fig:EQ\_CC\_phasefn\]. \[tab:poly\_ci\]

Figure \[fig:EQ\_CC\_devA\] shows percent of deviations of mean $A$ from $\bar{A}_{\rm UVI}(\alpha)$ at low and high latitudes as a function of time, where phase angles from 50 to 110$^{\circ}$ are represented according to the color bar at the top of each panel. The overall decline in the 365-nm albedo from 2006 to 2011 remains apparent. Additionally, relatively sharp albedo declines are observed at the end of 2009 and 2010, and at the beginning of 2013 that remain constant over short time periods ($\sim$months) before returning to the $\bar{A}_{\rm UVI}(\alpha)$ level. The robustness of the darker albedo conditions observed at the 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 [@Imamura01], and leaving behind only the bright newly formed aerosols that support the existence of Venus’ bright high latitudes.

Temporal variations of whole-disk albedo {#subsec:whole-disk-results}
----------------------------------------

We apply Eq. \[eq:dev\_A\] also for the whole-disk albedo $A_{\rm whole-disk}$ and use the mean phase curve (figure \[fig:diskint\_phase\] and Table \[tab:poly\_ci\]) to get the percent deviations of $A_{\rm whole-disk}$ from the 2016–2017 mean phase curve. Figure \[fig:wholeDiskA\_dev\] 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 20062[–]{}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 {#subsec:RTMresults}
------------------------

Because of the influence of the unknown absorber on solar heating rate near the cloud top level [@Crisp86; @Titov07; @Lee15b] owing to the broad absorption spectrum from UV to visible [@Crisp86; @Perezhoyos18], 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 \[subsec:RTM\]) 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 \[subsec:RTM\]), 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 the assumed initial absorption coefficient $Q_{\rm abs}$ taken from @Crisp86 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 $\alpha=88^{\circ}$: $(e, i)=(28^{\circ}, 60^{\circ}), (38^{\circ}, 50^{\circ}), (48^{\circ}, 40^{\circ}), (58^{\circ}, 30^{\circ}),$ and $(68^{\circ}, 20^{\circ})$, 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 \[subsec:photo-corr\]), and the mean value is used to find $f$ values that match the observed albedo.

Figure \[fig:f\] shows the relationship of $f$ and the mean value of calculated 365-nm albedo for the 88$^{\circ}$ phase angle conditions (Figure \[fig:pha78N88\]). 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$^{\circ}$ latitude, which is the middle of the low latitude bin, at local noon time. Figure \[fig:divfnet\] shows the net solar flux divergence spectrum as a function of altitude $z$, $-\nabla \cdot F_{\rm net}=-({dF_{\rm net}}/{dz})$, where $F_{\rm net}$ is net solar flux. The strongest influence of the unknown absorber appears around 400 nm, where decreasing absorption and increasing solar irradiance .

![Net solar flux divergence as functions of wavelength and altitude at local noon time at 15$^{\circ}$S. (a) f=0.65, (b) f=1.18, and (c) f=2.51 are multiplied the extinction coefficients of the unknown absorber ($Q_{\rm abs}$) in the 57–71 km altitude range and in the 0.3–0.8 $\mu$m range (see text for details). Spectral features are smoothed over 0.01 $\mu$m interval. \[fig:divfnet\]](divFnet.pdf){height="0.7\paperheight"}

Figure \[fig:HR\] shows the calculated local noon time solar heating rate at 15$^{\circ}$ latitude, derived from Figure \[fig:divfnet\]. This represents that the solar heating rate varied from 2[$-$]{}25% to +40% from the mean (2006–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 [@Lee15b], but the structure of the upper clouds at low latitudes from near infrared observations is shown to be rather stable during the time of Venus Express [@Ignatiev09; @Cottini15; @Fedorova16], 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 [@Molaverdikhani12; @Lee15a]. 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.

Discussion {#sec:discussion}
==========

Relationship between solar heating and zonal winds at the cloud top level {#sec:dis1_solarheating-winds}
-------------------------------------------------------------------------

The observed 365-nm albedo variations should in solar heating variations near the cloud top level as shown in Section \[subsec:RTMresults\]. We note that two long-term trends occurred in parallel; in the period when the 365 nm albedo declined, leading to increases in solar heating, the long-term cloud top zonal wind was observed to increase from 80–90 m/s to $\sim$110 m/s between 2007 and 2012 around local noon at 20$^{\circ}$S [@Khatuntsev13; @Kouyama13; @Hueso15]. Likewise, between 2014 and 2016, when the 365 nm albedo increased, leading to decreases in solar heating, the zonal wind speed was observed to slow down from 110 to 100 m/s [@Horinouchi18]. 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 [@GarateLopez18]. Starting from the reference simulation (see Appendix \[Appendix1\] 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$^{\circ}$S and 20$^{\circ}$S and at 30 mbar ($\sim$70 km, near the cloud top level) is plotted in figure \[fig:zonave\_T\_U\_dTsw\], 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 \[Appendix1\].

Possible causes of the observed albedo variations {#sec:dis2_reason}
-------------------------------------------------

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 [@Marcq13], and then subsequently increase around 2016 [@Encrenaz19]. 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 [@Jiang18] as shown in figure \[fig:UV\_Solar\_Cosmic\]. This resembles the correlation of Neptune’s reflectivity and the solar activity cycle [@Aplin16]. Since solar EUV radiation might affect photochemical reactions involving SO$_2$ that are necessary for aerosol formation on Venus [@Mills07; @Parkinson15], 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 [@Lovejoy04; @Kirkby07], and there are observations on H$_2$SO$_4$-H$_2$O ultrafine aerosols of less than 9 nm in diameter in Earth’s upper troposphere and stratosphere that were explained by the ion-induced nucleation [@Lee03]. The peak of ion production rate in the Venusian atmosphere due to galactic cosmic rays is predicted 62.5 km [@Nordheim15] 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 \[fig:UV\_Solar\_Cosmic\] shows comparisons of the 365-nm albedo at low latitudes, neutron cosmic ray detected from the Oulu station[^1], and Lyman-$\alpha$ 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-$\alpha$ 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-$\alpha$ 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 {#sec:dis3_comparisonwithotherplanets}
-----------------------------

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 [@Stephens14]. 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 $\pm$0.02 (corresponding to $\pm$2% changes in flux) at the blue (472$\pm$10 nm) and red (551$\pm$10 nm) filters over 1972–2014 [@Aplin16], and at Mars where the surface albedo varied by 10% at the red filters (575–-675 and 550–700 nm) from 1976–1980 to 1999–2003 [@Geissler05].

Further studies {#sec:dis4_furtherstudies}
---------------

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 may explain unbalanced net radiative energy below the clouds [@Lee17], 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 [@Haus17], 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 2[impact]{} Venus’ climate.

Summary {#sec:summary}
=======

We present the intense decadal variation of Venus’ 365-nm albedo between 2006 and 2017; the maximum albedo occurred in 2006–2007, the minimum in 2011–2014, and the recovery of albedo in 2016–2017 to the level in 2008–2009. 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 [@Shalygina15], and propose a new calibration correction factor for VMC in this study. The ranges of albedo variations are $\sim$0.2–0.4 at low latitudes, and $\sim$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 2[$-$]{}30% to +20% compared to the mean value in 2016–2017, meaning that the albedo variation occurred 2[on]{} a global scale.

Our one dimensional line-by-line radiative transfer model calculations 2[reveal]{} this level of albedo variation can alter solar heating rate from 2[$-$]{}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 2[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 $\sim$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.

Brief description on IPSL Venus GCM, and analysis of the simulation results. {#Appendix1}
============================================================================

Institut Pierre Simon Laplace (IPSL) Venus Global Climate Model (GCM) has successfully demonstrated the development of strong zonal wind near the cloud top level [@Lebonnois16]. Upgrading the GCM, [@GarateLopez18] employed latitudinally varying cloud structures [@Haus14]. The authors also used solar heating look-up table [@Haus15b] according to the latitudinal cloud structures. Thermal cooling is based on the net-exchange rate (NER) formalism [@Eymet09] 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 [@GarateLopez18].

In this latest IPSL Venus GCM, only solar heating has been reduced, the expected solar heating rate variations shown in this manuscript. Figure \[fig:zonave\_T\_U\_dTsw\] 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 \[fig:GCM\_result4\] at 30 mbar.

Looking at the different heating rates (figure \[fig:GCM\_result4\]a), 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 \[fig:GCM\_result4\]b). 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 \[fig:GCM\_result4\]c, except between 30 and 50$^\circ$ 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 \[fig:GCM\_result4\]d). Regardless 2[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 [@Migliorini12] or radio occultation temperature profiles from VeRa/Venus Express [@Tellmann12].

Y.J.L thanks T. Satoh, R. Lorenz, N. Ignatiev, R. Hueso, and C.C.C. Tsang for helpful comments. K.-L. J. was supported by NASA/PS program, grant NNX16AK82G, for the analysis of the 365 nm STIS data. S.P.-H. was supported by the Spanish MICIIN project AYA2015-65041-P (MINECOFEDER, 695UE), and Grupos Gobierno Vasco IT-765-13. S.Le. acknowledges support from INSU-PNP, IPSL VGCM computations were done thanks to the High-Performance Computing (HPC) resources of Centre Informatique National de l’Enseignement Supérieur (CINES) under the allocation n$^{\circ}$A0040110391 made by Grand Équipement National de Calcul Intensif (GENCI). J.P. the JAXA’s International Top Young Fellowship. Solar heating calculations were performed on a Supermicro SuperServer Intel(R) Xeon(R) CPU E5-2620 v4 funded through JAXA’s International Top Young Fellowship. T. H., M. T., and T. I. acknowledge JSPS grant 16H02231. S.Li. was supported by a NASA Participating Scientist in Residence Grant (NNX16AC79G). E.M acknowledges support from CNES and INSU-PNP for SPICAV data analysis.

[^1]: <http://cosmicrays.oulu.fi>

[^2]: <http://lasp.colorado.edu/lisird/data/composite_lyman_alpha/>
