On the temporal variability of the surface solar radiation by means of spectral representations

. This work deals with the temporal variability of daily means of the global broadband surface solar irradiance (SSI) impinging on a horizontal plane by studying a decennial time-series of high-quality measurements recorded at a BSRN ground station. Since the data have a non-linear and non-stationary character, two time-frequency-energy representations of signal processing are compared in their ability to resolve the temporal variability of the pyranometric signal. First, the continuous wavelet transform is used to construct the wavelet power spectrum of the data. Second, the adaptive, noise-assisted empirical mode decomposition is employed to extract the intrinsic mode functions of the signal, followed by Hilbert spectral analysis. In both spectral representations, the temporal variability of the SSI is portrayed having clearly distinguishable features: a plateau between scales of two days and two-three months that has decreasing power with increasing scale, a large spectral peak corresponding to the annual variability cycle, and a low


Introduction
The surface solar irradiance has been recognized by the Global Climate Observing System as one of the Essential Climate Variables, a restrained set of critical variables extensively used in science and policy circles for understanding climate evolution and guiding mitigation and adaptation efforts (Bojinski et al., 2014). As such, not only do long term time-series of the SSI play a key role in science and policy consulting, but they also strongly impact the engineering and financial sectors, providing support for, e.g., the selection of suitable sites for solar power plants or investment decisions, respectively (Schroedter-Homscheidt et al., 2006). Thus, better understanding the characteristic time-scales of variability of the solar radiation at ground-level, as recorded in the aforementioned long term time-series, is one of the goals of the present work.
However, in order to tackle the topic of the variability of the surface solar irradiance, the very concept of variability itself must first be properly addressed. Even though spatial variability is of course also of interest, as has been shown for example by Hoff and Perez (2010), in line with the scope of the study, only the temporal dimension of SSI variability will be considered.
One of the most common ways of quantifying variability is through the use of simple descriptive statistics, such as the mean and the standard deviation (Hammer and Beyer, 2012). Variability is then usually related to the time-scale at which analysis is carried out, with the fluctuations at the scale of interest being compared to the mean (or another low-pass filtered approximation) at the next higher scale. Thus, yearly sums of SSI are compared to their multiannual mean, daily sums to their monthly mean, minute values are compared to their hourly mean, etc. When, still, even finer (sub-second) scales are addressed, the second temporal derivative of irradiance is used to quantify its instantaneous variability (Yordanov et al., 2013). Another way to express SSI variability by using statistics is in terms of the probability density function of the irradiance, or by using a two-sample structure func-Published by Copernicus Publications. tion, which is the root-mean-square deviation of increments of irradiance, as described by Skartveit and Olseth (1992).
While adhering to the statistical approach and defining the root mean square of a high-pass filtered version of the data as the fluctuation factor of the time-series, Otani et al. (1997) allude to employing the power spectral density as a mean of describing variability. Hence, spectral methods provide a different approach of describing variability in that, by analysing the frequency contents of the data, a much clearer picture is rendered in terms of the scales, or the inverse of frequency, at which variations in the SSI occur. As such, Woyte et al. (2007) make use of wavelet periodograms to characterize fluctuations in solar irradiance, while Lave et al. (2013) use the wavelet transform to decompose a time series of SSI for their proposed Wavelet Variability Model. Calif et al. (2013) make use of the power spectrum of a time-series in order to model the multi-fractal properties of the SSI by means of, among other spectral methods, Hilbert analysis.
Having seen the variety of techniques used to describe the temporal variability of the SSI, the present work investigates the characteristic time-scales of variability embedded in long-term time-series of SSI, as revealed by two spectral methods in particular. Given the nonlinear and nonstationary character of solar radiation measurements (Zeng et al., 2013), the Hilbert-Huang Transform (HHT) is a natural choice since it has been specifically designed for the analysis of such data (Barnhart and Eichinger, 2011). The spectral picture rendered by the HHT is compared to the one obtained by the now classical Continuous Wavelet Transform (CWT) (Torrence and Compo, 1998), highlighting the commonly shared and contrasting features of the two representations. Since these methods are applied to long term time-series of the SSI, the sampling rate of the data has been fixed at one day.
The aims and intents of the study at hand can now be described as ascertaining what type of information about the temporal variability of the data the two spectral representations, the HHT and the CWT, may convey. Added to this, it is well known that the SSI, more specifically a decennial timeseries thereof, exhibits a characteristic time-scale of variability around a one year period which is easily accounted for by the orbital parameters of the Earth-Sun system. However, other factors controlling surface solar irradiance are also at play, such as atmospheric absorption or scattering, due to e.g. clouds, aerosols, etc. Thus, an ancillary locus of investigation is whether such processes also induce other characteristic scales of variability that may be embedded in the scrutinized signal.
In what follows, this study is organized as per the ensuing outline. Section 2 discusses the data and the analysis methods; the time-series under investigation is described in Sect. 2.1, the spectral methods employed are reviewed in Sect. 2.2. The obtained results are then presented in Sect. 3, while discussion thereof is deferred to Sect. 4. Finally, conclusions and outlook are presented in Sect. 5.

