Simulating clear-sky index increment correlations under mixed sky conditions using a fractal cloud model

Solar Energy 150 (2017) 255–264

Contents lists available at ScienceDirect

Solar Energy journal homepage: www.elsevier.com/locate/solener

Simulating clear-sky index increment correlations under mixed sky conditions using a fractal cloud model Gerald M. Lohmann a,⇑, Annette Hammer a,1, Adam H. Monahan b,1, Thomas Schmidt a, Detlev Heinemann a a b

Energy Meteorology Group, Institute of Physics, Oldenburg University, Germany School of Earth and Ocean Sciences, University of Victoria, Canada

a r t i c l e

i n f o

Article history: Received 11 January 2017 Received in revised form 13 April 2017 Accepted 18 April 2017

Keywords: Irradiance variability Fractal cloud edges Clear-sky index increments Spatial autocorrelation structures

a b s t r a c t Modelling clear-sky index increments (i.e. changes in normalized surface irradiance over specified intervals of time) and their spatial autocorrelation structures is important for the reliable grid integration of photovoltaic power systems. In order to capture increment correlation structures under mixed sky conditions, we apply a fractal cloud model. First, we use thousands of fish-eye sky camera and satellite images to estimate cloud edge fractal dimension D on scales between tens of meters and hundreds of kilometers. Box-counting analyses of cloud edges extracted from these images result in best-fit estimates of fractal dimension between 1.4 and 1.6, with satellite approximations being consistently higher than sky camera estimates. Contrary to previous studies, we find no evidence of any scale break and attribute the systematic discrepancy between the two estimates of D to intrinsic differences in the instrumentspecific techniques of cloud detection. Next, we synthesize small-scale fractal cloud fields, with D ¼ 1:5 and satellite-derived cloud images as input. Using cloud motion vectors from the same satellite images, we then translate the small-scale cloud fields across a set of model pyranometer locations to obtain corresponding clear-sky index time series. From the simulated time series we calculate increment correlation structures for distances between tens of meters and about ten kilometers, and compare them to observation-based results from 1 Hz measurements of up to 99 pyranometers. The observed isotropic, along-wind, and across-wind spatial autocorrelation structures of clear-sky index increments are captured well by the model, both in terms of overall value and shape. While there are some systematic differences between modelled and observed structures (e.g. underevaluation for very small time scales of a few seconds), the differences are essentially small. In general, the simulated correlations are not strongly sensitive to variations in the fractal model parameters. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction The world-wide installed capacity of photovoltaic (PV) power is expected to continue its recent fast-paced growth (e.g. an increase of  25% in 2015) (Solar Power Europe (SPE), 2016). PV power production is highly variable and increasing numbers of gridconnected PV systems will cause the associated problems (e.g. voltage and frequency stability) to consequentially increase as well (Perez et al., 2016b; Stetz et al., 2015). The relevant spatiotemporal scales of this variability range from days and hundreds of kilometers down to seconds and meters, and there is an ongoing need to characterize how PV power aggregation mitigates variability across these scales (Klima and Apt, 2015; Perez et al., 2016a; Widén et al., 2015). ⇑ Corresponding author. 1

E-mail address: [email protected] (G.M. Lohmann). These authors contributed equally to the paper.

http://dx.doi.org/10.1016/j.solener.2017.04.048 0038-092X/Ó 2017 Elsevier Ltd. All rights reserved.

In lieu of suitable PV power ramp data, clear-sky index increments (changes in normalized surface irradiance over specified intervals of time; also called ramp rates (e.g. Lave et al., 2012), step changes (e.g. Widén, 2015), or simply changes (e.g. Perez et al., 2012)) can be used to analyze PV variability characteristics of single points and spatial averages (Hoff and Perez, 2012). The presence of relatively high-magnitude increments is typical on all scales, with changes of up to tens of standard deviations away from the mean being recorded on a regular basis (Anvari et al., 2016). The probabilities of these large fluctuations decrease with increasing spatial and/or temporal scales (Lohmann et al., 2016). To assess variability reduction on different scales, it is useful to consider the spatial autocorrelation structures of clear-sky index increments (Lohmann et al., 2016). These structures can be used e.g. to predict power variability of large and/or distributed PV systems using a wavelet approach and a single irradiance sensor (Dyreson et al., 2014; Lave et al., 2013; Perpiñán and Lorenzo,

256

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264

2011; Perpiñán et al., 2013). Increment correlations (in both space and time) can depend on effective cloud speed (Hoff and Perez, 2010; Hoff and Perez, 2012), different daily variability categories (Perpiñán et al., 2013; Widén, 2015), or short-term sky types (Lohmann et al., 2016). In all cases, two-point correlation coefficients increase with increasing increment step size and/or decreasing distance (Hoff and Perez, 2012; Lohmann et al., 2016). Moreover, correlation coefficients can depend on the orientation of a sensor pair relative to the direction of cloud movement (Arias-Castro et al., 2014; Hinkelman, 2013). For along-wind correlations, the presence of relatively small negative peaks in the increment correlation structures has been suggested (Hinkelman, 2013; Lave and Kleissl, 2013; Perez et al., 2012; Widén, 2015). Several models have been proposed to reconstruct these spatial autocorrelation structures using empirical fits to measured data (Hoff and Perez, 2012; Lave and Kleissl, 2013; Lonij et al., 2013; Mills, 2010; Perpiñán et al., 2013; Widén, 2015) or simplified binary cloud shapes (Arias-Castro et al., 2014). Among these, a simple model based on large-scale satellite data (Hoff and Perez, 2012) can capture local small-scale correlations for overcast and clear skies, but fails to reproduce the structures for mixed sky conditions (Lohmann et al., 2016). These mixed sky increment correlation structures are primarily influenced by broken clouds (Hinkelman, 2013), and their representation can thus be improved through careful construction of synthetic binary cloud shadows and accounting for the anisotropic effects associated with their movement (Arias-Castro et al., 2014). Methods to generate simplified mixed-sky cloud shadow fields include the use of randomly positioned and subsequently blurred squares (Lave and Kleissl, 2013), unions of randomly positioned discs (Arias-Castro et al., 2014), and fractal models (Beyer et al., 1994; Cai and Aliprantis, 2013). The fractal approach has been applied to simulate irradiance and PV power time series for single points and spatial averages, without considering model-derived spatial correlations of increments (Beyer et al., 1994; Cai and Aliprantis, 2013). As for appropriate cloud edge fractal dimensions, previous studies have suggested the presence of a sub-kilometer scale break in box-counting dimension (Beyer et al., 1994) and cloud area-perimeter relation (Cahalan and Joseph, 1989; Gotoh and Fujii, 1998). However, these studies were based on very limited datasets (11 sky photographs (Beyer et al., 1994) and single high-resolution Landsat images (Cahalan and Joseph, 1989; Gotoh and Fujii, 1998)). Moreover, other studies have presented evidence for the absence of such a scale break, arguing for unified scaling in atmospheric variability (e.g. Lovejoy et al., 1993; Lovejoy and Schertzer, 2006; Sachs et al., 2002). The data used to assess the representation of increment spatial autocorrelation structures in the above-mentioned models were relatively coarsely resolved in space and time. For example, the most detailed set of irradiance measurements used in such analyses originates from a collection of 17 pyranometers deployed over an area of about 1 km2 (Arias-Castro et al., 2014; Hinkelman, 2013). More recent radiation data sets (Gagné et al., 2016; Lohmann et al., 2016) allow the analysis of irradiance variability on fine temporal and spatial scales of seconds and tens of meters, respectively. In this paper, we use a very extensive set of data from two dense pyranometer networks, a co-located fish-eye sky imager, and the Meteosat-10 satellite (described in detail in Section 2) in order to 1. revisit estimates of the fractal dimensionality of cloud edges under mixed sky conditions on scales between tens of meters and hundreds of kilometers (Section 3), 2. synthesize small-scale fractal cloud shadows, with satellitederived cloud cover fraction as input (Section 4), and

