Irradiance Variability Quantiﬁcation and Small-Scale Averaging in Space and Time: A Short Review

: The ongoing world-wide increase of installed photovoltaic (PV) power attracts notice to weather-induced PV power output variability. Understanding the underlying spatiotemporal volatility of solar radiation is essential to the successful outlining and stable operation of future power grids. This paper concisely reviews recent advances in the characterization of irradiance variability, with an emphasis on small spatial and temporal scales (respectively less than about 10 km and 1 min), for which comprehensive data sets have recently become available. Special attention is given to studies dealing with the quantiﬁcation of variability using such unique data, the analysis and modeling of spatial smoothing, and the evaluation of temporal averaging.


Introduction
The last decade has seen a noticeable increase in the installed capacity of photovoltaic (PV) power all over the world, and the corresponding fast-paced growth (e.g., 25% in 2015 [1]) is expected to continue in the future. From well over 200 GW in the beginning of 2016 [2], the global PV capacity is estimated to multiply more than tenfold until 2030, with capacity projections ranging between 3 TW and 10 TW [3]. As a result, problems related to the intrinsic variability of PV power production will considerably increase as well [4,5]. These challenges include the correct estimation of a PV system's yield [6], the proper dimensioning of energy storage [7], the balancing of generation and load [8], as well as the support of power quality (e.g., voltage and frequency stability) [9]. As PV power variability is mainly determined by weather-induced heterogeneity in fields of downwelling solar radiation [10], corresponding data-driven analyses and irradiance variability quantifications are essential to the successful outlining and stable operation of future power grids [11].
In this context, variability in irradiance itself is as interesting as variability in irradiance increments-that is, transitions from one point in time to another, also known as changes (e.g., [12]), step changes (e.g., [13]), or ramp rates (e.g., [14]). Irradiance variability mainly impacts a PV system's yield and the proper dimensioning of energy storage, while increment variability affects power quality as well as the maintenance of the generation load balance. The relevant spatiotemporal scales of variability span many orders of magnitude: from meters and seconds through hundreds of kilometers and multiple days, depending on the PV capacity and power grid size in question [15]. There is an ongoing need to deepen the understanding of misrepresentation of temporal variability caused by temporally coarse-resolution measurements [16,17], as well as how spatial smoothing (resulting from many PV systems dispersed over a large area) affects variability [18,19]. This is especially true for small sub-minute and sub-kilometer scales, which have only begun to receive increased scientific attention in the last couple of years, and for which comprehensive data sets have recently become available [10,12,15].
With an emphasis on the three last-mentioned research topics and on sub-minute scales below about 10 km, this short review article summarizes previous findings from the literature and concisely recapitulates the essentials of characterizing normalized time-scale-specific changes in irradiance (Section 2), as well as quantifying variability (Section 3) and averaging effects in space (Section 4) and time (Section 5). The conclusions finally provide a short outlook onto possible further research questions in light of the state-of-the-art (Section 6). The paper complements previous textbooks and literature reviews on solar energy forecasting and resource assessment [119][120][121][122][123][124][125][126][127][128], the estimation of variability in different renewable resources [15,129], the use of machine learning in solar power applications [130][131][132], and small-scale structures of irradiance fields [133].

