Issue 
A&A
Volume 626, June 2019



Article Number  A30  
Number of page(s)  9  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201935388  
Published online  06 June 2019 
Principal components of shortterm variability in the ultraviolet albedo of Venus
^{1}
Department of Complexity Science and Engineering, Graduate School of Frontier Sciences, The University of Tokyo,
Kashiwa,
Chiba Prefecture
2778561,
Japan
email: pushkarkopparla@gmail.com
^{2}
Institute of Space and Astronautical Science (ISAS), Japan Aerospace Exploration Agency (JAXA),
Sagamihara,
Kanagawa Prefecture
2525210,
Japan
Received:
2
March
2019
Accepted:
15
April
2019
We explore the dominant modes of variability in the observed albedo at the cloud tops of Venus using Akatsuki UVI 283 nm and 365 nm observations, which are sensitive to SO_{2} and unknown UV absorber distributions respectively. The observations, taken over the period Dec. 2016 to May 2018, consist of images of the dayside of Venus, most often observed at intervals of two hours, but interspersed with longer gaps. The orbit of the spacecraft does not allow for continuous observation of the full dayside and the unobserved regions cause significant gaps in the datasets. Each dataset is subdivided into three subsets for three observing periods, the unobserved data are interpolated, and each subset is then subjected to a principal component analysis to find six oscillating patterns in the albedo. Principal components in all three periods show similar morphologies at 283 nm, but are much more variable at 365 nm. Some spatial patterns and the timescales of these modes correspond to wellknown physical processes in the atmosphere of Venus such as the ~4day Kelvin wave, 5day Rossby waves, and the overturning circulation, while others defy a simple explanation. We also a find a hemispheric mode that is not well understood and discuss its implications.
Key words: planets and satellites: individual: Venus / planets and satellites: atmospheres / techniques: photometric / methods: statistical
© ESO 2019
1 Introduction
The atmosphere of Venus, as seen in ultraviolet wavelengths, is long known to have striking albedo patterns that are indicative of the dynamics and chemistry of the upper atmosphere (Ross 1928; Del Genio & Rossow 1990; Markiewicz et al. 2007). Several of these features have been identified and named, for instance: the Yfeature, associated with equatorial Kelvin wave (Del Genio & Rossow 1990; Peralta et al. 2015); midlatitude Rossby waves (Del Genio & Rossow 1990); and the equatorial and polar caps and bands (Rossow et al. 1980). The main absorbers in the near UV, SO_{2} and the unknown UV absorber, seem to vary on multiple timescales, from daily to multiyear (Encrenaz et al. 2012; Marcq et al. 2013; Lee et al. 2015) and their generation, destruction, and advection in the atmospheric flow are responsible for the albedo patterns.
Approximately 4–5 Earthday period variations are known to be related to the atmospheric background flow and short period Kelvin and Rossby waves (Khatuntsev et al. 2013; Kouyama et al. 2015; Peralta et al. 2015), but not fully understood since the relationship between the UV albedo contrast and physical quantities describing the wave field (e.g., temperature, pressure, and velocity) is not known. The causes for the longer period variability are uncertain, such as 255 days variability in zonal wind speed (Kouyama et al. 2013) and ~270 days in the 365 nm latitudinal contrast (Lee et al. 2015) and mesospheric SO_{2} gaseous abundance (Marcq et al. 2013). With sustained observations over a threeyear period now available from the Akatsuki Ultraviolet Imager (UVI) instrument (Nakamura et al. 2016; Yamazaki et al. 2018), we have a long enough baseline of observations to find the leading modes of oscillation in the albedo along with their periodicities and associated spatial structures. We use observations at 283 nm, which are correlated to SO_{2} abundance above the clouds (with some absorption by unknown UV absorber and ozone), and those at 365 nm, which is close to the peak absorption by the unknown UV absorber. With these data, we attempt to answer two main questions: first, whether we can resolve the variable UV patterns into a small number of recognizable components; and second, what these components can tell us about the physical processes active in the atmosphere of Venus.
The next section details the data reduction, including the handling of missing data and principal component analysis (PCA). The following section deals with the leading oscillations that result from the PCA, their physical interpretations, and the statistical significance of the results against noise. The final section summarizes the findings and speculates on future applications of such analyses for Venus climate studies.
2 Data and methods
2.1 Datareduction and normalization
The UVI on Akatsuki has two filters at 283 and 365 nm, which correspond to the absorption bands of SO_{2} and the unknown UV absorber. The field of view (FOV) is 12° × 12° and can observe the whole Venus disk except for about eight hours near the periapsis. The imaging area is composed of 1024 × 1024 pixels and has an angular resolution of 12 × 10^{−2} degrees per pixel, which corresponds to spatial resolutions of about 200 m and 86 km on the cloud top level in the observations from the altitudes of ~1000 km at the periapsis and 60 Venus radii at the apoapsis (Yamazaki et al. 2018). We use the Level 3 (L3) 283 nm data from the UVI, internal release version 20180901. The L3 data from Akatsuki have observed variables mapped onto a regular (equispaced) longitudelatitude grid. The resolution of the L3 longitudelatitude map is 0.125° × 0.125° (2880 × 1440) grids for 360° longitude and 180° latitude (Ogohara et al. 2017). The orbit of the spacecraft is elliptical with a revolution period of ~10.5 Earth days,a periapsis altitude of 1000–8000 km, and an apoapsis altitude of ~360 000 km (Nakamura et al. 2016). The data used in this study were taken on the apoapsis side, where the FOV of UVI can capture the whole Venusdisk. Since the major axis of the orbit is roughly fixed in the inertial coordinate, the observation condition changes with the revolution of Venus around the Sun. A visualization of the dataset is provided in Fig. 1.
To identify spatial structures and time series of oscillations, it is ideal to have continuous observations of the full dayside at high frequency over a long time period. However, owing to practical constraints, the fraction of dayside observed varies systematically and there are often long gaps in the time series. The majority of observations have a gap of 2 h between them, though the longest gap between successive observations is about 40 days, which is the time of missing global dayside images due to the highly elliptical equatorial orbit. Images that cover less than 80% of the dayside observed are excluded from this analysis and grid points on a given image where the dayside is not observed are referred to as missing data points. The dayside is the region between solar local time 600–1800 h, and the fraction of the dayside observed is calculated as the ratio of grid points in this region with observed values to the total number of grid points on the dayside. The exclusion threshold could be relaxed to 10− 60%, but doing so increased the ratio of missing data to observed data and introduced too many artificial patterns during the interpolation process described in Sect. 2.2.2. On the other hand, as the exclusion threshold increases, the number of images that can be included in the analysis decrease. As a result, the statistics become more uncertain and the time baseline shorter, limiting the longest period oscillation that can be studied. The 80% threshold offers a reasonable balance between these two competing effects.
We convert the radiances to relative albedos based on the Minnaert law described in Lee et al. (2017), using Eqs. (8) and (10) of that paper for 283 and 365 nm, respectively. The albedos in each observation thus derived are normalized by dividing by the spatial mean value in each map to give a relative albedo. We chose to normalize by the spatial mean instead of the commonly used spatial maximum, since the maximum is very unstable in fields with a large number of missing values. Alternately, another stable measure of the top values of the albedo distribution (such as the 90th percentile) could be used for normalization. We note that this normalization allows for values of albedo above 1 and removes trends in albedo with time, but that is not our primary concern at this time. We are interested in planetary scale patterns in albedo, for which the relative albedos are sufficient. The maps are then centered such that the subsolar longitude is always at 180° or 1200 local time. This is equivalent to plotting on a local solar time grid. Grid points where the cosines of the solar zenith angle or the viewer zenith angle are ≤0.2 are excluded, since the observing error at such points is rather large compared to natural variability (Fukuhara et al. 2017). The data are regridded to a 2° resolution (using the resize function from the Python library skimage.transform). The region of interest for our study is the dayside, which we define to be between 800–1600 h local time. The 283 nm dataset contains 665 images, each with 4476 grid points on the dayside. Of these 2985492 points, 8% are unobserved (“missing data”). The 365 nm dataset contains 652 images and has a similar fraction of missing data.
Fig. 1 Observation times of the 665 and 652 images used in this study for 283 nm (top left) and 365 nm (bottom left), respectively.The windows consist of 209, 244, and 212 images for 283 nm and 206, 237, and 209 images for 365 nm, respectively.The color of the data point indicates the fraction of the dayside that was observed in that image. Right: histograms showing the time gaps between successive images. The majority of the images are taken at intervals of 2 h, but longer intervals are not uncommon. 
2.2 Principal component analysis
A powerful dimensionality reduction technique, PCA is widely used to find oscillations in atmospheric and climate data on Earth (Wilks 2006). Our PCA follows this procedure:
 1.
