Multiscaling and joint multiscaling description of the atmospheric wind speed and the aggregate power output from a wind farm

We consider here wind speed time series and the aggregate output wind power from a wind farm. We study their scaling statistics in the framework of fully developed turbulence and Kolmogorov’s theory. We estimate their Fourier power spectra and consider their scaling properties in the physical space. We show that the atmospheric wind speed and the aggregate power output from a wind farm are intermittent and multifractal over a wide range of scales. The coupling between simultaneous data of the wind speed and aggregate power output is investigated through a joint multifractal description using the generalized correlation functions (GCFs). This multiscaling test is compatible with a linear relation between the wind speed and the aggregate power output fluctuations for timescales T > 103s' 15min.


Introduction
Increasing the wind energy contribution to electrical networks requires improving the tools to forecast the electrical power produced by wind farms, in order to proportion network lines.However the wind energy is a fluctuating energy ressource, due to the high variability of the wind at all spatial and/or temporal scales.In the atmospheric boundary layer, the Reynolds number (ratio of inertial to viscous force) can be as large as Re = 10 8 (Burton et al., 2001).Large values of the Reynolds number lead to a huge intermittency of wind speed fluctuations at all temporal or spatial scales, ranging from large-scale variations (years) to very small scale variations (few min down to seconds) (Stull, 1988).Small-scale intermittency remains a challenging problem for the turbulence community research (Peinke et al., 2004).Several ap-proaches can be used to consider the scaling intermittency of small-scale turbulence, the most classical one being structure functions analysis.In recent years several studies have been dedicated to the analysis of scaling laws and turbulent intermittency at small scales in the laboratory (Anselmet et al., 1984;She and Levêque, 1994) and in the atmospheric boundary layer (Schmitt et al., 1993(Schmitt et al., , 1994;;Katul et al., 1995;Schmitt, 2007;Böttcher et al., 2007;Morales et al., 2011;Calif and Schmitt, 2012).These studies have shown that atmospheric turbulent speed at small scales has multifractal scaling fluctuations and exhibits long-range power correlations.
However the knowledge of variations ranging from minutes to a few days -corresponding to 1-1000 km, i.e., the mesoscale range -is necessary to provide efficient tools for management and control of wind power generation.The studies concerning this scale range are fewer than those for the small-scale range, due to the possible nonuniversality of the power law slope in the mesoscale range.Recent works (Lauren et al., 1999;Muzy et al., 2010) have been dedicated to scaling and multiscaling properties of the atmospheric wind speed in the mesoscale range.
In this paper, the scaling properties of the atmospheric wind speed are provided for the small-scale and mesoscale ranges.Atmospheric wind speed data sampled at 20 and 1 Hz are analyzed through multifractal theory in order to characterize the wind speed fluctuations behavior for each scaling regime.In parallel, the aggregate power output data from a wind farm are analyzed within the same theoretical framework.The paper is organized as follows.In Sect.2, the data sets are described.In Sect.3, the theoretical framework, traditional spectral analysis and structure functions analysis are Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
presented.In Sect.4, in order to quantify the coupling of the relationship between the wind atmospheric and the aggregate power output, a joint multifractal description is presented.In Sect.5, the Fourier spectra are estimated for the atmospheric wind speed and the aggregate power output.To provide a full characterization of experimental fluctuations, at all intensities and at all scales, multifractal and joint multifractal approaches are used.

Description of the data
In order to investigate the spectral and multiscaling properties of the turbulent wind speed v in the atmospheric surface layer and the aggregate power output P , a data set is collected on the wind energy production site of Petit-Canal, an island in the French West Indies.This 10 MW production site, located at 16 • 15 N latitude and 60 • 30 W longitude, was positioned at approximately 60 m (197 ft) a.s.l., at the top of a sea cliff.The wind speed v is measured with a three-cup anemometer (model A100L2 from Vector Instruments) black having a response time of 0.15 s and an ultrasonic anemometer (model CSAT3).Both were mounted on a 40 m (131 ft) tall mast erected 20 m (66 ft) from the cliff edge, at 38 m (125 ft) from the ground.The two experimental databases for atmospheric wind speed are used here: (i) with the threecup anemometer, the sampling frequency and duration are respectively 1 Hz and 1 yr (January 2006(January -2007)); (ii) with the ultrasonic anemometer, the sampling frequency and duration are respectively 20 H z and 2 weeks (14-28 July 2005).Moreover, the aggregate power output produced by this wind farm was recorded with a sampling frequency 1 Hz, during 1 yr (January 2006(January -2007)).Table 1 gives a description of our data set with the mean and the standard deviation values.Additionally, two examples of simultaneous wind speed and aggregate power output data sequence for different timescales, are illustrated in Fig. 1: (i) Fig. 1a presents simultaneous time series sampled at 1 Hz; (ii) Fig. 1b presents the moving averages of simultaneous time series, over a time period T = 10 4 s. Figure 1a shows that for very short timescales of the order of seconds, the wind speed signal exhibits high fluctuations contrary to the wind power signal, while over large scales the observed moving averages have the same trends.Classically, a scale invariance can be detected by computing of power spectral density E(f ).For a scale-invariant process, the following power law is obtained over a range of frequency f : where β is the spectral exponent.According to some authors (Mandelbrot, 1982;Schertzer and Lovejoy, 1987;Marshak et al., 1994), it contains information about the degree of stationarity of the field: β < 1, the process is stationary, β > 1, the process is nonstationary, -1 < β < 3, the process is nonstationary with increments stationary.
It may also be considered as characterizing the degree of correlation (Ivanov et al., 2001;Telesca et al., 2003).
In the atmospheric surface layer, it is recognized that wind speed spectra possess three spectral regions (Pope, 2000): (1) the energy-containing range black (or large-scale turbulence) in which turbulent kinetic energy is produced, (2) the inertial subrange (small-scale turbulence) in which turbulent kinetic energy is transferred from large to small scales, and (3) the dissipation range (millimeters and smaller) where turbulent kinetic energy is converted to heat by the action of fluid viscosity.Figure 2 gives an illustration of the Kolmogorov quasi-equilibrium energetic cascade in the inertial range.In this zone, the Kolmogorov theory implies that the wind speed fluctuations possess a power law spectrum (Kolmogorov, 1941;Obukhov, 1941): This relation is written here for wave numbers, but can also be used for frequencies, involving Taylor's hypothesis to relate spatial fluctuations to temporal fluctuations.This power law has been verified many times for wind atmospheric data (Schmitt et al., 1993(Schmitt et al., , 1994;;Katul et al., 1995;Schmitt, 2007;Morales et al., 2011;Calif and Schmitt, 2012).In the energycontaining range where the scales extend from a few minutes to a few days, the properties of turbulent wind speed depend on their strong anisotropy and their dependence on the flow domain boundary characteristics (Katul and Chu, 1998;Katul et al., 2012).However in this study, the existence of universal power laws at low wave numbers for the power spectral density of the turbulent longitudinal velocity is examined theoretically and experimentally for the near-neutral atmospheric surface layer.Indeed, several experimental studies have found a −1 power law slope at production scales for the wind speed spectrum in the atmospheric surface layer in neutral conditions (Katul and Chu, 1998;Nickels et al., 2005;Katul et al., 2012): Moreover this power law slope was predicted by three theoretical approaches.More details are given in Katul and Chu (1998).However this power law does not present the same level of universality as β = 5/3 for turbulence in the inertial range: as mentioned in Lauren et al. (1999), in the mesoscale range, the power law slope can vary with the local topography and the atmospheric conditions.Let us note also that in the wind energy community some fits of speed power spectra are often used, based on von Karman's formula of the form E(f ) = a/(b + cf 2 ) 5/6 (Karman, 1948) or on Kaimal's formula E(f ) = A/(1 + Bf ) 5/3 (Kaimal et al., 1972), where a, b, c and A, B are constants.These two formulae have no theoretical ground: von Karman wanted to perform an inter-  (1972) proposed their formula on purely empirical grounds for the low-frequency part.In both cases the main idea is to capture in a single expression the injections scales and the inertial range with Kolmogorov spectrum.Indeed, these fits are based on Kolmogorov theory (Simiu and Scalan, 1978) which defines a general form of the fluctuating wind speed shown as follows (Zhang et al., 2008).In the present paper, we mainly focus on other scaling regimes and therefore do not need to consider the fits proposed by these authors.Finally, we must note also that the power spectrum density E(f ) is a second-order statistic (proportional to the square of the amplitude of a given frequency fluctuation) and its slope is not sufficient to fully specify a scaling process.Multifractal analysis is a natural generalization to study the scaling behavior of a nonlinear phenomenon, using qth-order structure functions, and to obtain a full characterization of wind speed fluctuations, at all intensities and all scales; this is discussed in the next subsection.

Intermittency and multifractal properties
Intermittency in turbulence has been a subject of research for almost 50 yr now, following Kolmogorov's 1962 seminal works (Kolmogorov, 1962).This intermittency in the inertial range has been modeled by cascade models in various studies for velocity (Novikov and Stewart, 1964;Yaglom, 1966;Mandelbrot, 1974).These models reproduce intermittency and multiscaling in the inertial range turbulence.Here we may define intermittency as the property of having large fluctuations at all scales, with a correlated structure, or, in other words, with various orders of singularities to be distributed over fractal sets with varying dimensions (rather than a unique dimension as for monofractal process); in the case of intermittency the scaling moment function ζ (q), which is introduced below, is nonlinear.Large fluctuations are much more frequent than what would be obtained for Gaussian processes (Frisch, 1995;Schertzer et al., 1997;Vulpiani and Livi, 2004).This is typically studied considering the probability density function (PDF), or more often using the moment of order q of these fluctuations, called "structure functions of order q": where X( t) = X(t + t) − X(t) and t is a time increment.We have here written the fluctuations in time, since below in this paper we deal with time series analysis.ζ (q) is the scaling exponent function.ζ (2) = β − 1 relates the secondorder moment to the power spectrum scaling exponent; the knowledge of the full (q, ζ (q)) curve for integer and noninteger moments provides a full characterization of wind speed fluctuations at all scales and all intensities.The parameter 1) is the Hurst exponent characterizing the nonconservation of the mean.Monofractal processes correspond to a linear function ζ (q) = qH , with H = 1/2 for the Brownian motion, 0 < H < 1 (H = 1/2) for a fractional Brownian motion (that can be defined as a fractional integration of order b (0 < b < 1) of a Gaussian noise, with H = 3 2 − b).The values of the function ζ (q) are estimated from the slope of S q ( t) versus t in a log-log diagram for all moments q.The function ζ (q) defines the types of scaling behavior; in other words, this exponent function is useful to characterize the statistics of the random process.If ζ (q) is linear, the statistical behavior is monoscaling; if ζ (q) is nonlinear and concave, the behavior is defined as mutiscaling, corresponding to a multlfractal process.The concavity of this function is a characteristic of the intermittency: the more the curve is concave, the more the process is intermittent (Frisch, 1995;Schertzer et al., 1997;Vulpiani and Livi, 2004).
Kolmogorov's 1941 model corresponds to a linear model for the exponent function (Kolmogorov, 1941): For a multifractal process, ζ (q) is nonlinear.Several models have been proposed in the literature to fit the scaling exponents ζ (q).Here we consider the classical lognormal model (Kolmogorov, 1962): with µ being the intermittency parameter.Yaglom (1966) has proposed a multiplicative cascade model compatible with Kolmogorov's ideas.Several other models have been proposed since in the literature, for example the log-stable model proposed by Schertzer and Lovejoy (1987): where H = ζ (1) is the Hurst parameter which defines the degree of smoothness or roughness of the field.The parameter C 1 is the fractal co-dimension of the set giving the dominant contribution to the mean (q = 1) and bounded between 0 and d (d the dimension space, here d = 1).It measures the mean intermittency characterizing the sparseness of the field: the larger C 1 , the more the mean field is inhomogeneous.The multifractal Lévy parameter α is bounded between 0 and 2, where α = 0 corresponds to the monofractal case and α = 2 coresponds to the multifractal lognormal case.The parameter α measures the degree of multifractility, i.e, how fast the inhomogeneity increases with the order of the moments.But here we consider the lognormal model which provides a reasonable fit up to q = 5; hence here the question of the best model (among the inifinitely divisible family of models) is not the topic of the present paper, and we consider here the lognormal fit as convenient for the joint analysis done in the next section.
4 Joint analysis for multivariate data

Test for data independence using second-order correlation
In order to test independence between two random process, second-order statistics are often considered in the wind energy community (Burton et al., 2001).More precisely the cross-correlation function in the time domain and the coherence function in the frequency domain, of two processes x(t) and y(t), are determined for highlighting possible correlations.We recall here the expressions for the cross-correlation function C xy and the coherence function H xy (Papoulis and Pillai, 2002): where < .> is the statistical average, τ is a time lag, σ x , σ y are respectively the standard deviations of processes x(t) and y(t) (Papoulis and Pillai, 2002), and These statistical tests are based on second-order statistics.In the following section a generalization of correlation coefficient is given using the multifractal framework.

Joint multifractal description for bivariate field:
generalized correlation functions (GCFs) and exponents (GCEs) We recall here that, whereas independence implies noncorrelation, noncorrelation does not imply independence.In order to better consider the relation between two scaling time series, we apply here a testing technique proposed in Seuront and Schmitt (2005).Instead of random variables x and y, we consider here the increments of two stochastic processes , and the normalization of the joint moments is given as (Seuront and Schmitt, 2005) When X h and Y g are independent, r(h, g) = 0 and c(h, g) = 1.In contrast, increasing values of c(h, g) would characterize increasing dependence between X h and Y g .The generalized correlation exponent (GCE hereafter), estimated as the slope of the power law of c(h, g) versus t in a log-log plot, is then expressed as: where ζ X (h) and ζ Y (g) characterize the multiscaling properties of the single fluctuations < X h >, and < Y g > and S(h, g) characterize the multiscaling properties of the joint fluctuations < X h Y g >.Both c(h, g) and r(h, g) are generalizations of correlation functions.This multiscaling test for independence between two stochastic processes takes into account the multifractal character of intermittent processes.It also allows to test for the phenomenology responsible for the high-intensity (rare and unexpected) fluctuations observed in intermittent distributions, considering their potential association with both high-and low-intensity fluctuations characterized respectively by high and low orders of moment.In Seuront and Schmitt (2005), GCEs are considered in special cases: -If X and Y are lognormal multifractal processes (Meneveau et al., 1990), then Consequently In this case, it is clear that r(1, 1) or r(2, 2) are enough to estimate the only needed parameter, namely the correlation coefficient σ , so that if r(1, 1) = 0 or r(2, 2) = 0, it can be concluded that the two processes are independent.
-If X and Y are independent, r(h, g) = 0; -if X and Y are proportional, i.e., X = aY , then Additionally the shape of the obtained surface is symmetric in the h-g plane.
-If a power law exists between X and Y , i.e., X = aY b , one has In this case r(h, g) > 0, but it is symmetric in the bh-g plane.

Spectral and multifractal analysis
In this section the Fourier power spectral densities are estimated for our database: wind speed data sampled at 20 and 1 Hz, and aggregate power output data sampled at 1 Hz.After this, the multifractal analysis is applied.

Spectral analysis
In this study, in order to estimate the power spectral densities of the wind speed and the aggregate power output from a wind farm, the discrete Fourier transform of the times series considered is computed.The expression of the power spectral density for a process x(t) is recalled here.A N point-long interval used to construct the value at frequency domain point f , X f (Bracewell, 1999): and the power spectral density is The Fourier spectra reveal two scaling ranges in the atmospheric wind speed and aggregate power output data.In Fig. 3, we have plotted the spectrum for the wind speed data sampled at 20Hz.This spectrum is obtained from an average of 2886 spectra computed with wind speed time series of length 1000 data points.This averaged spectrum highlights two spectral regimes.For the high frequencies 1 f 10 Hz, corresponding to timescales 0.1 T 1 s, in the inertial range, as expected the spectrum possesses a spectral slope β = 1.69 very close to 5/3 expected by Kolmogorov's theory.The slight difference is usually interpretated as coming from intermittency effects (Frisch, 1995;Schertzer et al., 1997).For the low-frequency part -i.e., f 1 Hz corresponding to timescales T 1 s in the energycontaining range -the spectrum follows a −1.28 spectral slope.This is consistent with Katul and Chu (1998) and recently Fitton et al. (2014).In the inset, we displayed the spectrum obtained from the wind speed data sampled at 1 Hz, exhibiting a −1.28 scaling for the low-frequency parti.e., 10 −7 f 1 Hz corresponding to timescales 1 T 10 7 s (10 7 s ≈ 4 months) -over seven decades.This con-firms the scaling obtained for frequencies f 1 Hz with the wind speed data sampled at 20 Hz.Furthermore this scaling is found also in Calif et al. (2014) using Hilbert spectral technique.Recent papers on longitudinal low-wave-number spectra provide evidence of a −1 power law in the atmospheric surface layer for neutral conditions, and a −5/3 power law for near-convective conditions (Katul et al., 1995;Katul and Chu, 1998;Kader and Yaglom, 1991).The power spectrum density E v (f ) presented in this paper seems to show a break around 1 Hz, in agreement with the results published in Lauren et al. (1999).The models of Kader and Yaglom (1991) and Katul et al. (1995) suggest a transition frequency of 0.3 Hz, and the Kaimal model a value of 1 Hz. Figure 3a shows that our present data are comparable with the latter value.
Figure 3b illustrates the average spectrum E p (f ) of the aggregate power output for the entire wind farm sampled at 1 Hz: the same slopes are observed.For frequencies 10 −4 f 0.5 Hz corresponding to timescales 2 T 10 4 s, E p (f ) displays a power law near the exact value 5/3, with β = 1.67.Previously published power spectra of wind generator power have shown power law regions of the power spectrum plot covering one or two decades of frequency (Sørensen et al., 2002;McNerney and Richardson, 1992).These studies have not provided a comparison between their data and the Kolmogorov spectrum.Apt (2007) has shown that the output wind power from a wind farm located follows a Kolmogorov spectrum over more than four orders of magnitude in frequency from 4.45×10 −6 to 3.33×10 −2 Hz, corresponding to timescales from 30 s to 2.6 days.Here, for frequencies f 10 −4 Hz corresponding to timescales T 10 4 s (approximately 3 h), a power law with β = 1.27 is observed for the first time.This highlights a scaling break approximately around 10 −4 Hz, for the power spectrum density E p (f ).
We obtain approximately the same value for the change of slope.Let us then remark here that wind speed and aggregate power output have the same type of regime with high-frequency −5/3 power law behavior and low-frequency −1.2 power law behavior.However the scale break for the regime change are not the same: T 0 = 1 s for the wind speed and T 0 = 10 4 s 3 h for the aggregate power output.We do not know whether this four-orders-of-magnitude change in the scale break is universal or is related to some spatial characteristics, such as local meteorology and spacing of wind turbine generators in the farm.
However, the power spectrum density E(f ) is a secondorder statistic (proportional to the square of the amplitude of a given frequency fluctuation) and its slope is not sufficient to fully specify a scaling process.Multifractal analysis is a natural generalization to study the scaling behavior of a nonlinear phenomenon, using qth-order structure functions; this is discussed below.

The structure functions scaling exponent ζ (q) for the atmospheric wind speed sequences
Here, we have analyzed two databases for atmospheric wind speed v sampled at 20 and 1 Hz.The structure function analysis is performed for wind speed increments v( t) = |v(t + t) − v(t)| for 0.05 t 500 s of data sampled at 20 Hz and for 1 t 6 × 10 5 s of data sampled at 1 Hz and for all moments between 0.15 and 5 with increments of 0.25.Figure 4a gives the scaling of the structure functions for a sequence of length 7000 data points for q = 1, 1.5, 2, and 2.5 in a log-log diagram for wind speed sampled at 20 Hz.The structure functions display two scaling tendencies: for 0.05 t 1 s corresponding to the inertial range and for 10 t 500 s. Figure 4b gives the scaling of the structure functions for a sequence for q = 1, 1.5, 2, and 2.5 in a log-log diagram for wind speed v sampled at 1 Hz, for 10 t 5 × 10 5 s.The straight lines in this figure indicate that the scaling of the relationship is well respected.The two scaling tendencies observed for the wind speed v correspond to the inertial range (small scales) and the energy-containing range (large scales).
We then estimate the scaling exponent functions ζ vH (q) and ζ vL (q) respectively for the inertial and the mesoscale ranges.The exponent functions ζ vH are estimated with 412 time series of length 7000 values sampled at 20 Hz and ζ vL with 5 time series of length 5 millions of values sampled at 1 Hz.The scaling exponent function ζ v (q) is the ensemble average obtained from all these scaling exponent functions ζ i (q):

Scaling exponent function ζ vH for the inertial range
Figure 5 presents the empirical scaling exponent function ζ vH , for the inertial range, compared to linear model K41 (ζ (q) = q/3) (Kolmogorov, 1941), and the lognormal model proposed by Kolmogorov (1962) that corresponds to the following relationship: where the parameter µ is called "intermittency parameter" and is close to 0.25 (Arneodo et al., 1996).The function ζ vH (q) obtained from our database is nonlinear and concave similar to the results obtained in previous studies (Schmitt et al., 1993;Schertzer et al., 1997), with the following properties: -0 < q < 3, ζ vH (q) > q/3 for small fluctuations, q > 3, ζ vH (q) < q/3 for large fluctuations.
The lognormal model provides a reasonable fit for the empirical exponent function ζ vH (q).The average of the fluctuations correspond to q = 1, and H = ζ vH (1) 0.36 is the socalled "Hurst" exponent characterizing the scaling nonconservation of the mean.The second moment ζ vH (2) = 0.68 is linked to the spectral exponent β 1 + ζ vH (2) = 1.68.

Scaling exponent function ζ vL for the energy-containing range
Figure 5 illustrates the empirical scaling exponent function ζ vL (q) obtained for the mesoscale range, compared to a linear model (ζ vL (1) 1/6) and a quadratic model ζ (q) = −0.023q 2 + 0.2q.The exponent function ζ vL is nonlinear and concave, characterizing a multifractal process.We obtain the value of H = ζ vL (1) 0.16 and the value ζ vL (2) 0.30, and consequently β 1.30 compatible with the value estimated by the Fourier analysis (β 1.27).The moment function K(q) = qH − ζ (q) is often used for characterizing a type of scaling.Lauren et al. (1999) performed a multiscaling analysis for characterizing the low-wavenumber statistical properties of surface layer winds.In this study, the data were measured with a cup anemometer at a sampling rate of 1 Hz.Then, the function K(q) was plotted for q ranging 0-3; for example K(1) and K(3) were estimated respectively to be 0 and 0.16.In our database the values of K(1) and K(3) were estimated respectively to be 0 and 0.13.

The structure functions scaling exponent ζ P (q) for the aggregate power output
Here we apply the structure function analysis to aggregate power output data from a cluster of wind turbine generators, for characterizing the type (monoscaling or multiscaling) and regimes of scaling for these data.We consider the increments of output power P ( t) = |P (t + t)−P (t)| for 1 t 10 6 s and for all moments between 0.25 and 5 with increments of 0.25.Figure 4b illustrates the scaling for the aggregate power output P for q = 1, 1.5, 2 and 2.5, in a log-log diagram.The structure functions display two scaling tendencies: for 1 t 6.5 × 10 4 s and for 6.5 × 10 4 t 10 6 s.Here, the straight lines black show also that the scaling of Eq. ( 2) is well verified.Then, we estimate the scaling ex- The scaling of the structure functions for wind power S P (q) illustrating two scaling regimes: for 1 t 6.5 × 10 4 s and for 6.5×10 4 t 10 6 s.(b) The empirical scaling exponent function ζ P H (q) ( ) estimated for 1 t 6.5 × 10 4 s and compared to the K41 model (black solid line) and the lognormal model (red dashed line).The empirical scaling exponent function ζ P L (q) ( • ) compared to a quadratic fit (blue dashed line).
ponent functions ζ P H (q) and ζ P L (q) from the slopes of the straight lines using a linear least regression, for each scaling regime.The scaling exponent function ζ P H (q) is estimated with i = 20 time series of length 10 6 values sampled at 1 Hz and ζ P L (q) with i = 2 time series of length 10 7 values sampled at 1H z.The scaling exponent functions ζ P (q) are the ensemble average obtained from all these exponents ζ i (q).
The empirical curves ζ P H (q) and ζ P L (q) for the aggregate power output P of this wind farm are illustrated in Fig. 6. ζ P H (q) is compared to the linear model K41 (q/3) and the lognormal model given in Eq. ( 6).For the timescales 1 t 6.5 × 10 4 s (≈ 18 h), ζ P H (q) is concave and nonlinear, characterizing a multifractal process.Moreover, the exponent function ζ P H is close to the lognormal model.ζ P H (q) presents the following properties similar to ζ vH (q): -0 < q < 3, ζ P H (q) > q/3 for small fluctuations, q = 3, ζ P H (3) = 1, q > 3, ζ P H (q) < q/3 for large fluctuations.
For the timescales t 6.5 × 10 4 s (∼ 18 h), ζ P L (q) is nonlinear and concave, characterizing a multifractal process.ζ P L (q) is fitted by a quadratic model.However, a representative estimate of ζ P L (q) must be performed with more data.
In short, the atmospheric wind data present two scaling regimes observed with spectral and structure function analysis.The scaling exponent functions ζ vH and ζ vL are fitted by an equation of type −Aq 2 + Bq with A characterizing the curvature indicating the degree of intermittency.Consequently, the turbulent wind speed is less intermittent in the inertial range (A = 0.013) than the energy-containing range (A = 0.023).The aggregate power output presents also two scaling regimes.For 1 t 6.5 × 10 4 s, the wind power is a multifractal process having the same properties as the wind speed in the inertial range, and for t 6.5 × 10 4 s the wind power is also a multifractal process less intermittent than the turbulent wind atmospheric in the mesoscale range.
6 Relation between the wind speed and the aggregate power output fluctuations

A scaling test for dependence: the coherence and the cross-correlation functions
The coherence function is often used in the wind energy community (Burton et al., 2001;Sørensen et al., 2002;Nanahara et al., 2004).In this work, the coherence function H vP between the measured wind speed v and the aggregate power output P from the wind farm considered is estimated from the co-spectrum E vp (f ) for the simultaneous data and density power spectrum E v (f ) and E p (f ) respectively for v and P .Figure 7a illustrates the coherence function H vP defined in Eq. ( 9), plotted for a frequency range 10 −7 f 1 Hz.Three regimes are observed: for f 10 −2 Hz corresponding to timescales T 100 s, H vP shows clearly the absence of correlation.For 10 −3 f 10 −2 Hz, corresponding to timescales 100 T 1000 s, the value of the coherence function increases from 0 to 0.7, highlighting the presence of a correlation between the wind speed v and the aggregate power output P for frequency range f 10 −3 Hz. Figure 7b illustrates the maximum of the cross-correlation function for a temporal range 1 T 10 7 s in a log-log plot.
We recall here the expression of scalar M vP computed for a given timescale T : where v T and P T are respectively the moving averages of wind speed and the power output data computer over a time period T , and σ v and σ P are the corresponding standard deviations.Figure 7b illustrates a the maximum M vP of the cross-correlation C vP in a log-log plot for the simultaneous wind speed/power output data.In the inset, we give examples of cross-correlation functions C vP for the simultaneous wind speed/power output data moving averages, computed for timescales T = 10 2 s, 10 3 , 10 4 , and 10 5 s.We can see that the scalar M vP is increasing for small timescales and reaches a plateau for large timescales with a transition around 10 2 −10 3 s; it shows the presence of a correlation between the wind speed and output power for timescales T 10 3 s.

A multiscaling test for dependence: generalized correlation functions (GCFs) and exponents (GCEs)
In the multiscaling framework, two stochastic processes, defined previously, v( t) and P ( t), characterized by their hth-and gth-order structure functions, are considered.The correlation between the two processes ( v( t)) h and ( P ( t)) g then becomes a function of the scale and of the statistical orders of moment h and g, expressed by the GCF.
Figure 8a illustrates the GCF c(h, g) plotted in log-log versus the time increment t 10 3 s, for simultaneous wind speed and output power measurements.The function c(h, g) represented is estimated for a constant value of the statistical order of moment g of wind power fluctuations (g = 3), and different values of statistical order of moment h of wind speed fluctuations (i.e.h = 1, 2, 3, from bottom to top).The function c(h, g) displays a scaling for timescales 10 3 < t 2.10 6 s.The GFEs r(h, g) are estimated from the linear regression slope of c(h, g), illustrated in Fig. 8b.We can observe that the correlation between the wind speed and the wind power fluctuations increases with increasing values of the statistical order.The iso-values of r(h, g) follow approximately a function of the form σ hg.Moreover, because of the symmetry of r(h, g) in the h-g plane, we expect to obtain a relation of proportionality between v( t) and P ( t). Figure 9 shows S(h, 0) versus S(0, g) indicating a proportional relation between v( t) and P ( t) of the form To illustrate this purpose, Fig. 10 presents the joint distribution P( v, P ) for wind speed fluctuations v and output power fluctuations P , for t = 10 4 s ≈ 3 h (a) and for t = 10 5 s ≈ 1 day (b), with the skeleton φ (red dashed line), indicating that a linear relation can be considered between the wind speed and the aggregate power output fluctuations at these timescales.The skeleton φ obtained with the joint pdf P( v, P ) is defined by φ = P max ( v, P ) = max P (P( v, P ).
For each time increment t 10 4 s, we estimate the coefficients a and b in the relation P ( t) = a v( t) + b.We find that b is always very small and that a has a small  shows the cross-correlations for the simultaneous wind speed/power output moving averages, computed for timescales T = 10 2 , 10 3 , 10 4 , and 10 5 s.t-dependence: a = A t B , with A = 10 2.2 and B = 0.04.Thus, We obtain the following expression to model the skeleton φ for t 10 4 s: This is not incompatible with the fact that the power P typically scales cubically with the wind speed v of the form P = cv 3 .This relation implies that dP = 3cv 2 dv.Hence to a first order, since v is larger, power increments P are proportional to wind speed increments v.

Discussions
The goal of this study was (ii) to characterize the coupling between the atmospheric wind speed and the aggregate power output measured simultaneously, in the multifractal framework.

The atmospheric wind speed
We have collected the wind speed data with a sampling rate of 20 and 1 Hz.The spectral analysis showed that the Fourier spectra E v (f ) follows a power law: E(f ) ∼ f −β for the small-scale and the mesoscale range.Thus, for the small timescales (0.05 T 1 s) corresponding to the inertial range, the Fourier spectrum was close to the Kolmogorov spectrum, with β = 1.68, in agreement with previous studies (Schmitt et al., 1993;Katul et al., 1995;Schmitt, 2007;Böttcher et al., 2007;Morales et al., 2011;Calif and Schmitt, 2012).The structure function analysis highlighted the concavity and the nonlinearity of the exponent function ζ vH (q).Furthermore, the theoretical quadratic relation for lognormal multifractals was well fitted (Böttcher et al., 2007;Morales et al., 2011;Schertzer et al., 1997).In Calif and Schmitt (2012), a lognormal continuous stochastic equation (Schmitt, 2003;Huang et al., 2008) is considered for modeling the atmospheric wind speed at these small temporal scales.
For the large timescales (1 T 10 7 s) corresponding to the mesoscale range, the Fourier analysis of the wind speed data sampled at 20 and 1 Hz showed a power law for the power spectral density E v (f ) ∼ f −β with β = −1.27.Contrary to the inertial range, the turbulence effects in the mesoscale range did not possess the same degree of universality.More recent papers on longitudinal low-wave-number spectra have provided evidence of a −1 power law in the atmospheric surface layer for neutral conditions, and a −5/3 power law for near-convective conditions (Katul et al., 1995;Katul and Chu, 1998;Kader and Yaglom, 1991).Several theoretical and phenomenological approaches have been devoted to the explanation of such −1 power law scaling at the low frequencies.Tchen (1953) proposed a theoretical analysis of an approximate spectral budget, based on interaction between the mean flow vorticity and the fluctuating vorticity.To model the energy production and energy transfer, phemenological analogies between molecular dissipation of turbulent kinetic energy, energy extraction from large to small scales and energy removal from the mean to turbulent flows have been developed using Heisenberg's turbulent viscosity approach (Heisenberg, 1948).Tchen (1953) established that (i) a spectrum close to −1 power law strong is likely to occur close to a rough surface because of a strong interaction between the mean flow vorticity and the fluctuating vorticity, and (ii) a spectrum close to a −5/3 power law is likely to occur far away from the rough surface because of a weak interaction between the mean flow vorticity and the fluctuating vorticity.This was verified experimentally in Klebanoff (1954) and Katul and Chu (1998).Recently Katul et al. (2012) also proposed a phenomenological theory based on Heisenberg's eddy viscosity approach for the explanation of the existence of −1 power law scaling.
The structure functions analysis for the wind speed indicated also scaling properties for 1 T 10 7 s.The exponent function ζ vL (q) was nonlinear and concave, and a quadratic relation was proposed for fitting ζ vL (q).However the exponent function ζ vL (q) was less concave than ζ vH (q): this indicated that the wind speed is more intermittent in the mesoscale range than the inertial range, as also shown in a previous study (Lauren et al., 1999).

The aggregate power output
We have also collected the aggregate output power from a wind farm with a sampling rate of 1 Hz.The spectral analysis showed two scaling regimes.For 2 T 9×10 3 s the slope was close to the spectral slope −5/3 as in a previous study (Apt, 2007) For timescales T 9 × 10 3 s, the Fourier analysis showed that the power spectral density E P (f ) of the aggregate power output from a wind farm has a spectral slope close to 1.2.The structure function analysis confirmed the obtained scaling, indicating the concavity and nonlinearity of the scaling exponent function ζ P L (q). ζ P L (q) was less concave than ζ P H (q), meaning a decrease of the intermittency with these timescales, as found also for the turbulent wind speed.However, for a more representative result, it is necessary to have more aggregate power output data.

The coupling between the wind speed and the output wind power fluctuations
We have collected simultaneously the wind speed v and the aggregate power output P , sampled at 1 Hz in order to highlight a possible dependence between both series.Generally, a scaling test, i.e., the coherence, is estimated in the wind energy community, for testing the dependence between the data.In this study, the coherence function H vP and the scalar M vP showed a coupling between the turbulent wind speed and the aggregate output wind power for timescales T 10 3 s.This study involved simultaneous analysis of turbulence time series and instantaneous power data, considering the multiscale correlations between both series.For that, a multiscaling test for independence between two stochastic processes was used.The generalized correlation functions c(h, g) and exponents r(h, g) were estimated.For timescales 10 3 T 10 5 s, the function c(h, g) presented a scaling relating a coupling between both series in agreement with the coherence function H vP and the scalar M vP .Moreover, their GFE r(h, g) indicated that larger fluctuations have stronger dependence.In fact, slow fluctuations (in the range of tens of minutes and hours) are mainly due to meteorological dynamics and are highly correlated among nearby wind farms.
A smoothing effect is operated by the size of the wind farm for short timescales (or short spatial timescales).This indicates that large fluctuations in wind speed over very short timescales (wind gusts) do not give rise to similar fluctuations in the aggregate power output.There is inertia in the wind turbine generator system, corresponding to a nonlinear transfer function from turbulent wind speed to aggregate power output production.This could also be caused by the spatial repartition combined with the number of wind turbines: turbulence at different wind turbines is weakly correlated or uncorrelated (Sørensen et al., 2002).Fast fluctuations have smaller spatial correlations (Nichita et al., 2002;Petru and Thiringer, 2002) and wind gusts are smoothed in the aggregate power output from a wind farm.Furthermore, the symmetry of r(h, g) in the plane (h-g) is in favor of a proportional relation between the wind speed fluctuations v( t), and the output power fluctuations P ( t) seem to be of the form P ( t) = K v( t) for these timescales.

Conclusions
This work analyzed the intermittency and the scaling properties of the wind speed v ant the aggregate power output P time series, through structure function analysis, at all scales and all intensities.The aggregate power output P possessed intermittent, multifractal properties and two scaling regimes in relation with the multifractal properties of the input wind speed v. Furthermore, we showed that the generalized correlation functions and exponents are a suitable tool for testing a dependence between two simultaneous intermittent processes.This multiscaling test allowed characterizing also the coupling of two simultaneously sampled intermittent processes by an analytical expression (power law, proportionality).Indeed, we proposed a proportionality expression for the coupling between the wind speed and the aggregate power output fluctuations.We believe that the actual relation is complex and has some nonlinear relations.We want here to understand the experimental result found about scaling relation for the power output.The relation given in relation ( 22) should be seen as a simplified model and not as a firm result.
We found here typical scales of 10 3 -10 4 s for the aggregate power output fluctuations.The scale of 10 4 s 3 h was found for the change of slope in aggregate power output spectral regime, and the scale of 10 3 s 15 min was found for the correlation between both series.This scale may be due to the local topography, the local meteorology, or the influence of the spatial disposition of the different wind turbine generators composing the wind farm.A better understanding of theses scales of 10 3 -10 4 s will be the subject of future works in order to consider a possible universality, or its relation with local conditions or wind farm composition.

Fig. 1 .
Fig. 1.Two examples of simultaneous wind power data: (a) simultaneous time series sampled at 1 Hz; (b) moving averages of simultaneous time series, over a time period T = 10 4 s 3 h.

Fig. 2 .
Fig. 2.An illustration of the Kolmogorov quasi-equilibrium energetic cascade in the inertial range.

Fig. 3 .
Fig. 3.The Fourier power spectra illustrating two scaling regimes: (a) averaged spectrum E v (f ) for the atmospheric wind speed sampled at 20 Hz displaying a −1.69 scaling for 1 f 10 Hz and a −1.28 scaling for f 1 Hz.The inset shows that the Fourier spectra of wind atmospheric sampled at 1 Hz display a −1.27 scaling over seven decades.(b) Averaged spectrum E p (f ) for the aggregate power output, for which a crossover between a −5/3 slope and −1.27 slope is also found.

Fig. 4 .
Fig. 4. The scaling of the structure functions of the wind speed S v (q) (a) for data sampled at 20 Hz and (b) for data sampled at 1 Hz.The structure functions display also two scaling regimes for 0.05 t 10 3 s for data sampled at 20 Hz and a scaling regime for 1 t 5 × 10 5 s.

Fig. 5 .
Fig. 5.The empirical scaling exponent function ζ vH ( ) for the inertial range compared to the linear model K41 (solid line) corresponding to homogeneous turbulence, and the lognormal model (dashed line) taking into account the intermittency effects.The empirical scaling exponent function ζ vL ( • ) for the mesoscale range, compared to a linear function ζ vL (1)q q/6 (blue solid line) and a quadratic function −0.023q 2 + 0.2q (green dashed line).
Fig. 6. (a)The scaling of the structure functions for wind power S P (q) illustrating two scaling regimes: for 1 t 6.5 × 10 4 s and for 6.5×10 4 t 10 6 s.(b) The empirical scaling exponent function ζ P H (q) ( ) estimated for 1 t 6.5 × 10 4 s and compared to the K41 model (black solid line) and the lognormal model (red dashed line).The empirical scaling exponent function ζ P L (q) ( • ) compared to a quadratic fit (blue dashed line).

Fig. 7 .
Fig. 7.The coherence function H vP (frequency domain) (a) and the maximum M vP (time domain) (b) of the cross-correlation function C vP in log-log plot for the simultaneous wind speed/power data: the insets in (a) shows the functions H vP in semilogx plot, and the inset in (b)shows the cross-correlations for the simultaneous wind speed/power output moving averages, computed for timescales T = 10 2 , 10 3 , 10 4 , and 10 5 s.

Fig. 8 .
Fig. 8. (a) The generalized correlation functions (GCFs) c(h, g) versus the time increments t in log-log plots for simultaneous wind speed and aggregate power output fluctuations.(b) The generalized correlation exponents (GCEs) r(h, g) estimated from the linear regression of c(h, g).

Fig. 9 .
Fig. 9. Plot of the generalized correlation exponents S(h, 0) versus S(0, g) This confirms the symmetry of the r(h, g) iso-values in the h-g plane.

Fig. 10 .
Fig. 10.The joint distribution P( v, P ) for wind speed fluctuations v and output power P , for (a) t = 10 4 s and (b) t = 10 5 s.

Schmitt: Multiscaling of the wind speed and wind power 381Table 1 .
Description of our database with the mean and the standard deviation values.
In contrast, if C xy and H xy are unity or close to unity, x and y are perfectly correlated; if C xy and H xy are close to −1, x and y are negatively correlated.
E xy (f ) is the Fourier co-spectrum, and E x (f ) and E y (f ) are the Fourier spectra of processes x(t) and y(t), respectively.If C xy and H xy are zero or close to zero, then x and y are uncorrelated.