Irradiance Normalization and Time-Scale-Specific Changes
It is well-known that variability in global horizontal irradiance (GHI) can generally be influenced by processes of both astronomical and atmospheric origin [119]. On the one hand, astronomical relationships cause seasonal and diurnal variations, which are accurately computable and relatively small on short sub-minute time scales. On the other hand, atmospheric phenomena can affect irradiance in complicated ways in a range of time scales (e.g., cloud dynamics can exert influence on the seasonal cycle as well as sub-minute variability). Some studies aim to focus on the actual magnitude of irradiance including all aforementioned influences, but occasionally employ a normalization to irradiance under standard test conditions (STC) G STC = 1000 Wm −2 [26,50]. In order to study the stochastic nature of weather-induced variability on short time scales, most authors take measures to remove deterministic trends from GHI time series G by estimating either the clearness index (a.k.a. transmission coefficient, e.g., [20]) or the clear-sky index The clearness index represents a normalization of G to the extraterrestrial solar radiation G extra (i.e., irradiance at a particular location without considering the atmosphere), while the clear-sky index relates G to a theoretical clear-sky radiation G clear (i.e., irradiance at a particular location assuming a cloud-free atmosphere). The extraterrestrial solar radiation G extra is exclusively determined by well-known astronomical interrelations, whereas the characteristics of G clear depend on parameters of atmospheric conditions and are thus model-specific. Many different clear-sky models have been proposed to date [134][135][136][137], and the use of any of them introduces unique uncertainties to the clear-sky index time series, which are not included in the source GHI observations. With G clear thus not being unambiguously defined, its values represent typical rather than effective clear-sky irradiance. Nevertheless, many authors have favored the clear-sky index over the clearness index, although nomenclature varies and some call k * as defined in Equation (2) "clearness index" (e.g., [100]). Figure 1 presents (a) an example time series of surface irradiance recorded with 1 s temporal resolution during the HD(CP)2 Observational Prototype Experiment (HOPE) campaign [138,139] along with simulated clear-sky irradiance derived by a basic model [140], and (b) estimated clear-sky index as per Equation (2). Several typical features of the clear-sky index are evident in the time series. For example, the lowest clear-sky index values are usually greater than zero, because irradiance is never attenuated completely, no matter how dark the cloud. Moreover, the upper boundary of k * often exceeds 1, mostly due to cloud enhancement (i.e., reflections from surrounding clouds add to the measurement) and to a lesser degree due to imperfectly modelled clear-sky conditions. Under broken cloud conditions, cloud enhancement has been reported to lead to single-point irradiance observations exceeding their clear-sky expectation by more than 50% on sub-minute time scales [17,64,65,[141][142][143][144]. However, cloud enhancement cannot be unambiguously identified using absolute values of k * , because these can be biased by the clear-sky model of choice. On these short time scales, the probability density function of clear-sky index time series is typically bimodal in nature, with the two peaks respectively representing states with and without cloud coverage [22][23][24]145,146]. Figure 2 correspondingly shows two estimates of the clear-sky index probability density function based on a histogram of the example time series shown in panel (b) of Figure 1, as well as a kernel density estimate (KDE) [147] of all available data from that sensor (about 4.5 · 10 6 s worth of 1 s data). While the histogram does not correspond to a perfectly bi-modal distribution due to the limited number of data available from the single-day example, the kernel density estimate clearly illustrates the two peaks representing overcast and clear skies.  Figure 1 (histogram, black lines), as well as a considerably longer series containing a total of about 4.5 · 10 6 s worth of single-sensor 1 s measurements collected during the HOPE campaign [138,139] (kernel density estimate, KDE, dashed blue line). The KDE was derived using Gaussian kernels with a smoothing bandwidth (i.e., smoothing kernel standard deviation) of 0.05.
The statistical moments associated with a given clear-sky index time series of length T, such as the sample arithmetic mean and the sample standard deviation can be used to characterize the time series' variability to a degree [51]. However, these measures do not consider the chronology of the observations, and are thus not suited to assess changes therein (i.e., temporally shuffling all clear-sky index data will affect neither the respective probability density function nor the corresponding moments). Instead, characterizations of clear-sky index variability can be geared towards changes over specified intervals of time τ by deriving statistics of k * increments which are a useful measure of intermittency [148].   (5) for three distinct short-term increment time steps τ = 5 s, τ = 10 s, and τ = 60 s. The underlying irradiance data were collected during the HOPE campaign [138,139].
The probability density functions of single-point clear-sky index increments generally exhibit a slim central peak (representing a high probability of occurence from relatively small changes), enclosed by widespread tails of decreasing probabilities from stronger variations. For increasing increment time steps, the numbers of high-magnitude changes also increase, leading to more pronounced tails of the corresponding distributions [78]. These general characteristics are illustrated for the same three time steps as before (τ = 5 s, τ = 10 s, and τ = 60 s) by means of kernel density estimates in Figure 4, using the same data as for the KDE in Figure 2 (i.e., about 4.5 · 10 6 s worth of 1 s clear-sky index data).
Density τ = 5 s τ = 10 s τ = 60 s Figure 4. Kernel density estimates of clear-sky index increments, using about 4.5 · 10 6 s worth of single-sensor 1 s data collected during the HOPE campaign [138,139]. The estimates were derived for increment time steps τ = 5 s, τ = 10 s, and τ = 60 s using Gaussian kernels with a smoothing bandwidth of 0.01.
Comparisons of spatially averaged increment distributions of multiple pyranometers (or PV plants of different capacities) indicate that the numbers of high-magnitude changes and, thus the width of the increment distribution, are reduced for an increasing number of sensors (or larger PV plants) [26,48,56]. However, higher-than-normal (i.e., compared to a Gaussian distribution) probabilities of comparatively high-magnitude changes remain common across time scales (e.g., increments on the order of tens of standard deviations can be recorded on a regular basis [44,104]). Similarly, temporal averaging on scales 1 s can also result in a narrowing of increment distributions, and hence a potential underestimation of actual variability [111]. Yet, this drawback can be acceptable, for example, when variability at different locations is compared by means of data featuring diverse temporal resolutions [50].