The analysed time-series and preprocessing
The time-series under investigation in this study was obtained from high-quality ground measurements of the SSI, at a station situated in Carpentras (44.083 • N, 5.059 • E), in southeastern France. The measuring station is part of the Baseline Surface Radiation Network (BSRN, König-Langlo et al., 2015). BSRN is a worldwide radiometric network providing accurate readings of the SSI at 1 min temporal resolution and with an uncertainty requirement at 5 W m −2 for SSI (Ohmura et al., 1998). The original 1 min measurements for the ten-year period 2001-2010 have been quality checked following Roesch et al. (2011). Next, daily means were computed from this raw data only if more than 80 % of the samples during daylight were valid. Isolated missing daily means were completed by linear interpolation applied to the daily clearness index, which is the ratio between the daily mean of SSI and the daily mean of the corresponding total solar irradiance received on a horizontal surface at the top of atmosphere. The resulting time-series is shown in Fig. 1.

Time-frequency-energy representations
For this study, two spectral methods have been employed in order to assess the variability of the SSI. Most notably, analysis techniques that reveal the temporal evolution of the scrutinized data have been privileged. Two such techniques, popular within the geophysics community (Tary et al., 2014), are the Continuous Wavelet Transform (CWT) an the Hilbert-Huang Transform (HHT). The methods are of interest mostly due to their innate ability to offer insight into the processes underlying the data by painting so called time-frequencyenergy representations, i.e. graphical representations depicting the distribution of power, or variance, occurring at specific time instants and specific scales or periods (assuming period 1/frequency).

The continuous wavelet transform (CWT)
The continuous wavelet transform is a mean to render a two parameter representation (time and scale, or frequency) of a signal of interest in terms of wavelets obtained by timeshifting and dilating a basic wavelet of fixed shape. Starting from this latter, elementary "mother" wavelet, any signal can be decomposed by measuring its degree of similarity to wavelets that are time-shifted and scaled versions of the original one. Consequently, this scaling of the basis functions enables the wavelet transform to achieve high frequency resolution when analysing slowly varying features and high temporal resolution when investigating bursts and transients.
To formalize this, consider a basic, or mother, wavelet, ψ(t), from which, by scaling with some factor λ ∈ R, λ > 0, and time-shifting by τ ∈ R, a family of wavelets can be defined: Then, the CWT can be regarded as measuring the similarity between a signal of interest, x(t) and the basis functions ψ λ,τ (t), as in (Rioul and Vetterli, 1991): From Eq.
(2) the wavelet transform is seen to be linear, since the analysis consists of computing the inner products between the wavelets and the signal of interest (Farge, 1992). Another point to be made is that the CWT is overly redundant, hence its amenability to quantitative analysis is arguable .
Since, due to its good trade-off between time and frequency resolution, the Morlet wavelet is frequently chosen as such a basic wavelet, following Kumar and Foufoula-Georgiou (1997), this study has also settled on the Morlet wavelet as "mother" wavelet, given in Eq. (3), where the parameter ω 0 ≥ 5.
The analysing wavelet is in general, as also in particular for Eq. (3), a complex-valued function; it follows that the continuous wavelet transform is also complexvalued. Hence, the wavelet power spectrum is defined as the square of the complex amplitude of the wavelet coefficients, i.e. |CWT x (λ, τ )| 2 , yielding a time-frequency-energy representation (Torrence and Compo, 1998). Nevertheless this power spectrum has been found to be biased towards large scales, thus rectification of the coefficients through division by the scale is undertaken in order to enable the comparison of the spectral peaks across scales (Liu et al., 2007).

The Hilbert-Huang transform (HHT)
In contrast with the CWT, which decomposes data into a set of a priori chosen basis functions, the Hilbert-Huang transform is completely data-driven and adaptive. In other words, the transform does not "measure" the similarity between the data and some predefined "rulers", such as trigonometric functions or wavelets, in the case of the Fourier analysis and the CWT, respectively. Rather, the method lets the signal itself be the driver of the decomposition. Furthermore, the HHT was designed with the analysis of non-linear, non-stationary signals in mind (Huang et al., 1998), closely resembling a previous proposal of non-linear analysis of non-stationary frequency-amplitude structure of long timescale solar activity (Nagovitsyn, 1997). Thus, owing to its ability of handling data issued from the non-linear interaction of physical processes, which are often under the influence of non-stationary external forcings , the method has seen extensive use in geophysical research (Huang and Wu, 2008).
The HHT consists of two separate steps. First, the oscillations embedded in the data, called Intrinsic Mode Functions (IMFs) are extracted by means of an algorithmic procedure, termed the Empirical Mode Decomposition (EMD), whose outline is sketched in algorithm (1). After all the IMFs have been extracted, the residue or trend, which cannot be mathematically thought of as an oscillation, is all that is left of the time-series.

