Variability of interconnected wind plants: correlation length and its dependence on variability time scale

The variability in wind-generated electricity complicates the integration of this electricity into the electrical grid. This challenge steepens as the percentage of renewably-generated electricity on the grid grows, but variability can be reduced by exploiting geographic diversity: correlations between wind farms decrease as the separation between wind farms increases. But how far is far enough to reduce variability? Grid management requires balancing production on various timescales, and so consideration of correlations reflective of those timescales can guide the appropriate spatial scales of geographic diversity grid integration. To answer ‘how far is far enough,’ we investigate the universal behavior of geographic diversity by exploring wind-speed correlations using three extensive datasets spanning continents, durations and time resolution. First, one year of five-minute wind power generation data from 29 wind farms span 1270 km across Southeastern Australia (Australian Energy Market Operator). Second, 45 years of hourly 10 m wind-speeds from 117 stations span 5000 km across Canada (National Climate Data Archive of Environment Canada). Finally, four years of five-minute wind-speeds from 14 meteorological towers span 350 km of the Northwestern US (Bonneville Power Administration). After removing diurnal cycles and seasonal trends from all datasets, we investigate dependence of correlation length on time scale by digitally high-pass filtering the data on 0.25–2000 h timescales and calculating correlations between sites for each high-pass filter cut-off. Correlations fall to zero with increasing station separation distance, but the characteristic correlation length varies with the high-pass filter applied: the higher the cut-off frequency, the smaller the station separation required to achieve de-correlation. Remarkable similarities between these three datasets reveal behavior that, if universal, could be particularly useful for grid management. For high-pass filter time constants shorter than about τ = 38 h, all datasets exhibit a correlation length &xgr; ?> that falls at least as fast as τ − 1 ?> . Since the inter-site separation needed for statistical independence falls for shorter time scales, higher-rate fluctuations can be effectively smoothed by aggregating wind plants over areas smaller than otherwise estimated.