Variability Quantification
Several different single-number metrics have been proposed to quantify variability in irradiance, normalized irradiance, and/or PV power (some of the underlying concepts can be applied to several of these quantities) within a well-defined period. These include specialized measures such as the variability index (VI) [96], which is tailored specifically to employ measured irradiance G and simulated clear-sky irradiance G clear , and is defined as with T denoting the number of consecutive observations and ∆T representing the averaging interval of the time series (e.g., ∆T = 1 for minute-averaged irradiance measurements). In contrast, the intra-hour variability score (IHVS) [98] and the daily aggregate ramp rate (DARR) [26] have also been proposed to be used with irradiance data, but can easily be adapted to other quantities, because they simply represent the sum of all absolute increment values during a given hour or day, respectively. For example, the DARR of any time series X can be defined as with T again denoting the number of consecutive observations and C representing an (optional) scaling constant (e.g., for irradiance, C = 1000 Wm −2 can be used to normalize to STC [26]). The variability score (VS) [50] is another versatile measure, which can be defined for any time step τ as with ∆X τ = X(t + τ) − X(t) denoting all increments of X using time step τ (i.e., analogous to Equation (5)), ∆X 0 consisting of values between 0 and max(|∆X τ |), and P(|∆X τ | > ∆X 0 ) representing the probability of absolute increment values greater than ∆X 0 to occur. In its original formulation, the VS is calculated based on increments in temporally averaged data, with averaging time scale ∆T = τ, and it is expressed as a percentage, with ∆X τ and ∆X 0 representing either increments in irradiance normalized to STC or increments in rated PV power capacity [50]. However, the concept of the variability score can generally also be applied to other time series, such as clearness index or clear-sky index. The variability score compares well to the variability index, as both yield comparable estimations of variability with a very high linear correlation between the two measures for a given time scale [16,50]. Finally, changes in the variability of any increment time series ∆X τ are also reflected by changes of the corresponding standard deviation with T denoting the number of increment values in the series, and ∆X τ representing the mean increment values (in practice, ∆X τ → 0). Increment standard deviation σ ∆X τ characterizes typical excursions from the mean, and has become a well-established measure to quantify variability in irradiance, clear-sky index, and PV power output [10,15,19,28,51,111]. For example, the standard deviation is particularly convenient when considering averaged k * time series across multiple pyranometers, as it will change as n −0.5 for a number of n uncorrelated locations, regardless of the underlying distribution [51]. Note, however, that the non-Gaussianity of irradiance-related increment statistics causes the standard deviation to not necessarily represent the size of extreme fluctuations appropriately. For example, the universal three-sigma rule, stating that a range of ± 3 σ ∆X τ around the mean will cover 99.73% of all values if ∆X τ are normally distributed [149], does not apply to increments in irradiance, clear-sky index, or PV power. Figure 5 compares daily values of the clear-sky index variability score VS k * τ (cf. Equation (8)) and the clear-sky index increment standard deviation σ ∆k * τ (cf. Equation (9)) for three time steps τ = 5 s, τ = 10 s, and τ = 60 s. Although increasing standard deviations generally correspond to increasing variability scores, there is considerable spread around the least-square fits shown in Figure 5. Consequentially, the coefficients of determination 0.60 R 2 0.65 associated with the fits only indicate a moderate correlation between the two quantities considered. As both measures are actively being used to quantify variability throughout the literature, it may be worthwhile to systematically compare the two in a more detailed manner in order to improve the interpretability of their values. However, such an analysis would go beyond the scope of this review.  (8)) and the clear-sky index increment standard deviation σ ∆k * τ (cf. Equation (9)) for three time steps (a) τ = 5 s, (b) τ = 10 s, and (c) τ = 60 s, using 105 days worth of single-sensor 1 s clear-sky index estimates derived from irradiance data collected during the HOPE campaign [138,139] (i.e., the same data as in Figure 4). Least-square fits are included as dashed lines, and the corresponding coefficients of determination R 2 are quoted in each panel.
Additionally, a number of studies have applied classification schemes in order to group sets of time series according to their variability characteristics. Typically, the resulting variability classes are made up of at least "high-variability", "medium-variability", and "low-variability" [97], although the nomenclature varies (e.g., high variability can be called "mixed-sky", medium variability "clear-sky", and low variability "overcast" [51]). Some authors split the range of a single-number metric into subsets using fixed thresholds [101], while others employ multivariate classification schemes [51, 100] or machine learning algorithms [97,99].