t) {Return IMFs and residue}
By construction, the EMD extracts oscillatory modes that have a well-behaved Hilbert transform (Huang et al., 1998), by means of which the instantaneous amplitude and frequency of each mode is computed as follows. Thus, following the notation from algorithm (1), for a real-valued IMF c k (t), its Hilbert transform can be written as: where P indicates the Cauchy principal value. Next, each IMF and its Hilbert transform are used to compute the complex-valued analytic signal: are the instantaneous amplitude and instantaneous phase, respectively. The instantaneous frequency can be then be thought of as the derivative of the instantaneous phase: The original data can now be written in terms of this amplitude modulation-frequency modulation (AM-FM) signal model: which is essentially a sum of zero-mean oscillations having symmetric envelopes that are riding onto the EMD trend. The second part of the method, Hilbert spectral analysis, is next used to construct the Hilbert energy spectrum. Huang et al. (2011) define this time-frequency-energy representation as "the energy density distribution in a time-frequency space divided into equal-sized bins of t × ω with the value in each bin summed [. . . ] at the proper time, t, and proper instantaneous frequency, ω".
For the purpose of this study, a modified version of the algorithm, the Improved Complete Ensemble Empirical Mode Decomposition, introduced by Colominas et al. (2014), is used, as it allows for an exact decomposition of the data and is more robust with respect to noise. A fast EMD routine, provided by Wang et al. (2014), has been employed to speed up computations.

Results
The continuous wavelet transform was applied to the timeseries of SSI, producing the wavelet power spectrum depicted in Fig. 2. It is noteworthy that the abscissa is the same as for the original data, from Fig. 1, providing the means to temporally resolve the various features of the spectrum. The ordinate, having a logarithmic scale this time, indicates the approximate period, expressed in days. Lastly, colour encoding is used to logarithmically express the amount of power Adv. Sci. Res., 13, 121-127, 2016 www.adv-sci-res.net/13/121/2016/ that the time-series contains at a certain point in time and at a certain period or scale. The colour bar at the top of the wavelet spectrum provides a scale assigning numerical values to the colour encoding. The white-out regions that run along the edges of the figure indicate the so-called cone of influence (COI), an area that is awash with edge effects which are a principal limitation of the CWT and that should be excluded from analysis. The COI edges correspond to the power of an impulse at the edge dropping by a factor e −2 (Torrence and Compo, 1998). As the scale increases, the COI progressively widens until the two edge-bound regions intersect in the middle of the spectrum, thus imposing a limitation with respect to the highest period at which the wavelet power spectrum is able to provide meaningful information.
The global wavelet spectrum in the right panel is the timeintegrated version, i.e. line-by-line sum, of the wavelet spectrum and indicates the amount of power at each scale, expressed in dB. Secondly, the empirical mode decomposition was used to extract IMFs from the data. By applying Hilbert spectral analysis to them, another time-frequency-energy representation of the same data is obtained, as illustrated in the left panel of Fig. 3. The axes of the Hilbert energy spectrum are analogous to the axes of the wavelet spectrum in Fig. 2, except for the colour coding, which now spans a different range. Edge effects are usually contained within a half-period of a component at data boundaries , an area that has been whitened out in similar way to wavelet COI from