3. assess the suitability of these artificial shapes for modelling spatial autocorrelation structures of clear-sky index increments, including along-wind and cross-wind structures (Section 5). A sensitivity analysis (Section 6) and conclusion (Section 7) finally complete the work. 2. Data All data sets used in this study have been presented in previous publications, which discuss in detail how the data were collected and processed (Hammer et al., 2003; Kühnert et al., 2013; Lohmann et al., 2016; Madhavan et al., 2016; Schmidt et al., 2016). Here, we will concisely recapitulate the aspects most relevant to the present study of 

1. clear-sky index k field data calculated from irradiance measurements of up to 99 pyranometers and a clear-sky model, 2. red-blue-ratio r maps for cloud detection derived from photographs of a fish-eye sky imager and cloud base height (CBH) estimations from a ceilometer, as well as 3. cloud index n maps and cloud motion vectors v based on data from the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) on board the Meteosat Second Generation (MSG) spacecraft. Fig. 1 provides illustrative examples of the data, and Table 1 summarizes the key variables used throughout the text. 2.1. Clear-sky index from pyranometer networks For the High Definition Clouds and Precipitation for Climate Prediction (HD(CP)2) Observational Prototype Experiment (HOPE), two pyranometer networks were deployed to measure global horizontal irradiance G with high spatio-temporal resolution. During the first campaign near Jülich, Germany (50.9°N, 6.4°E), 99 photodiode pyranometers were distributed over an area of about 80 km2 between April 2 and July 24, 2013 (see Fig. 1(a)). A subset of 50 of these pyranometers was later set up in a denser array over an area of about 4 km2 near Melpitz, Germany (51.5°N, 12.9°E) from September 3 until October 14, 2013 (see Fig. 1(b)). The original sample temporal resolution of 10 Hz was averaged to 1 Hz during post-processing, and simultaneous clear-sky irradiance Gclear was estimated using the model described by Fontoynont et al. (1998). For all solar elevation angles a > 15 , clear-sky index time series 

k ¼

G ; Gclear

ð1Þ

were then calculated for each sensor (see Fig. 1(e) for an example of  900 s worth of k field statistics). Madhavan et al. (2016) provided more information about the equipment, measuring inaccuracies, and quality flags of the irradiance measurements, and Lohmann et al. (2016) discussed in detail the processing to obtain the very  k time series used in this study. 2.2. Sky imager photographs and ceilometer During the daylight hours of the HOPE Jülich campaign, a sky imager (a digital CCD camera with fish-eye lens (Kalisch and Macke, 2008)) took hemispheric photographs every 15 s, realizing a resolution of 2592  1744 pixels and a 183° field of view. For the purpose of cloud detection, the color triplets (red; green; blue) of the original image were converted to values of red-blue-ratio 1

r ¼ red  blue , where low values indicate clear pixels and high values indicate cloudy ones.

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264

257

Fig. 1. Examples of available data sets: (a,b) show the relative positions of pyranometers during the field campaigns in (a) Jülich and (b) Melpitz, (c) shows a sky-imagerbased red-blue-ratio map, (d) shows a simultaneous SEVIRI cloud index image including individual (red arrows) and main (black arrow) cloud motion vectors, and (e) presents field statistics of clear-sky index, measured in Jülich during the 15 min following the two images. These statistics consist of the time-evolving probability  distribution of k across all pyranometers (greyscale), the ensemble mean (blue curve) and a single realization corresponding to a randomly-selected pyranometer (red curve). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1 Key variables used throughout the text. Name Clear-sky index Red-blue-ratio Cloud index Cloud motion vector

Symbol 

k r n

v

Description Irradiance normalized to clear-sky conditions (Lohmann et al., 2016) Sky imager product used for cloud identification in images (Schmidt et al., 2016) Satellite product representing cloud amount (Hammer et al., 2003) Displacement of similar cloud patterns found in two subsequent cloud index images (Kühnert et al., 2013)

Estimates of cloud base height were additionally available from a ceilometer located in the immediate vicinity of the sky imager. The ceilometer’s original single-point measurements were recorded every 20 s and then smoothed using a median filter (Astola and Kuosmanen, 1997) with a 600 s window. Using these smoothed CBH estimates, the red-blue-ratio images were undistorted and projected onto a Cartesian coordinate system, on which the values were resampled on a regular grid. The images we consider are the same as those used by Schmidt et al. (2016) (who also detailed all of the necessary steps to convert the fish-eye sky images to geolocated maps of r values) with the exception of the following changes: 1. generating a higher-resolution final image of 500  500 pixels with an edge length of 10 m per pixel, 2. masking all areas where the required horizontal pixel size of 100 m2 is exceeded in the original image (either due to distortion by the fish-eye lens or high CBH values), and 3. masking a 25° angular distance around the position of the sun, if the mean pixel intensity (i.e. gray scale value in the range between 0 and 255) in this area indicates an unobstructed view of the sun. This determination is made if the gray scale exceeds a specified threshold (taken to be 180).