Spatial Averaging
Compared to a single-point irradiance measurement, the larger the panel-covered area of a PV plant, the less variable its power production becomes [26], especially on short sub-minute time scales [55,56]. Similarly, the power output of a single high-capacity PV plant occupying a relatively small area is more variable than the aggregate power of a large number of small plants amounting to a similar total capacity, but occupying a considerably larger area [28]. In lieu of suitable PV power data, averages of irradiance or clear-sky index increments can serve as output variability estimates of a set of PV power plants located at different locations [10].
An illustrative example of spatial averaging is shown in Figure 6, where (a) single-sensor clear-sky index estimates are contrasted with the spatially averaged values of up to 99 synchronized pyranometers dispersed over an area of approximately 80 km 2 , and (b) 1 min increments are shown for both the single sensor and the spatial average. Pronounced spatial smoothing is evident in both panels. Single sensor Spatial average Figure 6. An example of spatial averaging: (a) the same 1 s single-sensor clear-sky index k * time series as in Figure 1 along with a corresponding spatially averaged 1 s clear-sky index computed over as many as 99 pyranometers dispersed over an area of about 80 km 2 ; (b) 1 min clear-sky index increments ∆k * τ derived according to Equation (5) for both the single-sensor time series and the spatial average. The underlying irradiance data were collected during the HOPE campaign [138,139].
As a characterization of the spatial structure of ∆k * τ fields, it is useful to consider spatial two-point correlation coefficients between locations i and j: which govern the process of spatial averaging (i.e., the standard deviation σ x+y of the sum of two random variables x + y with correlation coefficient ρ xy is σ x+y = σ x 2 + σ y 2 + 2ρ xy σ x σ y . In Equation (10), ∆k * τ,i (t) and ∆k * τ,j (t) are the two individual increment time series at the two locations i and j, while ∆k * τ,i and ∆k * τ,j are the corresponding arithmetic means (computed as in Equation (3)), and the quantity T denotes the number of data points in each of the two time series. These correlation structures can, for example, serve to directly model the generation variability of distributed and/or large PV plants based on a single pyranometer, using a wavelet approach [48,58,83,97].
Increment correlations (in both time and space) have been shown to depend on effective cloud speed [10,28], different daily variability categories [13,97], and short-term sky types [51]. Moreover, correlation coefficients can be influenced by the sensor pair orientation relative to the direction of cloud motion [49,52,68,71,103]. For all measures and approaches considered, increment correlation structures consistently decrease with increasing pair distance, and the rate of decrease becomes smaller (longer decorrelation distances) when considering larger-increment time steps. Along-wind increment correlation structures have also been suggested to feature negative peaks, with correlation coefficients reaching values between −0.5 < ρ ∆k * τ ij < 0 [12,13,49,70,71].
Several methods have been suggested to simulate increment correlation structures for specific time scales. Many authors have employed empirical fits to measured data [10,13,70,97,103,104,106], while others used simplified cloud shapes generated by randomly positioned and subsequently blurred squares [70], unions of randomly positioned discs [68], double-sigmoidal clear-sky index deviations [71], or fractal structures [52]. For example, k * increment correlation structures were extracted from hourly satellite-derived irradiance data, with pair distances in the range of 10 km ≤ d ij ≤ 300 km [10]. Taking increment time lags of 1 h ≤ τ ≤ 4 h and estimated relative cloud speed v into account, the correlation structure was estimated as  [70,83]. These models were obtained as curve fits from different data with relatively coarse spatiotemporal resolutions.    where differences between models range between −0.5 ρ ∆k * τ ij 0. For further increasing distances d ij > 100 m, each curve transitions to its own unique character of constant or increasing correlation coefficients, approaching 0.

Temporal Averaging
Similar to spatial averaging, irradiance variability is reduced when considering temporally averaged data. Figure 8 presents examples of such temporal averaging, based on the example period employed in the previous figures, using two different averaging periods of 1 min and 10 min. Panel (a) shows the clear-sky index time series for these averages, while panel (b) illustrates the averaging effect on 10 min increments in the clear-sky index. The misrepresentation of variability is evident for the longer averaging time in both panels. As indicated in Section 1, ground-based solar irradiance observations have often been averaged using a range of different temporal resolutions (from hours to fractions of a second), and there is no broad agreement on the ideal temporal resolution best-suited to record all relevant variability. Early studies showed temporal averaging on time scales larger than minutes to introduce considerable smoothing to the clear-sky index, and to affect its probability distribution [23,24]. More recent analyses have concluded that the temporal resolution required to determine irradiance variability across time scales may be as small as 0.1 s [17,63,65], 0.4 s [16], or 1 s [111]. The studies that argued for second-or-higher resolutions were based on 1. determining instantaneous irradiance variations for each of a few hundred days in spring and summer by calculating the second temporal derivative of each observation, considering the minimum (i.e., negative) value of a day's derivatives to represent the day's most severe fluctuation, and then computing an ideal averaging time by assuming the variations to feature parabolic shapes and accepting an error of 10 Wm −2 in the measurements [17]; 2. assessing the reduction of the standard deviation of an irradiance time series (measured during 7 h on a single summer's day) as a function of increasing averaging time scales [63]; 3. separately studying the variability index and the variability score of irradiance for seven selected days as functions of increasing averaging time scales [16]; and 4. characterizing the changes of k * and ∆k * τ standard deviations as a function of averaging time using thousands of hours worth of irradiance observations with raw temporal resolutions ranging from 0.01 through 1 s [111].
In addition to these differences in methods and quantity of data, some of the studies were limited to specific geographic areas (i.e., Southern Norway [17], Southern Finland [63], and Eastern Canada [16]), whereas another compared data from four different regions in the Northern Hemisphere (i.e., Eastern Canada, Germany, Hawaii, and Arizona) [111].
In practice, the effectively necessary temporal resolution of data strongly depends on the spatial scale considered. On the one hand, fluctuations of up to ±50% from one second to the next (and changes of more than 90% within 20 s) have been documented in the feed-in of a relatively small 48 kWp PV plant [56]. On the other hand, considerably larger multi-megawatt utility-scale PV plants may not require highly resolved measurements on the order of seconds for monitoring purposes [150], while minute-averaged observations may be resolved too coarsely [26].
In view of the relatively small number of high-resolution datasets available to date, several methods have been proposed to downscale more easily-accessible low-resolution data to smaller scales. For example, Markov-model-based approaches have been applied to 15-min PV power [118], hourly irradiance [112], hourly weather observations [69,109], or daily clearness index [151], in order to simulate variability on temporal scales well below those of the input data. Other downscaling approaches include providing a library of representative high-resolution irradiance samples [152] as well as simulating changes in the clear-sky index differently for different classes of variability [100,117]. A number of other sophisticated methods have been proposed to simulate relatively high-resolution data in time and/or space (using, for example, copulas or Markov chains), but without conditioning to large-scale information [27,41,61,105,153,154]. While these last-mentioned approaches were not specifically designed to downscale coarsely resolved real-world measurements, they can likely be adapted for this purpose.

Conclusions
This article concisely reviews recently published essentials from the literature regarding the quantification and small-scale averaging of irradiance variability in time and space. Complementing previous textbooks and literature reviews [15,[119][120][121][122][123][124][125][126][127][128][129][130][131][132][133], the paper emphasizes relatively small sub-minute scales below about 10 km. Despite the many articles touching on the subject of irradiance variability, suitable high-resolution measurement data are still relatively scarce. While some small-scale irradiance data have been publicly released [e.g., [72][73][74][75][76], there is still a need for more high-resolution measurements to robustly validate previous findings [52]. Granting open access to such data and corresponding models is considered best practice in order to foster scientific discussion and facilitate knowledge-based energy policies [155].
Beyond the collection of more suitable in-situ data, and with new satellite generations such as Meteosat Third Generation (MTG) in the pipeline, future efforts may also be geared towards inferring characteristics of small-scale variability from large-scale spaceborne observations. Two initial studies show great promise in this regard [67,99]. Both can be developed further in order to estimate complex climatological variability characteristics [67], and to eventually nowcast the character of short-term volatility in real time [99].
With regards to the applicability of irradiance variability research to PV power systems, there is a need for considering variability in irradiance on tilted surfaces as well as its effects on the processes of energy conversion. While most of the research presented in this review has been based on global horizontal irradiance, PV systems typically feature tilted modules. On a daily basis, irradiance variability has been shown to be higher on an inclined plane than on the horizontal [6,156], but corresponding analyses of high-resolution sub-minute data are still lacking. Likewise, an extensive validation of models separating direct and diffuse irradiance (a necessary step to predict irradiance on an inclined plane as a function of GHI) is not yet possible on sub-minute scales for lack of suitable data (it has, however, recently been performed using minute-scale data [157,158]). In general, small-scale specific phenomena such as cloud enhancement may call for future adaptations and extensions of well-established large-scale-based methods and conclusions [159].
The studies covered by this review have unanimously identified clouds as the dominant source of short-term fluctuations in surface irradiance and photovoltaic power. Beyond that, aerosols have been shown to considerably affect irradiance variability as well, albeit on relatively large spatiotemporal scales and with an emphasis on concentrated solar power applications [160][161][162][163][164][165][166][167][168][169][170][171]. In consideration of extreme dust events, such as forest fires or sand storms (which can quickly generate large numbers of aerosols [172]), aerosol-induced small-scale irradiance variability could be non-negligible and further research is necessary in this regard. However, corresponding high-quality data sets of both aerosol optical depth and surface irradiance are relatively difficult to obtain [160], especially on the small scales considered in this review.