The twodimensional (latitude longitude grid) images containing relative albedos are flattened into onedimensional arrays.
 2.
A data matrix is created with p variables(here the number of spatial gridpoints, 4476) as rows and n observations (here the number of images) as columns. An element x_{i,j} of this matrix corresponds to the ith gridpoint and jth observation.
 3.
The time mean (rowwise mean) of each variable is subtracted to generate the anomalies from the mean as follows: (1)
 4.
The covariance matrix is calculated with the dimension p × p. An element of this matrix c_{k,l} is the covariance between the observations at gridpoints k and l, which is calculated as (2)
 5.
The eigenvalues and eigenvectors of the covariance matrix are calculated (we use the eigsh function from the Python library scipy.sparse.linalg). The eigenvectors are called the principal components (PCs) or empirical orthogonal functions (EOFs), which give the spatial structure of the leading modes of variability. Each PC has a dimension of p × 1 and there are p PCs in total.
 6.
The dot product of the meanremoved data with the PCs give the PC loadings (hereafter referred to as loadings). The loadings give the time series of contribution of each mode of variability to the data, each of which has a shape 1 × n, and there are p loadings in total. The loadings are derived from the EOFs, ɛ, and dataanomalies, x^{a} as follows: (3)
 7.
The original data can be fully reconstructed from the PCs and loadings thus, (4)
The PCs corresponding to the largest few eigenvalues and their loadings give the spatial patterns and the time series of the leading oscillations of the dataset, respectively. The fraction of total data variance explained by each eigenvector is given by the ratio of its eigenvalue to the sum of all eigenvalues.
2.2.1 Calculation of the covariance matrix
The covariance of albedo values between two pixels (true covariance or population covariance) has to be approximated using the limited data sample available (sample covariance). However, the existence of missing data means that the sample covariance matrix cannot be calculated immediately from the data. The simplest approaches for dealing with missing data in the calculation of the sample covariance matrix are deletions; that is, removing variables with missing data points (known as listwise deletion) or the calculation of covariances by considering only data points for which data exist for both variables of interest (pairwise deletion) (e.g., Carter 2006; Nakagawa & Freckleton 2008). Listwise deletion is undesirable, since it discards a grid point entirely even if that point has missing data at only a single time step. Also, the use of pairwise deletion introduces undesirable properties to the sample covariance matrix, such as the existence of negative eigenvalues. This indicates that the sample covariance is not a good estimate of the true covariance of the system being observed and the sample covariance matrix lacks the properties of a true covariance matrix such as positive semidefiniteness (only positive or zero eigenvalues) (Pourahmadi 2011). Using such bad estimates often results in the first few eigenvalues being biased high.Meanwhile, the existence of negative eigenvalues makes the calculation of the variance associated with each eigenvalue hard to interpret, since negative variances do not have a clear physical meaning (Dray & Josse 2015).
A better way to deal with missing values is imputation, either by interpolation from available data or the imposition of additional constraints on the structure of the covariance matrix (regularization). Several different techniques have been proposed to regularize the properties of such illconditioned sample covariance matrices such as ridge regressions/Tikhonov regularization, banding, and tapering (Bickel et al. 2008; Warton 2008). These methods are primarily focused on eliminating bad behavior in the covariance matrix by introducing constraints such as smoothness or sparsity, where the constraints may or may not be determined from the dataset. Other approaches, often used in climate studies, fill in missing data through various interpolations, such as nearest neighbor regressions followed by smoothing [DCTPLS] (Garcia 2010; Wang et al. 2012); optimal interpolation (Burgess & Webster 1980); singular spectral analysis (Kondrashov & Ghil 2006); or iterative techniques like regularized estimation maximization [RegEM] (Schneider 2001), data interpolating empirical orthogonal functions [DINEOF] (Beckers & Rixen 2003), and several other such variants of PCA based interpolations (Ilin & Raiko 2010). For our dataset, we need an interpolation technique that must be able to handle a large percentage of missing data and continuous data gaps. It must also be able to interpolate effectively in three dimensions (as opposed to purely spatial interpolation). The method must derive the regression parameters from the data itself, i.e., be nonparametric. Since the missing values lie on a smooth map, there is a roughness constraint on the interpolated values. This can be enforced in many ways: creating interpolations using loworder PC truncations (like DINEOF), ridge regularizing the covariance matrix (RegEM), or enforcing a roughness penalty by explicitly specifying a smoothness parameter (DCTPLS).
Since our ultimate goal is to perform a PCA on this dataset, two of the above methods, RegEM and DINEOF, iteratively create imputed datasets that converge on a set of PCs. But RegEM uses a 2D spatialonly linear regression that is unsuited for our data since missing values are primarily located near the edges of the region of interest rather than in the center; therefore we use the DINEOF method (Beckers & Rixen 2003). This method is closely related to other spatiotemporal methods like singular spectral analysis and optimal interpolation and produces comparable results (Allen & Smith 1996; AlveraAzcárate et al. 2005).
Fig. 2 Distribution of the 20 largest eigenvalues in the DINEOF interpolated datasets for 283 nm. The 365 nm datasets show similar trends. The red vertical line shows the point at which the eigenvalue curve has an inflection. The PCs corresponding to eigenvalues left of this line are considered significant by the Scree test. 
2.2.2 Data Interpolating Empirical Orthogonal Functions
The Data Interpolating Empirical Orthogonal Functions (DINEOF) method is based on the descriptions from Beckers & Rixen (2003) and AlveraAzcárate et al. (2005) and is implemented as follows:
 1.