2.3. MSG-based cloud index and motion vectors Based on quarter-hourly SEVIRI measurements on board the MSG satellite, cloud index estimates and cloud motion vectors have been routinely processed and archived since 2004. The cloud index relates measured reflectance to approximated surface albedo values and takes values between 0:2 (completely clear) and 1:2 (completely cloudy). The data processing used to calculate cloud index images is presented in detail by Hammer et al. (2003). We use both their original satellite product (with an irregular spatial resolution of about 1:2  2 km at the HOPE campaigns’ locations), as well as resampled images with a regular grid of 500  500 pixels with an edge length of 2.5 km per pixel. Cloud motion vectors v were estimated by matching similar cloud patterns found in two subsequent cloud index images, as explained by Kühnert et al. (2013). Single motion vectors can be unreliable, so we derive a main cloud motion vector v m from the 9 vectors closest to the centroid of the sensor locations’ distribution by calculating the respective medians of the vectors’ East and North components. Fig. 1(d) presents an example of a cloud index image section and corresponding individual cloud motion vectors (red arrows) along with the derived main vector (black arrow). 2.4. Mixed sky subset selection

These conservative modifications of the processing procedure ensure high-resolution red-blue-ratio maps, and minimize cloud shape deformation as well as misclassification of bright but clear pixels as cloudy around the position of the sun. Fig. 1(c) shows an example of a red-blue-ratio map including a masked region at the bottom of the image.

From the large dataset available, we pick only those data that are classified as mixed sky conditions (as defined in Lohmann et al., 2016). For each 900 s window following a SEVIRI image, the standard deviation of all pyranometer-based clear-sky index field samples is calculated. If the standard deviation exceeds a cer-

258

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264

tain threshold (taken to be 0.18), the respective period and its preceding SEVIRI image are regarded as mixed sky data. This results in a total of 1688 such quarter-hour periods (1463 originating from the Jülich campaign and 225 from Melpitz), with as many SEVIRI images and roughly 1:5  106 s worth of corresponding  pyranometer-based k field measurements forming the data basis  of all analyses. These mixed sky k values are bimodally dis tributed, with a local minimum around the median P 50 ðk Þ ’ 0:8 separating the two peaks within which the first quantile   P25 ðk Þ ’ 0:45 and the third quantile P 75 ðk Þ ’ 1 are located (not shown). These peaks correspond to instantaneous cloud-covered and cloud-free states, respectively. Most absolute values of the associated cloud motion vectors are relatively small, with the first quantile P 25 ðjjv m jjÞ ¼ 2:4 m s1 , the median P50 ðjjv m jjÞ ¼ 5:4 m s1 , and the third quantile P75 ðjjv m jjÞ ¼ 11:0 m s1 . For the estimation of cloud edge fractality in Section 3, an additional set of sky-imager-based red-blue-ratio maps are available every 15 s for two months of the Jülich campaign (April and May 2013). Among the thousands of sky images taken during mixed sky windows, only 199 were taken simultaneously (i.e. within 10 s) with a satellite image. In order to strengthen the red-blueratio data basis without processing an unreasonably large number of images (two consecutive sky images do not usually differ appreciably), we consider an extra 1000 randomly selected r-field snapshots originating from time windows classified as mixed sky. 3. Cloud edge fractal dimensions From each map of cloud index and red-blue-ratio, we extract the shapes of cloud edges by 1. generating a threshold-based binary cloud image (using respective empirical thresholds nth ¼ 0:3 and r th ¼ 0:85), and 2. defining a cloud-covered pixel in this binary image as making up part of a cloud edge if it is immediately surrounded by at least 1 cloud-free pixel. The thresholds were obtained by analysing histograms of n and r during days of variable cloud cover, such that nth and r th effectively separate the two populations respectively corresponding to cloud-covered and cloud-free values (see e.g. Schmidt et al., 2016; Skartveit and Olseth, 1992). To characterize their fractal properties, we subject all resulting 500  500 pixel cloud edge images to a box counting analysis (Barnsley, 2014) using 45 different box sizes (with edge lengths l decreasing by 5 pixels from lmax ¼ 225 pixels down to lmin ¼ 5 pix-

els). For each scaled box size  ¼ l  L1 (with L ¼ 500 pixels denoting the edge length of an entire image), we use 5 image rotations (by steps of 72°) around each of 5 randomly positioned box origins to obtain 25 samples of each image. Based on initial analyses of a few representative images, we found the accuracy of the fractal dimension estimates to not appreciably increase for more rotations and/or box origins (not shown). We then count the numbers of non-overlapping boxes necessary to cover the rotated cloud edges separately for each of these samples and subsequently derive the average box counts N  for each . Considering the relationship between  and N  we can then estimate the cloud edge boxcounting dimension D as the slope m of a least-squares linear regression between lnðN  Þ and lnð1 Þ. For each data source, we use the ImageJ (Abràmoff et al., 2004; Rasband, 2016; Schneider et al., 2012) plugin FracLac (Karperien, 2015) to both determine the box-counting dimensions separately for each image and globally for all images. Fig. 2 summarizes the results of this analysis. Panels (a,b) show examples of cloud edges derived from the images of cloud index

and red-blue-ratio previously shown in Fig. 1 (the shaded area in panel (a) indicates the image section shown in Fig. 1(d)). Panels (c,d) present the linear regressions of lnðN  Þ and lnð1 Þ using all cloud edge images from the satellite and the sky imager, respectively, also quoting the corresponding fractal dimension estimates. Panel (e) presents histograms of the single-image dimension estimates from sky-imager- and satellite-derived values, and panel (f) displays a scatter plot of simultaneously estimated dimension values from the two data sources. The box-counting dimensions estimated from satellite data are consistently higher than those based on the sky imager, both for the single-image analyses (satellite: 25th percentile P25 ¼ 1:55 and 75th percentile P 75 ¼ 1:65; sky imager: P25 ¼ 1:30 and P75 ¼ 1:46), and for the average over all images (satellite: m ¼ 1:6; sky imager: m ¼ 1:4). The single-instrument regression lines fit the data points very well on all scales. We attribute the small slope differences between instruments to intrinsic differences in the cloud images obtained by the sky imager and SEVIRI (e.g. the sky imager can more easily happen to see the sides of clouds), and the different techniques of cloud detection (cf. Section 2). Likewise, the variability in single-image estimates is likely to be caused by sampling, rather than there being systematic differences in fractal dimension for different cloud scenarios, because the simultaneously estimated slopes from the two data sources appear uncorrelated (Fig. 2(f)). In contrast to previous suggestions of a sub-kilometer fractal scale break (Beyer et al., 1994; Cahalan and Joseph, 1989; Gotoh and Fujii, 1998), our results show no evidence of such a break, and are in line with findings of uniform scaling (e.g. Lovejoy et al., 1993; Lovejoy and Schertzer, 2006; Sachs et al., 2002). While we cannot exclude the possibility of a scale break in the small range of length scales between those covered by the two instruments (roughly between 2 and 10 km), the subkilometer scales are well-represented by the sky imager and show no evidence of a scaling break. 4. Cloud shadows and time series generation Synthetic clear-sky index time series can be generated by producing fractal cloud images and then moving them in the direction of the mean wind. Our approach is based on two previous publications (Beyer et al., 1994; Cai and Aliprantis, 2013) that used the diamond-square algorithm (Fournier et al., 1982) to synthesize cloud shapes with specified cloud edge fractal dimension. We modify the method, in particular by using SEVIRI cloud index data to specify values of cloud cover and cloud motion vectors, and by using a Hurst exponent (Mandelbrot and Wallis, 1968) H ¼ 0:5 on all scales. This value of H corresponds to a fractal dimension D ¼ 1:5 (i.e. H ¼ 2  D), which is the average of the box-counting dimensions estimated from the sky imager (1.4) and satellite (1.6) in Section 3. The fractal construction (denoting fractal values as x) is carried out on square boxes, so the relevant cloud domain (determined by the area under consideration and the motion of clouds) must be represented as a set of non-overlapping squares. For each SEVIRI image, we generate such a synthetic cloud strip and model 900 s worth of clear-sky index time series according to the following procedure (see Fig. 3 for corresponding illustrations): (A) Observations and data preparation 1. For a given SEVIRI image, the coordinates of specified locations are defined to represent a set of virtual ground-based irradiance sensors. Here, we use the pyranometers’ coordinates from the HOPE campaigns. 2. The image and sensor locations are transformed into along-wind and cross-wind coordinates (relative to the main cloud motion vector v m ) by rotation around the sensor locations’ centroid.

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264