Introduction
Low CO 2 emission footprints make wind and solar power attractive choices for future electricity needs. However, their natural variability is challenging for the electric grid, which requires instantaneous matching of generation and load on all time scales from one ac cycle, through operational scheduling horizons (a day or two), out to planning horizons of more than a decade (von Meier 2006). Variability can be reduced and the grid-integration challenge lessened by interconnecting renewable electricity generators distanced enough that their variation is not fully correlated (Thomas 1945, Kahn 1979. At a great enough distance, which depends on variability time scale as we see here, the wind-speed variations approach statistical independence: an ensemble of wind plants separated at least this distance from each other is 'geographically diverse,' and at the associated time scale, variability of the ensemble's summed power output will be reduced compared to that of an individual plant.
Better plans for future grids could be made if this effect were better understood. Many previous studies have relied on an empirical approach where the degree of smoothing that might have occurred over a specific region was estimated from historical meteorological or reanalysis data or even calculated from actual renewable generation (Molly 1977, Palutikof et al 1990, Landberg 1997, Milligan and Factor 2000, Wiemken et al 2001, Archer and Jacobson 2003, 2007, Holttinen 2005, Sinden 2007, Kiss and Jánosi 2008, Kempton et al 2010, Fertig et al 2012, Cosseron et al 2013, Fisher et al 2013, Huang et al 2014, Louie 2014. However, it is often difficult to generalize the findings from the study scenario to realistic planning scenarios. As a first step, some have sought to model how correlation falls off with distance, using a variety of different forms: exponential, both with (Haslett and Raftery 1989, Landberg 1999, Kiss and Jánosi 2008, and Katzenstein et al 2010 and without (Giebel 2000, Simonsen and Stevens 2004, Holttinen 2005, Achberger et al 2006, and Kempton et al 2010 a nugget (non-unity correlation at zero separation), stretched exponentials (Hasche 2010, Louie 2014), Gaussian forms (Buell 1972) and Cauchy forms (Buell 1972, Julian andThiebaux 1975); see table 1. A few have gone further to propose forms predicting probability distributions of aggregated power as a function of the region size (Justus and Mikhail 1978, Carlin and Haslett 1982, Haslett and Raftery 1989, Hasche 2010. How correlations depend on the variability time scale is equally important to how they depend on distance, since, for a given magnitude operational power excursion, the faster the excursion generally the more costly its regulation (Kirby 2004). Without suggesting a particular functional form for the dependence, studies by Ernst et al (1999) and by Mills and Wiser (2011) have found for wind and solar, respectively, that faster variations become uncorrelated at smaller spatial separations than slower variations. Variability time scale τ or frequency f can be inserted into functional forms for correlation ρ versus distance r by introducing a characteristic velocity v, and replacing correlation length ℓ by vτ or v/f, as proposed by Davenport (1961). Beyer et al (1990Beyer et al ( , 1993, McNerney and Richardson (1992), Nanahara et al (2004), and Sørensen et al (2008) have taken this approach for wind-speed correlation. Calling r/(vτ) a 'dispersion factor', Hoff and Perez (2010) introduced this concept to the study of solar variability and identified velocity v with cloud motion; it has since been used in solar correlation studies by Marcos et al (2012), Lave andKleissl (2013), andHinkelman (2013). None of the studies based on characteristic-velocity functional forms have considered regions much larger than 100 km in extent, and more typically have been limited to the size of a single PV plant.
Here we revisit the methods used by Ernst et al (1999) and by Mills and Wiser (2011), apply them to both wind-speed data and wind generation data from three different geographical regions, and find a single quantitative relationship between correlation length and time scale that parsimoniously characterizes behavior on time scales ranging from an hour to a quarteryear and over distances of a kilometer to a continent.
In section 2, we highlight the unique features of the datasets used in this analysis as well as novel methods for investigating correlations between sites and quantifying correlation length scales. Section 3 describes the results found in our investigation and section 4 sums up our conclusions and suggestions for future work. In the provided supplementary data (stacks.iop. org/ERL/10/044004/mmedia), we describe, in detail, the filtering methodology used to eliminate non-stationarities in the time series and present the correlation results from all datasets.

Datasets
Although previous studies have characterized correlations within an individual region, the present study seeks universal behavior across multiple regions and datasets: an Australian wind generation dataset (AUS), a Canadian wind speed dataset (CAN) and a Bonneville Power Administration wind speed dataset (BPA). The three datasets presented here bring different features to the study, such as wind power production data (AUS), great extent (CAN) or fine time resolution (AUS, BPA). More details on these datasets, including maps of the stations used, appear in the supplementary data (stacks.iop.org/ERL/10/044004/mmedia).

Australian generation dataset (AUS)
The Australian Energy Market Operator (AEMO) provides a 1 year (October 2013-October 2014) dataset of five-minute electricity generation data from 32 wind farms across Southeastern Australia including the provinces of New South Wales, South Australia, Tasmania and Victoria (table 2). AEMO reported in 2013 that they might curtail wind power systems in the future (Australian Energy Market Operator 2013); therefore we suspect few curtailment effects in this dataset. We utilize a subset of 29 wind farms, combining production for different wind farm stages and selecting only the larger-producing member of wind farm pairs within 5 km of one another. We make no correction for wake effects of neighboring wind farms (Fitch et al 2013) in this analysis.
The terrain elevations of these wind farms range from 3 to 900 m above sea level (asl). Plant output data are normalized by each wind plant's generation capacity, which range from 20 to 420 MW. The closest pair of farms is separated by 8 km and the farthest by 1274 km.

Canadian wind speed dataset (CAN)
The National Climate Data Archive of Environment Canada provides a 54 year dataset of hourly surface wind-speeds from 117 stations ranging across more than 5000 km of Canada (table 2). While the entire dataset spans 1953 through 2006, missing data are less frequent in later years, so we use a 45 year subset of data, from 1962 to 2006. During the data collection, the measurement heights of the stations ranged from 5 to 32 m above ground level (agl), but Wan et al (2010) homogenized these data to the standard 10 m height assuming a logarithmic wind profile. The terrain elevations within this dataset vary from sea level to 1084 m asl. Separation distances between sites vary from a minimum of 19 km to a maximum of 5344 km.

BPA wind speed dataset
The BPA provides a dataset of five-minute windspeeds from meteorological towers in Washington and Oregon in the Northwestern US (table 2). For this analysis, we use data from 14 sites from 2010 to 2014 due to the higher temporal resolution of these data, although several stations provide data for longer durations with lower time resolution (BPA). These 14 towers in the BPA network are generally located along the Columbia River Gorge on the border of Washington and Oregon or along the Pacific Coast. The elevations of these stations range from 19 to 1261 m asl. The anemometer heights range from 9 to 53 m agl. No 'standard height' homogenization methodology is used on this dataset; correlation coefficient calculations between sites would be unchanged if either a logarithmic wind profile was assumed or the wind power law was used to homogenize the measurements to a 10 m height. Separation distances between sites vary from a minimum of 13 km to a maximum of 354 km.

Data processing 2.2.1. Pre-processing
Although both deterministic variability, such as diurnal or seasonal cycles, and random or stochastic variability will figure in any complete accounting of wind power variability, it is the stochastic phenomena that are less understood and not as predictable, so we focus on them here. The analysis of underlying stochastic components is complicated by the presence of temporal periodicities, as noted by Haslett and Raftery (1989), Gunst (1995), Robeson and Shein (1997), Achberger et al (2006), andHill et al (2012). To remove the diurnal cycle in order to more clearly characterize stochastic spatial correlations, we make a 'local' estimate of the amplitudes of the first four daily harmonics by least-squares fitting using a 90 day moving window to allow for seasonal variation in the daily cycle (Baïle et al 2011). This cycle, estimated once for each of the days in the dataset, is subtracted from that day's wind-speed or power data, as further described in the supplementary data. By also subtracting the 90 day moving average, we remove lowfrequency variability, including the seasonal or annual cycle.
As shown by previous investigators (Ernst et al 1999), rapid wind-speed or power variations become uncorrelated at smaller spatial separations than do slow variations. In order to further investigate this effect, we prepare versions of each time series with trends slower than a chosen time-constant removed by calculating the average of the wind speed or power data over a segment centered at each time point and subtracting the segment average from the value at the center point. By varying the width of the segment ('window width'), we control the corner frequency of this high-pass filtering operation. Our process, described further in the supplementary data (stacks.iop.org/ ERL/10/044004/mmedia), differs from that used by Ernst et al (1999), Mills andWiser (2011), andHinkelman (2013) in that it subtracts from the original data an averaged version of that data, rather than calculating 'ramps' by differencing one block average value from the next. Thus, the number of data points in the filtered versions of our time series is reduced from the original only by the width of filter window, rather than being decreased by a factor of the window width. We analyze the effects of high-filtering for window widths τ ⩽ 1049 h, and denote by τ = 2160 h data from which the 90 day seasonal bias term has been subtracted but no further filtering applied.

Correlation length
After pre-processing the data as described above and discussed in more detail in the supplementary data (stacks.iop.org/ERL/10/044004/mmedia), we calculate correlations between each pair of stations for each high-pass filtered version of the data. The scatter plot in figure 1, for example, shows the correlation data for CAN at filter window τ = 65 h. Given the high degree of scatter shown here, estimating a single correlation length for each dataset at each filtering window τ is problematic since a functional form of correlation versus distance appropriate for fitting is not known a priori. Wishing to avoid the functional-form conundrum, we sought a 'non-parametric' distance measure. For a given averaging window width τ we order the correlation coefficient values ρ ij by increasing station separation distance, with increasing values of separation distance signified with single index as r k , r 1 denoting the smallest separation (the closest pair of stations) and r N denoting the largest separation (the farthest pair). We then numerically 'integrate' the correlation data over distance using the trapezoid rule: ( 1 ) n k n k k r r 1 1 k k1 In the familiar case ρ(r i ) = exp(−r i /ℓ) + ε i , where ε i is mean-zero noise, this procedure gives ξ = ℓ in the limit of large r. For our observations, we take r 0 ≡ 0, and estimate ρ r 0 as the intercept of a local regression curve (see section 2.2.3 below). As shown by the light blue line in figure 1, ξ(r n ) varies smoothly with r n and approaches a near-constant value at the largest separations (here 139 km at r N = 5300 km) as the ρ r k scatter there around zero. We refer to the value at the greatest separation ξ(r N ) (at a given high-pass filter window width τ hereinafter ξ τ ) as the correlation length; it has a meaning similar to Bretherton's 'correlation radius' (1999) or the term 'integral scale' in turbulence (Batchelor 1953). Our methodology bears some resemblance to that of Şen's cumulative semivariogram (Şen and Şahin 1997, Şen 1989).

Local regression curves
To facilitate graphical comparisons of correlation versus distance behaviors that might otherwise be obscured by the high degree of scatter in the individual ρ(r k ) values, we also calculate ρ(r) curves using local regression techniques (Cleveland 1979), where, to find the curve value at a given site-separation distance r′, we fit a 2nd-order polynomial to the fraction α of the correlation data points that are closest to the given distance (for example in the CAN dataset and choosing α = 0.05, the 339 correlation data points out of 6786 total having the smallest values of − ′ r r ) . This gives results as exemplified by the thick blue line in figure 1.

Importance of pre-processing
The high-pass filtering and diurnal-cycle removal (outlined in detail in the supplementary data (stacks.iop.org/ ERL/10/044004/mmedia) have significant effects on the correlation behavior, as shown in figure 2, which, to facilitate comparison, portrays correlation versus distance as local regression curves. With no pre-processing (raw data, dashed brown curve), correlation doesn't fall below about 0.04 at the largest site separations, while with the diurnal cycle and seasonal bias removed (solid brown curve), correlation falls to zero. For raw data high-pass-filtered with τ = 33 h (aqua), the correlation does not fall below about 0.10 (dashed curve) for separations smaller than 4000 km, while if the high-pass filter is applied to data from which the diurnal cycle has first been removed, correlation falls to near zero for separations of 600 km (solid curve). High-pass cutoffs of τ = 25-37 h give similarly high correlation floors on raw data while the height of the floor decreases noticeably for cutoffs smaller than 21 h or larger than 41 h.

Correlations
Despite the differences in methodology, the behavior of correlation with distance and high-pass filtering window width qualitatively resembles the previous wind-power results of Ernst et al (1999) and insolation work of Mills and Wiser (2011): correlation falls off more rapidly with site separation the smaller the window-width τ, as seen for the CAN data in figure 3. Interestingly, we observe that for some high-pass filter window widths (for example, a 65 h window width in CAN), the correlations actually become negative between 500 and 1000 km separation distances, perhaps for reasons similar to those identified for solar correlations by Hoff and Perez (2010) and Hinkelman (2013).
Substantial scatter of the correlation versus distance values is a prominent feature of data from all three regions at all but the shortest filter widths (see supplementary figures S5-7). We explored several potential causes of this large scatter as well as the negative correlations by separating stations by region (North versus South as well as East versus West) and by azimuthal bearing. East/West differences (figure 4(a)) and North/South differences ( figure 4(b)) explain neither the scatter nor the slight anti-correlations. It is clear from previous work (Buell 1972, Ramanathan et al 1973, Julian and Thiebaux 1975) that the longitudinal and transverse horizontal wind components have different correlation behaviors, with the transverse correlation falling faster and exhibiting more negative values, similar to the behavior of fully developed turbulence (von Kármán 1948, Batchelor 1953. Unfortunately, the CAN dataset provides wind speed but not wind direction. Nevertheless, some of the scatter in the CAN correlations can be attributed to the azimuthal bearing of each station pair, as stations separated along a line 10°N orth of West-East are systematically less correlated and more often anti-correlated than those perpendicular to that line ( figure 4(c)). Šaltytė Benth and Šaltytė (2011) also observed directional anisotropy in wind-speed correlation fall off. We presume the anisotropy we observe arises from the prevailing westerly winds produced by the large-scale circulation in this region, although, contrary to our expectation, we see higher correlations for stations separated along the crosswind direction.  Finally, given previous work that has identified a connection between the interannual climate oscillation of the El Niño Southern Oscillation and Canadian wind resources (St. George and Wolfe 2009), we separated the time series into intervals with strong climate signals (1988-1989La Niña and 1982-1983El Niño, 1964-1965La Niña and 1965-1966El Niño, 1970-1971La Niña and 1972-1973, and calculated each site-pair's correlation coefficient over the El Niño interval and over the La Niña interval. Comparison of all these periods gave results similar to that seen in figure 4(d) for the 1964-1966 periods. No difference in correlations between El Niño and La Niña periods emerges ( figure 4(d)). The pairwise differences between the '65-66 El Niño and the '64-65 La Niña correlations have mean |μ| = 0.0011 and standard deviation 0.044 (see figure S8); the non-parametric sign test of null hypothesis H 0 that μ = 0 against alternative hypothesis H 1 that μ ≠ 0, fails to reject H 0 with a p-value of 0.12. Figure 5 shows the results of our 'integral-scale' calculations for CAN, AUS and BPA. The 5300 km geographic extent of Canada clearly exceeds the observed correlation lengths, allowing the integration to fully saturate for even the widest filter, giving ξ 2160h = 273 km. The AUS data yield slightly larger values, with ξ 2160h = 368 km. The geographic extent of the BPA region is not sufficient for the correlation integration to saturate, making the ξ 2160h = 89 km value, obtained an underestimate of the full correlation length. Figure 6 shows the variation of correlation length ξ τ with high-pass filter cutoff τ for the three regions, the ξ values for each region normalized by the maximum to facilitate comparison between datasets. The envelope surrounding the correlation length values as a function of τ indicates an empirically-found range of behavior common to all three datasets:

Discussion
The characteristics of the correlation length-scale metric introduced in equation (1) and used above deserve further explication. This measure has the advantage that it is essentially free of assumptions about the functional form of the correlation versus distance data. It quantifies the site-separation distance needed to reduce average inter-station correlation to a small value. However, it is important to remember, especially with data with such large scatter as observed here, that ξ τ is just one of multiple possible measures of 'how far is far enough'. For example, as seen in figure 1 for the CAN data, numerically integrating the scattered ρ(r) points according to equation (1) gives ξ 65h = 139 km (righthand end of light-blue curve), while integrating the dark-blue robust local-regression curve for the same window width gives the significantly larger value of 172 km (integral not shown). This occurs because the residuals of the local-regression fit are negatively skewed, and the 'robust' fitting process (Cleveland 1979) pulls the fit towards the median by assigning small weights to points determined to be outliers; here more often below the curve than above. The integral-scale metric ξ τ measures the separation distance required for correlation to fall to a value small compared to unity; this can be substantially less than the distance ℓ over which it falls to a fraction (say, 1/e) of its initial value (at r = 0) if the initial value is small. This discrepancy between the integral-scale metric and the 1/e distance ℓ can be seen from the parameters of fitting of an exponential form to the data, as shown by the dashed-lines in figures 5(b) and (c). Here we fit ξ(r) = βℓ(1 − e −r/ℓ ) as would be the case for correlation falling as βe −r/ℓ with distance.
Fitting results are shown in table 3. Due to the large nugget effect (Matheron 1963) observed in all three regions, the best-fit correlation at zero site-separation β is substantially smaller than unity, proportionately reducing the area ξ under the ρ(r) curve. For the CAN data this gives ξ(r N ) ≈ βℓ, but in the AUS and BPA regions ξ(r N ) < βℓ since the regions sizes do not permit separations r large enough to bring ρ(r) to zero and hence to bring ξ to saturation.
The nugget effect has a particularly strong influence on ξ τ values at small τ. Had we constrained ρ(0) = 1 rather than estimating the value at the origin using local regression, figure 6 would have shown the decline in ξ τ / ξ 2160 ending at a floor value limited by our trapezoidrule integration to half the separation of the closest site pair (0.03, 0.01, and 0.07 for CAN, AUS, and BPA,

Conclusions
Variability in wind-generated electricity can be reduced by interconnecting wind farms across large regions, as distant wind farms are only weakly correlated. We investigate the geographic extent needed for this aggregation to be effective using three extensive data sets. Using a 'non-parametric' estimator ξ based on the area under the correlation versus distance curve, we found for the slowest variabilities ξ = 273 km max for 10 m wind speeds in Canada, 368 km for wind-plant generation in Southeastern Australia, and 89 km for tower wind-speeds in the Northwestern US. Since the Australia and BPA regions are small enough that even for the most distant sites correlation never drops to zero, our ξ values for the widest filter window-widths are underestimates of the extent needed for fully effective aggregation. Quantities more representative of the extent needed for smoothing the slowest variabilities are given by the larger βℓ values of 447 and 130 km respectively, in table 3, obtained by fitting ξ(r) to the form βℓ(1 − e −r/ℓ ) and extrapolating ξ to r = ∞. These values can be compared to like values from previous work listed in table 1. Although the regional length scales have different magnitudes, we find a dependence of ξ on variability time-scale that is remarkably similar across the three regions, as seen in figure 6. At time scales τ shorter than 38 h, ξ falls at least as fast as τ −1 , while at longer scales it is essentially constant. Thus, on time scales longer than a day or so, the variability-reduction benefit of aggregating wind power over a region of a given size will be independent of time scale. For time scales shorter than a day, the faster the variability, then the more smoothing that region could provide. It is the shrinking of correlation length with time scale that gives high-frequency spectral slopes a larger magnitude for power aggregated over a region compared to a single site (Beyer et al 1990, McNerney andRichardson 1992 Our findings help disentangle the effects on variability reduction of generator number, region size, and variability time scale. In general, aggregating the outputs of N uncorrelated generators should reduce variability magnitudes by a factor of N . A geographic region of area A has roughly N ≈ A/(2ξ) 2 potential uncorrelated sites; thus, for a fully populated region (Hasche 2010) variability could be attenuated by a factor up to ξ A /(2 ). For time scales shorter than τ ≈ 1.5 days, the number of potential uncorrelated sites within a fixed area, and hence achievable attenuation, grows at least as fast as 1/τ. If the actual number of generators in a region is less than A/(2ξ) 2 , the lesser number will determine attenuation.
Further work to better understand whether the high degree of scatter in correlation versus distance is a stable manifestation of some unidentified geographic process or is just persistent, random temporal variation that would average away over longer records would improve model utility. Additionally, analysis of solar data over a large region could determine if time and length scale of solar variability are linked in a way similar to what we found here for wind.