All missing values in the timemean subtracted dataset are filled in with zeros and the covariance matrix is calculated.
 2.
The PCs, eigenvalues, and loadings of this zero filled covariance matrix are calculated as described in Sect. 2.1.
 3.
The missing values are imputed using a reconstruction consisting only of the first leading PC and corresponding loadings as follows: (5)
 4.
The imputed dataset is again used to calculate a new covariance matrix and the previous step is repeated to find better imputation values. This procedure is iterated until the imputed values converge. Convergence is defined to occur when the root mean square difference in imputed albedo values from two successive iterations differ by less than 10^{−4} or 10% over three successive iterations, computed as (6) where RMS_{(n,n−1)} is the root mean square difference in imputed values between iterations n and n − 1.
 5.
The next set of iterations then uses two PCs for the imputation (7) and is repeated to convergence.
 6.
Thus wecan arrive at a converged imputed dataset for any given number of eigenvectors.
Since it is possible to create an imputed dataset with any k PCs, where 1 ≤ k ≤ p, the optimal number of PCs is determined through a generalized crossvalidation procedure. Following AlveraAzcárate et al. (2005), we choose 1% of the existing data points randomly and set them as artificially missing. This number of points is often used for crossvalidation procedures (e.g., Wilks 2006). The imputed datasets are reconstructed using 1,2,3 … 100 PCs using the convergence criterion above. The artificially missing data is not included in the convergence calculation as described in Step 4 above. For each converged imputed dataset, the root mean square error of the reconstruction of the artificially missing data is calculated. The optimal imputation is that for which this RMS error is minimum. For our dataset, this minimum occurs at N ~ 25 PCs.
3 Results and discussion
3.1 Dominant modes and statistical significance
The question of how many PCs are significant and should be retained has been widely discussed, and many different stastical and ruleofthumb approaches exist (Jackson 1993; PeresNeto et al. 2005; Cangelosi & Goriely 2007). Of these we use a very simple criterion that looks for an inflection point in the eigenvalue spectrum to find the break between significant and noisy PCs. This is known as Catell’s Scree test (Cattell 1966) and is a commonly used heuristic significance metric. From Fig. 2, we see that our interpolated datasets have an inflection point around 5 or 6, indicating a dimension of 6. We therefore retain only the first 6 spatial PCs, and these are shown in Figs. 3 and 5.
The LombScargle periodograms of the time series corresponding to each PC are shown in Figs. 4 and 6, with a minimum period of 1 day and a maximum of half the total length of the window (typically each window is about 40 days). The periodograms show several peaks, but their significance must be first confirmed. The spectral analyses of atmospheric data are affected by red noise (e.g., Allen & Smith 1996; Meinke et al. 2005), where the noise contributionincreases with period as opposed to a flat spectrum for white noise. This leads to spurious long period peaks in periodograms. This noise is usually approximated by an autoregressive process with a time lag of one (AR1) (Gilman et al. 1963), however irregular time sampling prevents estimation of the autocorrelation coefficient directly from the data. We therefore use existing methods in the climate science literature (Mudelsee 2002; Schulz & Mudelsee 2002) for variable spaced time series data. We do not use the programs made available by these studies, but only use the methods as described in the text with some minor modifications. Our implementation is summarized below:
 1.