259

Fig. 2. Results from cloud edge box-counting analyses: (a,b) cloud edges based on the images in Fig. 1 (the shaded area in (a) corresponding to the domain shown in Fig. 1(d)), (c,d) log-log plots of number of boxes vs inverse of normalized box size along with best-fit linear regression lines, (e) histograms of single-image regression slopes, (f) scatter plot of simultaneously estimated slopes using the two instruments (199 points).

3. The number of squares necessary so that all virtual sensor locations will be covered by an along-wind strip of squares moving with the velocity v m for 900 s is calculated using a cross-wind square edge length of 513 pixels with a resolution of 6 m (Melpitz) or 25 m (Jülich) each (see Fig. 3(a)). 4. Within this along-wind strip, the mean cloud index n is calculated. (B) Fractal and time series generation 1. An empty model domain of 513 cross-wind pixels and the appropriate number of along-wind pixels (such that the domain has the same dimensions as the strip used for the observations) is initialized. 2. Initial corner values x0 of all big squares in the domain are randomly generated (at the positions denoted ‘‘0” in Fig. 3 (b)), based on a Gaussian distribution with mean 0 and standard deviation r0 ¼ 1. 3. Fractal values x of the remaining pixels are constructed stage by stage, with each stage consisting of two steps, namely (a) square center value generation, and (b) diamond center value generation. According to the diamond-square algorithm, each value is compiled by adding the arithmetic mean of the respective corner values to a realization of a Gaussian random variable with mean 0 and standard deviation