Discussion
The wavelet power spectrum in Fig. 2 reveals that the most energetic component of the SSI time-series is confined between scales of approximately 256 and 512 days, as indicated by the thick orange power-band near the bottom of the plot. In CWT, spectral smearing and leakage occur because of the finite support of the analysing basic wavelet (Huang et al., 1998). Power "spills over" to neighbouring scales and is assigned to this wider orange band on the plot instead of a distinct peak. This is also apparent in the support of this feature in the global wavelet spectrum. This feature can be interpreted in terms of the yearly seasonal variation of the SSI with a period of around 365 days; as the amount of sunlight reaching the ground is markedly different between summer and winter, the apparent cycle has a strong power contribution at this scale. Nevertheless, this band also contains a bright yellow strip centred on 365 days, which is thinning beginning in 2003 and reaching the smallest extent in scale space around 2008, indicative of a minimum in the eleven year cycle solar cycle (Hathaway, 2015). Is is also interesting to note, that because of the linear nature of the CWT, non-linearities of the yearly cycle are also identifiable in the green band at the scale of about 180 days (the second harmonic of the fundamental yearly cycle), which is dotted with blue patches that indicate low power during boreal winter. The power of this feature is such that it also makes a slight impression on the global wavelet spectrum. Furthermore, the grouping of yellow-orange "filaments" at scales of between 64 and 2 days, that correspond to an increase in SSI during boreal summer, could also be interpreted as high-order harmonics. However, the same feature may as well be indicative of amplitude modulation of these bands by the yearly cycle (similar to Liu et al., 2007, Fig. 1), since this cycle can be inferred just from the alternation of blue-green (winter lows) and orange-yellow (summer highs) hues in this band. Power in this high-frequency band decreases with increasing scale, as can be inferred from the global wavelet spectrum, where a minima is also noticeable at around 128 days. This minima, centred in a larger depression, is a consequence of the low-powered regions occurring in the 150 to 64 days band in the left panel, denoted by dark blue patches, most notably between the second half of 2003 and 2005, the second half of 2005 and 2007, and from the second half of 2009 onward. These dark blue patches correspond to minimal levels of correlation between the signal and the analysing wavelets at those scales and temporal localizations. A final remark must be made here by reiterating that owing to the overly redundant nature of the CWT, even after rectification, the global spectrum is best interpreted qualitatively and not quantitatively.
The Hilbert spectrum depicted in Fig. 3 shows similar features. As in Fig. 2, the most energetic component is found oscillating around the one year period, in the form of a yellow trace, with the notable exception that in this representation the frequency modulation of the component can be clearly distinguished, i.e. there are slight downward period shifts, with respect to 365 days, during winter and upward shifts during summer. Also of notice is the fact that the domain of scale modulation of the yellow trace is somewhat restrained with respect to the thick power band between 256 and 512 days from Fig. 2, as can also be inferred from the support of this feature on the marginal spectrum where it accounts for the greatest power peak. The decrease in frequency modulation between 2007 and 2009 corresponds, once again, to a period of solar activity at its minimum. Furthermore, two low power components can be identified as having periods between 4 years and 1 year, with two distinct spectral peaks showing on the marginal spectrum at about 1200 and 800 days, respectively. The physical meaning of these peaks is ambiguous since most of the power of these components lies in the COI area. Another feature of the spectrum that stands out is the lack of power in the 2-3 months to one year band, which is best viewed on the marginal spectrum. This was also hinted at by the depression at about the same scales on the global wavelet spectrum in Fig. 2, however since the EMD does not detect any spectral features in this region, the Hilbert spectrum, owing to sparsity, assigns a negligible amount of power in this band. Another feature is the band of quasi-constant power in the high-frequency region, between 2 days (the Nyquist limit) and 2-3 months. As already indicated by the CWT global spectrum in Fig. 2, here too the power of this plateau decreases with increasing scale. Modulation in amplitude of this band by the yearly cycle is apparent here too, and is most visible for the years 2005 and 2007 where the blue hues occurring during winter turn green-yellow in the high-SSI regime of summer. A final note regarding the Hilbert spectrum is that binning, both in the colour and in the scale space (see Sect. 2.2.2), may yield some features as continuous lines while others are rendered as point-like, especially where rapid frequency modulation takes place, such as the high-frequency plateau.

Conclusion and outlook
To sum up, it has been shown that both the CWT and the HHT are able to capture the characteristic time-scales of variability of the data. However, to our opinion, the HHT is able to pick up more subtle features, such as the yearly period modulation, while also yielding a somewhat more sparse representation. This makes it a useful tool for the analysis of the temporal variability of the surface solar radiation in forthcoming studies. Nevertheless, both transforms provide a time-frequency-energy representation of the signal where it can be clearly seen that the quasi-annual oscillations have the highest energetic contribution. A power depression in the spectrum can also be identified, which is a particular feature of the analysed data and may prove useful in modelling the SSI for this specific site, i.e. any such model should avoid injecting power into this frequency band, or filter it out at the very least. A slightly slanted plateau between 2 days and 2-3 months can also be identified, whose power declines towards lower scales. Preliminary investigations suggest that this high-frequency region is noise-like and modulated by the yearly cycle, and is an indication of cross-scale coupling of phase-amplitude (Paluš, 2014;Huang et al., 2016); further research is needed for a definitive answer. If this plateau is noise-like, it would impact numerical modelling of the SSI in the context of forecasting, as it would imply that, for this specific site at Carpentras, there exists a deterministic yearly component, a gap in the spectrum at medium scales, and a stochastic component modulated by the highly energetic yearly cycle.