The AR1 process is modeled as (8)
where is the time series value at time t, is the time step, τ is persistence timescale of the time series, and ɛ is the Gaussian noise component.
 2.
The value τ is estimated by a leastsquares fit of the AR1 equation to the loadings (using the function curve_fit from the Python library scipy). There is one τ for each PC loading in each window.
 3.
One thousand artificial time series are generated using the AR1 process described above; the initial value and the time steps are the same as the PC loading. The Gaussian process is taken to have a variance of (9)
 4.
Each artificial time series is scaled so that the periodogram has the same area under the curve as the original time series. The mean, ⟨G_{am}(f)⟩, and the standard deviation, ⟨G_{as}(f)⟩, of the setof all the periodograms of the artifical time series are calculated, where f is the frequency.
 5.
The analytical expression for the power spectrum in the periodogram of an AR1 process is given by (10)
where f_{max} is the highest frequency in the periodogram (often taken to be the Nyquist frequency, but in this case set to 1 day^{−1}), G_{o} is the mean spectral power, ρ is the mean autocorrelation coefficient calculated as ρ = exp(−Δt_{mean}∕τ), where Δt_{mean} is the arithmetic mean of the time steps in the time series.
 6.
Correction factors are calculated for the deviations of the mean artificial periodogram from the analytical value (11)
These correction factors account for biases resulting from taking periodograms of irregularly spaced data.
 7.