(

pffiffiffi

r0  2Hi  2; for square centers ri ¼ r0  2Hi ; for diamond centers;

ð2Þ

depending on the stage number i. In Fig. 3(b), an example of the first stage i ¼ 1 is shown with steps ‘‘a” (red arrows)

and ‘‘b” (blue arrows), along with the first step of the second stage ‘‘2a”. 4. Based on the empirical cumulative distribution function (CDF) obtained from all x values, a threshold value xth is obtained, such that the probability that x < xth is numerically equal to the mean satellite-derived cloud index n, linearly mapped to the range between 0 and 1 (see Fig. 3(c)). 5. Around the value xth , a transition zone is defined using the parameter p ¼ 0:15, such that all fractal values (a) x < xth  p are set to 0.2, (b) x > xth þ p are set to 1.2, and (c) xth  p < x < xth þ p are linearly mapped to the range between 0.2 and 1.2. The resulting artificial high-resolution cloud shapes feature smooth transitions between clouds and sky, as well as fractal edges with box-counting dimensions of  1:5. By construction, averaging over the entire domain results in the same mean cloud index as obtained from the original SEVIRI image (although the cloud realizations are different). The parameter p is a new element of the model, i.e. previous applications of this algorithm took p ¼ 0 (though Cai and Aliprantis (2013) employed a so-called ‘‘multi-layer rendering technique” to blur the cloud edges). This new parameter is introduced to allow a continuous transition between cloudy and cloud-free pixels in the modelled cloud field. Section 6 discusses the sensitivity of the model results to the value of p. 6. The modelled cloud index values in the full domain are then converted to clear-sky index using the empirical relationship (Rusen et al., 2013)

260

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264



Fig. 3. An illustrative example of the algorithm used to generate artificial high-resolution cloud shadows and clear-sky index k time series, using the same SEVIRI image and 900 s period as in Fig. 1: (a) rotated SEVIRI image of cloud index n with the main cloud motion vector v m (red), Jülich sensor coordinates (black dots), and a sufficient set of along-wind boxes (black dashed lines); (b) fractal field generated with the diamond-square algorithm, overlaid with the first steps of the model sequence from random initialization (0), via the first stage of square center (1a, red arrows) and diamond center (1b, blue arrows) generation, to the beginning of the second stage (square center, 2a); (c) cumulative density function (CDF) of fractal field values with indication of the threshold and transition range; (d) final high resolution cloud index field; (e) enlarged area  to emphasize the smoothness of the cloud edge transitions; (f) field statistics of modelled k (as in Fig. 1(e)). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

8 1:2; > > < 1  n;

n 6 0:2 0:2 < n 6 0:8  k ¼ 2 > 1:1661  1:7814  n þ 0:7250  n ; 0:8 < n 6 1:05 > : 0:09; n > 1:05: ð3Þ

While other models to convert cloud index to clear-sky index have been developed (e.g. Perez et al., 2002), Eq. (3) is routinely employed in satellite-based irradiance retrieval using the same SEVIRI data considered in this study. Due to SEVIRI’s quarter-hourly temporal resolution and its spatial footprint on the order of a few square kilometers, Eq. (3) does not account for temporary cloud enhancement (Piacentini et al., 2011; Yordanov et al., 2013) and other local short-term phenomena (it compares best to hourly-averaged ground-based pyranometer measurements when used with a mean cloud index derived from 3  5 pixels and 4 – 5 consecutive SEVIRI images). Although the absolute accuracy of Eq. (3) may consequently be expected to decline in our high-resolution application, this is of secondary importance when assessing the spatial structures of increment correlations. In  Fig. 3(d) the resulting k field is presented; a magnified section is shown in Fig. 3(e). 7. Clear-sky index time series are obtained by virtually mov ing the k shadow maps over the initially defined sensor  locations (where k can be sampled) with the satellitederived cloud speed (see Fig. 3(f)).

When comparing the variability characteristics of the modelled clear-sky index time series in Fig. 3(f) to those of the measured time series previously shown in Fig. 1(e), there are both evident differences and similarities. As for the former, the abso lute ranges of modelled and measured k values are not identical (i.e. vertical axes are scaled slightly different), because the relation between cloud index and clear-sky index given by Eq. (3) is based on large-scale satellite observations and does not perfectly represent local short-term conditions. The range of  all modelled values is confined to exactly 0:09 < k < 1:2 by design, while the observed ranges are influenced by local short-term phenomena such as cloud enhancement and differ slightly between sensors (e.g. due to errors of measurement, such as marginally tilted or imperfectly calibrated pyranometers), leading to a less definite overall range in the measured  k field statistics. As for agreement between simulated and  observed k variability, the randomly selected modelled sensor is dominated by a similar on/off behaviour as its measured counterpart, and even the occasional short-duration changes of intermediate magnitudes are reproduced well. The modelled field statistics also resemble the measured ones, with both fea turing times of wide and narrow k distributions. Finally, the variable character of the measured and modelled sensor ensemble mean compares very well: both values drift slowly (relative  to the 900 s window) in the range 0:4 < k < 1:0 with relatively little short-term noise.

261

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264

5. Increment spatial correlations 

From the observed and synthetic k time series, we derive clearsky index increments

Dks ðtÞ ¼ k ðt þ sÞ  k ðtÞ; 





ð4Þ

for different time lags s (measured in seconds). Based on these increment time series we subsequently calculate the spatial twopoint Pearson correlation coefficients between locations i and j Dks ij

q

PT









t¼1 ðDks;i ðtÞ  Dks;i ÞðDks;j ðtÞ  Dks;j Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; PT   2 PT   2 t¼1 ðDks;i ðtÞ  Dks;i Þ t¼1 ðDks;j ðtÞ  Dks;j Þ

ð5Þ

where T ¼ 900  s denotes the number of Dks samples in time. By separately selecting only those sensor pairs within a range of 22:5° around the direction of the cloud motion vector and the perpendicular range of directions, we additionally assess along-wind 

Dk

and cross-wind effects on qij s . Fig. 4 compares modelled spatial autocorrelation structures of  Dks with results derived from measurements for three distinct increment time steps of s ¼ 1 s, s ¼ 10 s, and s ¼ 100 s, using 10 logarithmically spaced bins (3 leftmost columns). For each increment value, the distribution of simulated and measured crosscorrelations for all pairs of stations (combined from both Melpitz and Jülich), for each 900 s window, are shown along with the median value at each distance. The figure also includes the root-meansquare difference (RMSD) and the signed mean difference (SMD) (i.e. bias) (e.g. Lorenz et al., 2009) of these bins’ medians as functions of increment time step s to assess model accuracy (2 rightmost columns). The top row shows results of an isotropic analysis (in which correlation structures are computed as functions of distance, over all directions), while the two bottom rows differentiate between along-wind and cross-wind correlations as described above. The spread around the medians of the observation-derived results is consistently larger than that of the model runs in panels (a–i). This broader range could be caused by different cloud types (e.g. stratiform vs. cumuliform; different cloud base and thickness) all combined together in the mixed cloud category, but with different effects on irradiance. While the fractal cloud model does not

include such effects, it reproduces the median spatial autocorrelation structures well. The apparently large deviations between modelled and observed correlation values for s ¼ 1 s are differences between very small correlation values, and in both cases indicate negligible spatial structure for this short increment step. Disregarding the smallest time lag, the model captures the magnitude and shape of the correlation structures particularly well – the upwards curvature for s ¼ 10 s, as well as the inflection structure for s ¼ 100 s (concave down for small distances, then concave up). This is a rather subtle feature of the correlation structures, and it is non-trivial that the model captures it. Judging by the RMSD of the isotropic analysis (top row), model performance fluctuates slightly for increasing s, with best performance (i.e. low RMSD) around s ¼ 4 s and s ¼ 80 s, and worst performance (i.e. high RMSD) around s ¼ 1 s; s ¼ 20 s, and s ¼ 500 s. The RMSD of the along-wind and cross-wind analyses exhibit similar characteristics. The overall range of RMSD values is small in all cases (always < 0:08), suggesting a good model accuracy. The SMD plots demonstrate that the model underestimates correlations for s < 10 s and overestimates them for s > 10 s in the isotropic and cross-wind analyses, while along-wind correlations are underestimated for all timescales, except between about s ¼ 5 s and s ¼ 10 s. Here also, the range of values is very narrow (0:04 < SMD < 0:03), indicating a good overall model accuracy. Some previous studies have reported small negative peaks in Dk

along-wind correlation structures (qij s J  0:4 on the short time scales investigated in our study). These studies were based on analyses of single-sensor-derived virtual networks (Perez et al., 2012; Widén, 2015), a small multi-sensor pyranometer network (Hinkelman, 2013), and a simple cloud simulator (Lave and Kleissl, 2013). While some of the simulated and observed sensor pairs in Fig. 4 exhibit negative along-wind correlation coefficients, the respective medians do not become appreciably negative. This difference from previous findings may be because (1) the virtual networks assume a perfect along-wind virtual sensor configuration, while we allow for a range of 22:5° around the direction of the cloud motion vector, (2) the characteristics of increment correlation structures may differ considerably between different locations (Hinkelman (2013) used data from Hawaii, where trade wind cumulus cloud fields are typical throughout the year; while HOPE was conducted in Germany, where cloud structures are more



Fig. 4. Measured and modelled spatial autocorrelation structures of clear-sky index increments Dks as functions of sensor pair distance (combined from both Melpitz and Jülich) for increment time steps (a–c) s ¼ 1 s, (d–f) s ¼ 10 s and (g–i) s ¼ 100 s. The medians of 10 logarithmically spaced bins are shown along with shaded indications of the respective underlying distributions. The corresponding (j–l) root mean square differences (RMSD) and (m–o) signed mean differences (SMD) of medians (computed by averaging across all spatial scales) are presented as functions of time scale. All results are conditioned on relative wind directions (top row: isotropic, middle row: along-wind, bottom row: cross-wind). The vertical axes’ scales are set to only differ between columns.

262

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264

Fig. 5. Model sensitivity to changes of the parameter p (representing cloud-edge smoothness, top row), the cloud edge fractal dimension D (middle row), and the main cloud motion vector v m (bottom row) corresponding to the isotropic analysis in Fig. 4. Each parameter is shown using its default (red circle), an increased (green upward-triangle),  and a decreased (blue downward-triangle) value. The panels show measured and modelled median spatial autocorrelation structures of clear-sky index increments Dks as functions of sensor pair distance (combined from both Melpitz and Jülich) for increment time steps (a–c) s ¼ 1 s, (d–f) s ¼ 10 s and (g–i) s ¼ 100 s, as well as the corresponding (j–l) root mean square differences (RMSD) and (m–o) signed mean differences (SMD) of medians as functions of time scale. The model sensitivities for alongwind and cross-wind analyses are qualitatively similar (not shown). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

diverse), and (3) the fractal cloud edges are inherently scaleinvariant and do not feature any predominant length scales of cloud size or cloud separation distance, while the more regular shapes used by Lave and Kleissl (2013) include fixed length scales leading to negative correlation coefficients by design. 6. Sensitivity analysis We assess the sensitivity of the modelled spatial autocorrelation structures of clear-sky index increments to the predefined 1. parameter p, representing the bluriness of the synthesized cloud edges, 2. fractal dimension D used by the diamond-square algorithm, and 3. main cloud motion vector v m , based on which the localized time series are generated. Sensitivity to these particular parameters is considered because p is without a priori constraints, and estimates of D and v m are uncertain. The parameter variations are chosen within the ranges of 0 < p < 0:3; 1:4 < D < 1:6, and 0:5 v m < v m < 1:5 v m , with the default values respectively amounting to p ¼ 0:15; D ¼ 1:5, and directly-estimated v m . Fig. 5 accordingly presents modelled structures for time steps of s ¼ 1 s, s ¼ 10 s, and s ¼ 100 s in panels (a–i), based on the default parameter values as well as the minimum and maximum of these ranges. The observation-based median is also provided in these panels, and the RMSD and SMD of the medians are presented in panels (j–o) to assess accuracy as a function of increment time. Results are shown for the isotropic analyses only, as along-wind and cross-wind sensitivities are similar. A systematic increase of correlation coefficients in panels (a–i) can be observed for 1. increasing values of p, which widen the transitional zone between ‘‘cloud” and ‘‘no-cloud” in the model, thereby intro ducing more intermediate values of k and smoother cloud edge transitions in the overall cloud images, 2. decreasing fractal dimensions that result in less complex fractal shapes and thus more regular overall cloud images, and 3. increasing cloud speeds, because the covariability is then dom inated by larger structures in the simulated k field. In particu-

lar, two locations are more likely to be affected by the passing of the same cloud edge within a certain time for larger speeds. The respective opposites apply to decreasing correlation coefficients. For the shortest time lag s ¼ 1 s (panels (a–c)), the default model settings do not provide the best results. Instead, an increased p (corresponding to smoother cloud edge transitions), a decreased D (as obtained from the sky image box-counting analysis), or an increased v m match the observational results slightly better. However, for s ¼ 10 s (panels (d–f)) and s ¼ 100 s (panels (g–i)), the default parameters yield the best results, while the correlation values from the varied parameters consistently overestimate or underestimate the observations. Based on the RMSD (panels (j–l)) and SMD (panels (m–o)), a cross-over time scale can be identified in the order of a few seconds, starting at which the default model settings generally perform better than the other values considered (i.e. lowest RMSD and SMD closest to zero). In general, the simulated correlation structures are not strongly sensitive to variations in these parameters. 7. Conclusions In the first part of this study, we revisited the issue of cloud edge fractal dimension under mixed sky conditions using thousands of fish-eye sky camera and satellite images. Box-counting analyses of cloud edges extracted from these images resulted in best-fit estimates of fractal dimension between 1.4 and 1.6 for these conditions, with satellite approximations being consistently higher than sky imager estimates. However, we found no evidence of any scale break and attribute the systematic discrepancy between the two estimates of D to intrinsic differences in the instrument-specific techniques of cloud detection. Using D ¼ 1:5 as a compromise estimate, we then synthesized corresponding fractal cloud fields and translated these across a set of model pyranometer locations. From the resulting clear-sky index time series we calculated spatial autocorrelation structures of clear-sky index increments, and compared these to results originating from two extensive pyranometer measurement campaigns. The observed isotropic, along-wind, and across-wind correlation structures are captured well by the model, both in terms of overall value and

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264

shape. While there are systematic differences between modelled and observed correlation structures (e.g. underevaluation for very small s; cf. Fig. 4), the differences are generally small. Neither the previously suggested scale break in cloud edge fractal dimensions (Beyer et al., 1994; Cahalan and Joseph, 1989; Gotoh and Fujii, 1998), nor the earlier anticipated slight negative peak in along-wind correlations (Hinkelman, 2013; Lave and Kleissl, 2013; Perez et al., 2012; Widén, 2015) are confirmed by our findings. Both proposals were based on very limited data sets compared to those used in the present study. These differences emphasize the ongoing need to collect extensive high-resolution multi-sensor irradiance measurements in order to correctly characterize irradiance field variability for use in model development and assessment. In particular, the extensive experimental datasets that have recently become available can be used for a systematic evaluation of other recently published irradiance variability correlation studies and models (Arias-Castro et al., 2014; Hinkelman, 2013; Hoff and Perez, 2010; Hoff and Perez, 2012; Lave and Kleissl, 2013; Lonij et al., 2013; Mills, 2010; Perpiñán et al., 2013; Widén, 2015). The present study contributes to the greater goal of synthesizing irradiance fields with a correct representation of increment correlation structures on a range of spatiotemporal scales, including very small ones. However, the fractal model developed using midlatitude observations may not necessarily be applicable to different climatic regions, where specific cloud types prevail (e.g. trade wind cumulus in Hawaii). Without considering wind, synthetic but realistic irradiance fields are already useful to study instantaneous solar irradiance; when coupled with wind vector information, they can produce realistic simulations of photovoltaic power generation time series and probability distributions of discretionary numbers of arbitrarily distributed PV systems. Ultimately, this allows for detailed studies of how the structures of distributed PV generation networks affect power quality and electrical grid stability. With increasing numbers of grid-connected PV power systems, this is essential to a successful grid operation in the future.

Acknowledgements We thank Andreas Macke at the Leibniz Institute for Tropospheric Research TROPOS (Leipzig, Germany) for sharing the pyranometer network datasets of the HD(CP)2 Observational Prototype Experiment (HOPE), and Jan Kühnert for providing cloud motion vector data for Germany. Also, we acknowledge helpful comments from two anonymous referees. The time and effort invested in developing and maintaining R (version 3.2.5) by the R Core Team (2016) and the active community of package authors is gratefully appreciated. The same holds true for the developers and maintainers of ImageJ (Abràmoff et al., 2004; Rasband, 2016; Schneider et al., 2012) and its plugin FracLac (Karperien, 2015). The model simulations were performed at the HPC Cluster HERO, located at the University of Oldenburg (Germany) and funded by the DFG through its Major Research Instrumentation Programme (INST 184/108-1 FUGG) and the Lower Saxony Ministry of Science and Culture. This research was partially funded by the Lower Saxony research network ’Smart Nord’, which acknowledges the support of the Lower Saxony Ministry of Science and Culture through the ’Niedersächsisches Vorab’ grant program (grant ZN 2764/ ZN 2896). It was also partly funded by the ’Performance Plus’ research project through the European Union’s Seventh Framework Program for research, technological development and demonstration (grant agreement No. 308991), and received support from the German Federal Ministry for Economic Affairs and Energy BMWi (Bundesministerium für Wirtschaft und Energie – grant 03ET4027B).

263

We also acknowledge funding support from the Natural Sciences and Engineering Research Council of Canada. References Abràmoff, M.D., Magalhães, P.J., Ram, S.J., 2004. Image processing with ImageJ. Biophoton. Int. 11 (7), 36–42. Anvari, M., Lohmann, G., Wächter, M., Milan, P., Lorenz, E., Heinemann, D., Tabar, M. R.R., Peinke, J., 2016. Short term fluctuations of wind and solar power systems. New J. Phys. 18 (6), 063027 http://stacks.iop.org/1367-2630/18/i=6/a=063027? key=crossref.69a7aa46f7b7166259cd8984f57b01a9. Arias-Castro, E., Kleissl, J., Lave, M., 2014. A Poisson model for anisotropic solar ramp rate correlations. Sol. Energy 101 (Mar.), 192–202 http://linkinghub. elsevier.com/retrieve/pii/S0038092X13005549. Astola, J., Kuosmanen, P., 1997. Fundamentals of Nonlinear Digital Filtering, vol. 8. CRC Press. Barnsley, M.F., 2014. Fractals Everywhere. Academic Press. Beyer, H.G., Hammer, A., Luther, J., Poplawska, J., Stolzenburg, K., Wieting, P., 1994. Analysis and synthesis of cloud pattern for radiation field studies. Sol. Energy 52 (5), 379–390 http://www.sciencedirect.com/science/article/pii/0038092X949 0115I. Cahalan, R.F., Joseph, J.H., 1989. Fractal statistics of cloud fields. Monthly Weather Rev. 117 (2), 261–272 http://dx.doi.org/10.1175/1520-0493(1989)117<0261: FSOCF>2.0.CO;2. Cai, C., Aliprantis, D., 2013. Cumulus cloud shadow model for analysis of power systems with photovoltaics. IEEE Trans. Power Syst. 28 (4), 4496–4506. Dyreson, A.R., Morgan, E.R., Monger, S.H., Acker, T.L., 2014. Modeling solar irradiance smoothing for large PV power plants using a 45-sensor network and the wavelet variability model. Sol. Energy 110 (Dec.), 482–495 http:// linkinghub.elsevier.com/retrieve/pii/S0038092X14004678. Fontoynont, M., Dumortier, D., Heinemann, D., Hammer, A., Olseth, J., Skarveit, A., Ineichen, P., Reise, C., Page, John, H., Roche, L., Beyer, H.-G., Wald, L., 1998. SatelLight: A WWW server which provides high quality daylight and solar radiation data for Western and Central Europe. 9th Conference on Satellite Meteorology and Oceanography, vol. EUM-P-22. American Meteorological Society Ed., Boston, Massachusetts, USA, Paris, France, pp. 434–437 https:// hal-mines-paristech.archives-ouvertes.fr/hal-00472444. Fournier, A., Fussell, D., Carpenter, L., 1982. Computer rendering of stochastic models. Commun. ACM 25 (6), 371–384 http://dl.acm.org/citation.cfm?id= 358553. Gagné, A., Turcotte, D., Goswamy, N., Poissant, Y., 2016. High resolution characterisation of solar variability for two sites in Eastern Canada. Sol. Energy 137 (Nov.), 46–54 http://linkinghub.elsevier.com/retrieve/pii/ S0038092X16303000. Gotoh, K., Fujii, Y., 1998. A fractal dimensional analysis on the cloud shape parameters of cumulus over land. J. Appl. Meteorol. 37 (10), 1283–1292 http:// journals.ametsoc.org/doi/abs/10.1175/1520-0450(1998)037%3C1283% 3AAFDAOT%3E2.0.CO%3B2. Hammer, A., Heinemann, D., Hoyer, C., Kuhlemann, R., Lorenz, E., Müller, R., Beyer, H.G., 2003. Solar energy assessment using remote sensing technologies. Remote Sensing Environ. 86 (3), 423–432 http://www.sciencedirect.com/science/ article/pii/S003442570300083X. Hinkelman, L.M., 2013. Differences between along-wind and cross-wind solar irradiance variability on small spatial scales. Sol. Energy 88 (Feb.), 192–203 http://linkinghub.elsevier.com/retrieve/pii/S0038092X12004021. Hoff, T.E., Perez, R., 2010. Quantifying PV power output variability. Sol. Energy 84 (10), 1782–1793 http://www.sciencedirect.com/science/article/pii/S0038092 X10002380. Hoff, T.E., Perez, R., 2012. Modeling PV fleet output variability. Sol. Energy 86 (8), 2177–2189 http://linkinghub.elsevier.com/retrieve/pii/S0038092X11004154. Kalisch, J., Macke, A., 2008. Estimation of the total cloud cover with high temporal resolution and parametrization of short-term fluctuations of sea surface insolation. Meteorol. Zeitschr. 17 (5), 603–611 http://openurl.ingenta. com/content/xref?genre=article&issn=0941-2948&volume=17&issue=5&spage= 603. Karperien, A., 2015. FracLac for ImageJ. . Kühnert, J., Lorenz, E., Heinemann, D., 2013. Satellite-based irradiance and power forecasting for the German energy market. In: Solar Energy Forecasting and Resource Assessment. Boston, pp. 267–297. Klima, K., Apt, J., 2015. Geographic smoothing of solar PV: results from Gujarat. Environ. Res. Lett. 10 (10), 104001 http://stacks.iop.org/1748-9326/10/i=10/a= 104001?key=crossref.acd1bfffc3bd66821eb229ff87c4e014. Lave, M., Kleissl, J., 2013. Cloud speed impact on solar variability scaling – application to the wavelet variability model. Sol. Energy 91 (0), 11–21 http:// www.sciencedirect.com/science/article/pii/S0038092X13000406. Lave, M., Kleissl, J., Arias-Castro, E., 2012. High-frequency irradiance fluctuations and geographic smoothing. Sol. Energy 86 (8), 2190–2199 http://linkinghub. elsevier.com/retrieve/pii/S0038092X11002611. Lave, M., Kleissl, J., Stein, J.S., 2013. A wavelet-based variability model (WVM) for solar PV power plants. IEEE Trans. Sustain. Energy 4 (2), 501–509 http:// ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=6269913. Lohmann, G.M., Monahan, A.H., Heinemann, D., 2016. Local short-term variability in solar irradiance. Atmos. Chem. Phys. 16 (10), 6365–6379 http://www.atmoschem-phys.net/16/6365/2016/.

264

G.M. Lohmann et al. / Solar Energy 150 (2017) 255–264

Lonij, V.P., Brooks, A.E., Cronin, A.D., Leuthold, M., Koch, K., 2013. Intra-hour forecasts of solar power production using measurements from a network of irradiance sensors. Sol. Energy 97 (Nov.), 58–66 http://linkinghub.elsevier.com/ retrieve/pii/S0038092X13003125. Lorenz, E., Hurka, J., Heinemann, D., Beyer, H.G., 2009. Irradiance forecasting for the power prediction of grid-connected photovoltaic systems. IEEE J. Select. Top. Appl. Earth Observ. Remote Sensing 2 (1), 2–10 http://ieeexplore.ieee.org/ lpdocs/epic03/wrapper.htm?arnumber=4897348. Lovejoy, S., Schertzer, D., 2006. Multifractals, cloud radiances and rain. J. Hydrol. 322 (1), 59–88 http://www.sciencedirect.com/science/article/pii/S002216940 500106X. Lovejoy, S., Schertzer, D., Silas, P., Tessier, Y., Lavallée, D., 1993. The unified scaling model of atmospheric dynamics and systematic analysis of scale invariance in cloud radiances. In: Annales Geophysicae, vol. 11. pp. 119–127. Madhavan, B.L., Kalisch, J., Macke, A., 2016. Shortwave surface radiation network for observing small-scale cloud inhomogeneity fields. Atmos. Measure. Tech. 9 (3), 1153–1166 http://www.atmos-meas-tech.net/9/1153/2016/. Mandelbrot, B.B., Wallis, J.R., 1968. Noah, Joseph, and operational hydrology. Water Resour. Res. 4 (5), 909–918 http://dx.doi.org/10.1029/WR004i005p00909. Mills, A., 2010. Implications of wide-area geographic diversity for short-term variability of solar power. . Perez, R., David, M., Hoff, T.E., Jamaly, M., Kivalov, S., Kleissl, J., Lauret, P., Perez, M., 2016a. Spatial and temporal variability of solar energy. Found. TrendsÒ Renew. Energy 1 (1), 1–44 http://www.nowpublishers.com/article/Details/REN-006. Perez, R., Ineichen, P., Moore, K., Kmiecik, M., Chain, C., George, R., Vignola, F., 2002. A new operational model for satellite-derived irradiances: description and validation. Sol. Energy 73 (5), 307–317 http://www.sciencedirect.com/science/ article/pii/S0038092X02001226. Perez, R., Kivalov, S., Schlemmer, J., Hemker Jr, K., Hoff, T.E., 2012. Short-term irradiance variability: preliminary estimation of station pair correlation as a function of distance. Sol. Energy 86 (8), 2170–2176 http://www. sciencedirect.com/science/article/pii/S0038092X12000928. Perez, R., Rábago, K.R., Trahan, M., Rawlings, L., Norris, B., Hoff, T., Putnam, M., Perez, M., 2016b. Achieving very high PV penetration – the need for an effective electricity remuneration framework and a central role for grid operators. Energy Policy 96 (Sep.), 27–35 http://linkinghub.elsevier.com/retrieve/pii/ S0301421516302452. Perpiñán, O., Lorenzo, E., 2011. Analysis and synthesis of the variability of irradiance and PV power time series with the wavelet transform. Sol. Energy 85 (1), 188– 197 http://www.sciencedirect.com/science/article/pii/S0038092X10002811. Perpiñán, O., Marcos, J., Lorenzo, E., 2013. Electrical power fluctuations in a network of DC/AC inverters in a large PV plant: relationship between correlation, distance and time scale. Sol. Energy 88 (0), 227–241 http:// www.sciencedirect.com/science/article/pii/S0038092X12004197.

Piacentini, R.D., Salum, G.M., Fraidenraich, N., Tiba, C., 2011. Extreme total solar irradiance due to cloud enhancement at sea level of the NE Atlantic coast of Brazil. Renewable Energy 36 (1), 409–412 http://www.sciencedirect. com/science/article/pii/S096014811000265X. R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. . Rasband, W.S., 2016. ImageJ. . Rusen, S.E., Hammer, A., Akinoglu, B.G., 2013. Estimation of daily global solar irradiation by coupling ground measurements of bright sunshine hours to satellite imagery. Energy 58 (Sep.), 417–425 http://linkinghub.elsevier.com/ retrieve/pii/S036054421300491X. Sachs, D., Lovejoy, S., Schertzer, D., 2002. The multifractal scaling of cloud radiances from 1m to 1km. Fractals 10 (03), 253–264 http:// www.worldscientific.com/doi/abs/10.1142/S0218348X02001385. Schmidt, T., Kalisch, J., Lorenz, E., Heinemann, D., 2016. Evaluating the spatiotemporal performance of sky-imager-based solar irradiance analysis and forecasts. Atmos. Chem. Phys. 16 (5), 3399–3412 http://www.atmos-chemphys.net/16/3399/2016/. Schneider, C.A., Rasband, W.S., Eliceiri, K.W., 2012. NIH Image to ImageJ: 25 years of image analysis. Nat. Meth. 9 (7), 671–675 http://dx.doi.org/10.1038/nmeth. 2089. Skartveit, A., Olseth, J., 1992. The probability density and autocorrelation of shortterm global and beam irradiance. Sol. Energy 49 (6), 477–487 http:// www.sciencedirect.com/science/article/pii/0038092X92901554. Solar Power Europe (SPE), 2016. Solar Market Report & Membership Directory – 2016 Edition. . Stetz, T., von Appen, J., Niedermeyer, F., Scheibner, G., Sikora, R., Braun, M., 2015. Twilight of the grids: the impact of distributed solar on Germany’s energy transition. Power Energy Mag., IEEE 13 (2), 50–61. Widén, J., 2015. A model of spatially integrated solar irradiance variability based on logarithmic station-pair correlations. Sol. Energy 122 (Dec.), 1409–1424 http:// linkinghub.elsevier.com/retrieve/pii/S0038092X15005940. Widén, J., Carpman, N., Castellucci, V., Lingfors, D., Olauson, J., Remouit, F., Bergkvist, M., Grabbe, M., Waters, R., 2015. Variability assessment and forecasting of renewables: a review for solar, wind, wave and tidal resources. Renew. Sustain. Energy Rev. 44 (Apr.), 356–375 http://linkinghub.elsevier.com/retrieve/pii/ S1364032114010715. Yordanov, G., Midtgärd, O.-M., Saetre, T., Nielsen, H., Norum, L., 2013. Overirradiance (cloud enhancement) events at high latitudes. IEEE J. Photovolt. 3 (1), 271–277.