The unbiased periodogram of the original series is calculated as G_{ts}(f)∕c(f), where G_{ts}(f) is the power spectrum of the original time series. The standard deviation is scaled similarly and plotted in Fig. 4. Assuming a Gaussian distribution, the 95% confidence level is taken to be two standard deviations from the mean. The periods of the significant peaks are labeled.
The results are shown in Figs. 4 and 6.
Fig. 3 First six modes of variability from the DINEOF interpolated UVI dataset for the three windows at 283 nm. The first row is lined up in order of the PCs, while the second and third are rearranged so that similar patterns line up in the same column. The first mode shows a pattern similar to a Hadley circulation, while the second show a hemispheric oscillation. The third, fourth, and fifth columns appear to show a combination of short period atmospheric waves, while the sixth shows the well known Yfeature. 
3.2 Physical interpretations of the oscillations at 283 nm
The first column in Fig. 3 shows a (roughly) longitudinally uniform latitudinal gradient from the equator to the poles. This likely reflects the Hadley circulation, which lifts SO_{2} from the lower atmosphere to the cloud tops at the equator from which it is advected to the poles while being lost to photodissociation and conversion to sulfuric acid (Marcq et al. 2013) and polysulfur (Mills et al. 2007). The corresponding periodicities in Fig. 4 show a major periodicity of around 10 days, which is likely due to the spacecraft orbital period (Nakamura et al. 2016). The synchronization with the orbital motion might indicate a solarlocked component, i.e., a contribution of thermal tides. The approximately 4day period can be attributed to the atmospheric rotation around the planet, and the 8day peak is probably a subharmonic of this period. However, GCM studies indicate the existence of barotropic or baroclinic waves with periods ranging from 6 to 23 days (Lebonnois et al. 2016) in the cloud layer, which is another possibility for the 8day signature. An ~ 8day period wavewas observed a little deeper than the cloud top level in nearinfrared observations (~1.7 μm) and is also thought to result from the interaction of Kelvin and Rossby waves with the mean flow (Hosouchi et al. 2012).
The second column shows a pattern broadly described as a hemispheric oscillation, showing an oscillation with opposite directions ofchange between the two polar regions. A polar asymmetric brightening or darkening was observed by a groundbased telescope (Dollfus 1975), and it was thought to be caused by polar clouds independently evolving in either hemisphere. The pattern identified in this work is somewhat different in that the two hemispheres do not evolve independently; rather they evolve together in opposite directions, that is, one brightens while the other darkens. Notably, Fig. 4 shows no significant periods associated with this pattern other than the spacecraft orbital period. This can mean either that the phenomenon is aperiodic, or more likely, that the period is longer than the observational window considered in this work (approximately 40 days). Such a large scale poletopole pattern could indicate the existence of an atmospheric teleconnection. Another interpretation might be that the period comparable to the spacecraft orbit might indicate a contribution of asymmetric components of thermal tides. Alternatively, this pattern may also arise from the nonequatorially symmetric component of the meridional circulation. However, further observationsare required for confirmation, as completely symmetric polar features were observed from the midinfrared groundbased observations (Sato et al. 2014). A very recent study of IR1 (900nm) images from Akatsuki revealed an asymmetric pattern in the middle clouds at low latitudes with a periodicity of 4–5 days (Peralta et al. 2018). It is unclear at this time if this IR pattern is directly related the hemispheric mode we find at the cloud tops. Analysis of midinfrared images (Longwave InfraRed Camera of Akatsuki; Fukuhara et al. 2011) from the same time period as the UV images studied in this work, using PCA, may confirm this oscillation and offer clues to robust physical interpretations in near future. Notable hemispheric dichotomies have also been observed in CO concentrations below the clouds (Arney et al. 2014; Marcq et al. 2006, 2008). Higher concentrations were found to variably occur either on the northern and southern hemispheres during observations separated by several months. SO_{2} hemispheric dichotomies were also studied, but an unambiguous detection was only made once in 2010. It is unclear at this time how these longterm variations are related to the shortterm mode we find.
We interpret the third, fourth, and fifth columns to be combinations of atmospheric wave patterns. They often show shortperiodvariabilities of approximately 5 days. Waves with periods of a few days have been observed (Del Genio & Rossow 1990; Kouyama et al. 2015; Imai et al. 2016) and simulated in Venus GCM models (Yamamoto & Takahashi 2006). In particular, the 5day wave has been interpreted as a Rossby wave and the 4day wave as an equatorial Kelvin wave (Rossow et al. 1990; Kouyama et al. 2012; Imamura 2006). It must also be noted that several subplots do not show any notable periodicity even similar spatial patterns are associated with waves in other windows. For example, in column three, row one shows no periodicity except for the spacecraft orbit and some high frequency noise around the 1–2 day range. But rows two and three clearly show 5 and 4day periods, respectively.It is interesting that similar spatial patterns appear to form from different wave configurations in different periods.
The sixth column clearly shows the famous Yfeature (Boyer & Camichel 1961; Rossow et al. 1980), and is clearly associated with a strong periodicity of about four days. This is also consistent with the interpretation of the major cloud features of the Venus atmosphere interpreted as a trapped Kelvin wave (Del Genio & Rossow 1990; Kouyama et al. 2012; Peralta et al. 2015). It is interesting that the order of PCs in each window is slightly variable and needs rearrangement to align similar patterns under each column. This indicates that different processes dominate albedo variability in different periods, which is consistent with the previous finding that the dominant periodicity of atmospheric waves in the Venus atmosphere varies on a timescale of several months (Del Genio & Rossow 1990; Kouyama et al. 2013; Imai et al. 2016).
A joint PCA of data in all windows combined together was also attempted. The periodograms of that study were dominated by peaks at about 220 days and its harmonics at 110, 55 days and so on, which are functions of the changing orbital orientation of the spacecraft relative to the Venus dayside, the observational window frequency, and other factors rather than interesting atmosphericphenomenon. Also, the dominance of these observational periodicities in the data caused peaks from the 4 and 5day waves to become statistically insignificant. As such, results from that analysis are not discussed in this paper since they are dominated by systematics.
Fig. 4 LombScargle periodograms for the loadings corresponding to the first six PCs for three windows at 283 nm. The numbers given for the spectral peaks indicate the associated periods. The solid red line is the average spectral power of red noise, while the dotted red line shows the 95% confidence limit of the distribution from 1000 Monte Carlo AR1 red noise processes. The explained variance fraction listed in the title of each subplot is calculated as the ratio of the eigenvalue of the mode to the sum of all eigenvalues from the covariance matrix, expressed as a percentage. 
Fig. 5 Same as Fig. 3 but for 365 nm. The patterns are much more variable between the windows, but can still be broadly classified into four groups. The first mode shows a pattern similar to a Hadley circulation, while the second shows a hemispheric oscillation. The third and fourth groups show signatures of the midlatitude Rossby waves and Yfeature, respectively. Not all groups are represented in all windows. 
3.3 Physical interpretations of the oscillations at 365 nm
The 365 nm PCs are comparatively much more variable across the three windows than the 283 nm PCs as shown in Fig. 5. However, they can still be approximately classified into a few broad categories. Two PCs in each window can be placed into the first category, which is similar to the latitudinal gradients seen at 283 nm. The periodicities, as seen in Fig. 6, are also similar to those seen at 283 nm. This is generally consistent with the understanding that the overturning circulation is responsible for increased UV absorber concentrations leading to low albedos at the equator at 365 nm (Titov et al. 2008; Molaverdikhani et al. 2012). However, latitudinal gradients are much weaker than those seen at 283 nm, suggesting that choatic variability in morphology is perhaps more important at this wavelength. Spatial distributions are in general far more complicated at 365 nm and have significant ambiguities in classification. This can be understood on the condition that the 365 nm observations are sensing an altitude level slightly below the cloud top, while 283 nm is at a higher altitude (Horinouchi et al. 2018). Therefore, the 365 nm features are affected by both the absorber and bright sulfuric acid cloud aerosols, which are formed through photochemical process (Mills et al. 2007; Parkinson et al. 2015), whereas the 283 nm would reveal the absorbers above the clouds so atmospheric flow patterns are more apparent. The latter would be also affected by photodissociation of SO_{2} (Mills et al. 2007) and the upper haze (Luginin et al. 2016), but our results suggest that these influences may be less effective than at 365 nm.
The second category is the hemispheric asymmetry, which is represented by two PCs each in the first two windows, but appears to be absent in the third. In window 3 PC 3 shows something like a hemispheric mode, however it is associated with a strong wave peak of around 4 days in Fig. 6, which is characteristic of the Yfeature. Window 3, PC 4 may also qualify for this category, but is somewhat ambiguous because the pattern is not purely hemispheric. The hemispheric mode typically does not show any strong periodicities apart from the spacecraft orbital period, as seen at 283 nm, although hints of the 4 and 5day periods are seen when the asymmetry is strong (window 1, PC 5, and window 2, PC 6). The large variability between windows could be attributed to concealment of gradients by transient cloud features, as for the previous category.
The third category shows a very clear transition region around 50° N/S, from dark lowtomiddle latitudes to the bright polar hood (Titov et al. 2012), which is sometimes symmetric across both hemispheres. The structure is associated with drastic changes in cloud tracked 365 nm winds, which slow down poleward of about 50° (Kouyama et al. 2012) in cloud top altitudes and also decrease poleward (Ignatiev et al. 2009), and in thermal structure, which present colder temperatures near the cloud top level (Tellmann et al. 2009). This means that winds, cloud top structure, and thermal structure are closely linked with the 365 nm patterns, implying overlapped vertical locations each other. The morphology bears a resemblance to model simulations of Rossby wave patterns at the cloud top (Kouyama et al. 2015). The periodicities always showa peak at around the 5day mark confirming the signature of midlatitude Rossby waves. This can be interpreted as features arising from the perturbation of the midlatitude bright bands or the “cold collar” (Titov et al. 2008, 2012) by Rossby waves. Interestingly, such a clear Rossby wave signature is not apparent in the 283 nm data.
The fourth category is the Yfeature, associated with a period of ~ 4 days with some other waves occasionally contributing at 4.5, 9, and 12.3 days. The absence of the Yfeature in window 2 does not mean that this feature did not occur during this time, since it is clearly seen at 283 nm for the same window. It maybe that the feature was sufficiently obscured at 365 nm that it is not represented in the first six PCs considered in this work.
4 Conclusions
We performed a PCA of the Akatsuki UVI 283 nm and 365 nm data taken over about 1.5 yr, subdivided into three observational windows each. The first six PCs over each window are considered significant and show similar morphologies at 283 nm, but are more variable at 365 nm. The difference can be understood if the 365 nm observations are sensitive to altitudes below the cloud top that areaffected by transient cloud variability, while 283 nm probes the atmosphere above the cloud tops. Additionally, since the unknown UV absorber is the result of unidentified chemical reaction chains, the kinetics of those reactions may also be responsible for some of the observed differences. The signatures of the overturning circulation, Rossby, and Kelvin waves are apparent from the spatial patterns and associated periodicities. We also note that similar spatial patterns are sometimes associated with different waves and the relative importance of different waves changes across the different observing periods as seen by the changing order of PCs with similar spatial patterns.
A hemispheric asymmetry mode is also apparent from this analysis. To be sure of its existence and to understand its dynamics, the same mode will need to be studied in other atmospheric observations of Venus. The analysis technique described in this work is general, has potential for use in other gridded atmospheric datasets from Akatsuki, Venus Express, and Pioneer Venus, and can be applied to other planetary atmospheric image analysis with a longterm monitoring dataset.
Acknowledgements
We thank Takehiko Satoh and the anonymous referee for comments that improved the paper. PK acknowledges generous funding support from the grantsinaid program of JSPS. Codes used in this study are available on request and Akatsuki data are publicly available at https://darts.isas.jaxa.jp/pub/pds3/staging/.
References
 Allen, M. R., & Smith, L. A. 1996, J. Clim., 9, 3373 [NASA ADS] [CrossRef] [Google Scholar]
 AlveraAzcárate, A., Barth, A., Rixen, M., & Beckers, J.M. 2005, Ocean Model., 9, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Arney, G., Meadows, V., Crisp, D., et al. 2014, J. Geophys. Res. Planets, 119, 1860 [NASA ADS] [CrossRef] [Google Scholar]
 Beckers, J.M., & Rixen, M. 2003, J. Atm. Ocean Technol., 20, 1839 [NASA ADS] [CrossRef] [Google Scholar]
 Bickel, P. J., & Levina, E. 2008, Ann. Stat., 36, 199 [CrossRef] [Google Scholar]
 Boyer, C., & Camichel, H. 1961, Ann. Astrophys., 24, 531 [NASA ADS] [Google Scholar]
 Burgess, T. M., & Webster, R. 1980, J. Soil Sci., 31, 315 [CrossRef] [Google Scholar]
 Cangelosi, R., & Goriely, A. 2007, Biol. Direct, 2, 2 [CrossRef] [Google Scholar]
 Carter, R. L. 2006, Res. Pract. Assess., 1, 4 [NASA ADS] [Google Scholar]
 Cattell, R. B. 1966, Multivariate Behav. Res., 1, 245 [Google Scholar]
 Del Genio, A. D., & Rossow, W. B. 1990, J. Atm. Sci., 47, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Dollfus, A. 1975, J. Atm. Sci., 32, 1060 [NASA ADS] [CrossRef] [Google Scholar]
 Dray, S., & Josse, J. 2015, Plant Ecol., 216, 657 [CrossRef] [Google Scholar]
 Encrenaz, T., Greathouse, T., Roe, H., et al. 2012, A&A, 543, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fukuhara, T., Taguchi, M., Imamura, T., et al. 2011, Earth Planets Space, 63, 1009 [NASA ADS] [CrossRef] [Google Scholar]
 Fukuhara, T., Futaguchi, M., Hashimoto, G. L., et al. 2017, Nat. Geosci., 10, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Garcia, D. 2010, Comput. Stat. Data Anal., 54, 1167 [CrossRef] [Google Scholar]
 Gilman, D. L., Fuglister, F. J., & Mitchell, Jr J. M. 1963, J. Atm. Sci., 20, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Horinouchi, T., Kouyama, T., Lee, Y. J., et al. 2018, Earth Planets Space, 70, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Hosouchi, M., Kouyama, T., Iwagami, N., Ohtsuki, S., & Takagi, M. 2012, Icarus, 220, 552 [NASA ADS] [CrossRef] [Google Scholar]
 Ignatiev, N., Titov, D., Piccioni, G., et al. 2009, J. Geophys. Res. Planets, 114 [Google Scholar]
 Ilin, A., & Raiko, T. 2010, J. Mach. Learn. Res., 11, 1957 [Google Scholar]
 Imai, M., Takahashi, Y., Watanabe, M., et al. 2016, Icarus, 278, 204 [NASA ADS] [CrossRef] [Google Scholar]
 Imamura, T. 2006, J. Atm. Sci., 63, 1623 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, D. A. 1993, Ecology, 74, 2204 [CrossRef] [Google Scholar]
 Khatuntsev, I., Patsaeva, M., Titov, D., et al. 2013, Icarus, 226, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Kondrashov, D., & Ghil, M. 2006, Nonlinear Process. Geophys., 13, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Kouyama, T., Imamura, T., Nakamura, M., Satoh, T., & Futaana, Y. 2012, Planet. Space Sci., 60, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Kouyama, T., Imamura, T., Nakamura, M., Satoh, T., & Futaana, Y. 2013, J. Geophys. Res. Planets, 118, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Kouyama, T., Imamura, T., Nakamura, M., Satoh, T., & Futaana, Y. 2015, Icarus, 248, 560 [NASA ADS] [CrossRef] [Google Scholar]
 Lebonnois, S., Sugimoto, N., & Gilli, G. 2016, Icarus, 278, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, Y., Imamura, T., Schröder, S., & Marcq, E. 2015, Icarus, 253, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, Y., Yamazaki, A., Imamura, T., et al. 2017, AJ, 154, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Luginin, M., Fedorova, A., Belyaev, D., et al. 2016, Icarus, 277, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Marcq, E., Encrenaz, T., Bézard, B., & Birlan, M. 2006, Planet. Space Sci., 54, 1360 [NASA ADS] [CrossRef] [Google Scholar]
 Marcq, E., Bézard, B., Drossart, P., et al. 2008, J. Geophys. Res. Planets, 113 [Google Scholar]
 Marcq, E., Bertaux, J.L., Montmessin, F., & Belyaev, D. 2013, Nat. Geo., 6, 25 [Google Scholar]
 Markiewicz, W., Titov, D., Limaye, S., et al. 2007, Nature, 450, 633 [NASA ADS] [CrossRef] [Google Scholar]
 Meinke, H., DeVoil, P., Hammer, G. L., et al. 2005, J. Clim., 18, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Mills, F. P., Esposito, L. W., & Yung, Y. L. 2007, Geophysical Monograph Series (Washington: AGU Publications) [Google Scholar]
 Molaverdikhani, K., McGouldrick, K., & Esposito, L. W. 2012, Icarus, 217, 648 [NASA ADS] [CrossRef] [Google Scholar]
 Mudelsee, M. 2002, Comput. Geosci., 28, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, S., & Freckleton, R. P. 2008, Trends Ecol. Evol., 23, 592 [CrossRef] [Google Scholar]
 Nakamura, M., Imamura, T., Ishii, N., et al. 2016, Earth Planets Space, 68, 75 [Google Scholar]
 Ogohara, K., Takagi, M., Murakami, S.y., et al. 2017, Earth Planets Space, 69, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Parkinson, C. D., Gao, P., Schulte, R., et al. 2015, Planet. Space Sci., 113, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Peralta, J., SánchezLavega, A., LópezValverde, M., Luz, D., & Machado, P. 2015, Geophys. Res. Lett., 42, 705 [NASA ADS] [CrossRef] [Google Scholar]
 Peralta, J., Iwagami, N., SánchezLavega, A., et al. 2018, Geophys. Res. Lett., 46, 1 [Google Scholar]
 PeresNeto, P. R., Jackson, D. A., & Somers, K. M. 2005, Comput. Stat. Data Anal., 49, 974 [CrossRef] [Google Scholar]
 Pourahmadi, M. 2011, Stat. Sci., 369 [CrossRef] [Google Scholar]
 Ross, F. E. 1928, AJ, 68, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Rossow, W. B., Del Genio, A. D., Limaye, S. S., Travis, L. D., & Stone, P. H. 1980, J. Geophys. Res. Space Phys., 85, 8107 [NASA ADS] [CrossRef] [Google Scholar]
 Rossow, W. B., Del Genio, A. D., & Eichler, T. 1990, J. Atm. Sci., 47, 2053 [Google Scholar]
 Sato, T. M., Sagawa, H., Kouyama, T., et al. 2014, Icarus, 243, 386 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, T. 2001, J. Clim., 14, 853 [NASA ADS] [CrossRef] [Google Scholar]
 Schulz, M., & Mudelsee, M. 2002, Comput. Geosci., 28, 421 [NASA ADS] [CrossRef] [Google Scholar]
 Tellmann, S., Pätzold, M., Häusler, B., Bird, M. K., & Tyler, G. L. 2009, J. Geophys. Res. Planets, 114 [Google Scholar]
 Titov, D. V., Taylor, F. W., Svedhem, H., et al. 2008, Nature, 456, 620 [NASA ADS] [CrossRef] [Google Scholar]
 Titov, D. V., Markiewicz, W. J., Ignatiev, N. I., et al. 2012, Icarus, 217, 682 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, G., Garcia, D., Liu, Y., De Jeu, R., & Dolman, A. J. 2012, Environ. Model. Softw., 30, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Warton, D. I. 2008, J. Am. Stat. Assoc., 103, 340 [CrossRef] [Google Scholar]
 Wilks, D. S. 2006, Statistical Methods in the Atmospheric Sciences, International Geophysics Series (Cambridge: Academic Press), 91 [Google Scholar]
 Yamamoto, M., & Takahashi, M. 2006, J. Atm. Sci., 63, 3296 [NASA ADS] [CrossRef] [Google Scholar]
 Yamazaki, A., Yamada, M., Lee, Y. J., et al. 2018, Earth Planets Space, 70, 23 [CrossRef] [Google Scholar]
All Figures
Fig. 1 Observation times of the 665 and 652 images used in this study for 283 nm (top left) and 365 nm (bottom left), respectively.The windows consist of 209, 244, and 212 images for 283 nm and 206, 237, and 209 images for 365 nm, respectively.The color of the data point indicates the fraction of the dayside that was observed in that image. Right: histograms showing the time gaps between successive images. The majority of the images are taken at intervals of 2 h, but longer intervals are not uncommon. 

In the text 
Fig. 2 Distribution of the 20 largest eigenvalues in the DINEOF interpolated datasets for 283 nm. The 365 nm datasets show similar trends. The red vertical line shows the point at which the eigenvalue curve has an inflection. The PCs corresponding to eigenvalues left of this line are considered significant by the Scree test. 

In the text 
Fig. 3 First six modes of variability from the DINEOF interpolated UVI dataset for the three windows at 283 nm. The first row is lined up in order of the PCs, while the second and third are rearranged so that similar patterns line up in the same column. The first mode shows a pattern similar to a Hadley circulation, while the second show a hemispheric oscillation. The third, fourth, and fifth columns appear to show a combination of short period atmospheric waves, while the sixth shows the well known Yfeature. 

In the text 
Fig. 4 LombScargle periodograms for the loadings corresponding to the first six PCs for three windows at 283 nm. The numbers given for the spectral peaks indicate the associated periods. The solid red line is the average spectral power of red noise, while the dotted red line shows the 95% confidence limit of the distribution from 1000 Monte Carlo AR1 red noise processes. The explained variance fraction listed in the title of each subplot is calculated as the ratio of the eigenvalue of the mode to the sum of all eigenvalues from the covariance matrix, expressed as a percentage. 

In the text 
Fig. 5 Same as Fig. 3 but for 365 nm. The patterns are much more variable between the windows, but can still be broadly classified into four groups. The first mode shows a pattern similar to a Hadley circulation, while the second shows a hemispheric oscillation. The third and fourth groups show signatures of the midlatitude Rossby waves and Yfeature, respectively. Not all groups are represented in all windows. 

In the text 
Fig. 6 Same as Fig. 4 but for 365 nm. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.