- Email: [email protected]
Description and distribution fitting of transformed sodar wind observations
ARTICLE IN PRESS
Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100 www.elsevier.com/locate/jastp
Description and distribution fitting of transformed sodar wind observations Isidro A. Pe´rez, M. Luisa Sa´nchez, M. A´ngeles Garcı´ a, Beatriz de Torre Department of Applied Physics, Faculty of Sciences, University of Valladolid, c/Prado de la Magdalena, s/n 47071 Valladolid, Spain Received 11 July 2007; received in revised form 2 October 2007; accepted 17 October 2007 Available online 30 October 2007
Abstract Sodars are useful devices in atmospheric research beyond the first tenths of metres, the management of their observations proving equally important as the data itself. In this paper, a 3-year database of wind direction and speed was used, although in this case both variables have been transformed. The lowest measurement level was considered as a reference and differences in wind direction, ratios in wind speed and the power law exponent were calculated up to the 500 m level. First, the Ekman spiral was considered in order to fix its range as well as daily and seasonal evolutions. Diurnal cycles were then determined to establish the sharp contrast between day and night, and the difference in transitions in morning and afternoon. Histograms of wind direction differences, ratios in wind speed and the power law exponent were also calculated. Their shape was similar, evidencing only one mode with a rapid decrease from it, though the histogram was right skewed for transformed wind speed whereas it was symmetrical for the other two variables. These histograms were fitted by an original three-parameter distribution function using simple linear regressions, considering the least root mean square error as fitting criterion. The cumulative distribution function was also theoretically obtained. A close-up expression for the interquartile range as a function of the distribution function parameters was proposed with satisfactory results. Robust kurtosis behaviour was also observed. A bivariate distribution was calculated and satisfactorily contrasted with experimental values. Finally, the influence of wind speed on the distribution parameters was considered. r 2007 Elsevier Ltd. All rights reserved. Keywords: Sodar data; Wind statistics; Wind analysis
1. Introduction To date much is known about the lower atmosphere both in experimental as well as theoretical terms. However, experimental devices such as masts, and theoretical developments such as the similarity theory, are only valid up to a few tenths of Corresponding author. Tel.: +34 983 42 30 00x4187; fax: +34 983 42 30 13. E-mail address: [email protected] (I.A. Pe´rez).
1364-6826/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2007.10.004
metres, in the so-called surface layer. Beyond this layer, availability of observations is low and knowledge of its behaviour even scarcer. Sodars are devices that have been used in meteorological research over the last 30 years. At present, sodar observations are used in air pollution applications (Gariazzo et al., 2005; Venkatram et al., 2005) or meteorological studies (Martano et al., 2005; Puygrenier et al., 2005). Intense experimentation has yielded much robust evidence and the observations to emerge have proven useful
ARTICLE IN PRESS 90
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
to investigate the lower atmosphere beyond the limits of other devices. However, obtaining good observations is as important as the analysis of this information if satisfactory results are to be acquired (Pe´rez et al., 2006). The first step towards gathering information on a variable is by means of its histogram, the distribution function being its natural extension. This is an easy way of describing all the data. In this paper, a broad database was used to gain an insight into wind in the low atmosphere, which might prove useful in further wind energy and air pollution studies. Two variables are considered: wind direction, d, and wind speed, v. Wind direction analysis is infrequent due to its angular nature. Hence the interest in accurately describing it, since treatment of it must differ from that of linear variables (Farrugia and Micallef, 2006). By contrast, wind speed studies are required in wind energy applications. The behaviour of this variable is widely described by the Weibull and Rayleigh distributions (Bivona et al., 2003; Zhou et al., 2006). However, discrepancies between observations and the frequently used theoretical distributions have recently led to the use of different functions (Carta and Mentado, 2007; Jaramillo and Borja, 2004). Sodar data may be difficult to understand due to the random component of atmospheric variables, which is particularly important in wind. For this reason, the transformation of variables is a procedure that should be investigated. The main emphasis of this paper thus lies not in the analysis of wind direction and speed, but in related variables. The lowest level measured (40 m) has therefore been considered as a reference, due to the higher data availability. Wind direction has been used by differences, d– d40, and wind speed by ratios, v/v40. The power law exponent, a, has also been studied. The power law expression is common in wind speed applications (Al-Nassar et al., 2005), where the a value frequently used, 1/7, is a simple but serious restriction, requiring further detailed analysis since daily evolution has been observed (Rehman and AlAbbadi, 2005). Firstly, the Ekman spiral is obtained. This spiral is occasionally analysed in bodies of water, yet infrequently in the atmosphere. Daily cycles are then revisited with the new variables used. Histograms are then also considered. This paper shows that d–d40, v/v40 and a have similar distributions, which are described by a central maximum and frequencies decaying extremely quickly from it.
Consequently, treatment of d– d40, v/v40 and a may be simplified greatly, since their distributions may be considered by the same expression. An original three-parameter function is used, whose parameters are calculated and whose height variation is also explored. Some theoretical considerations concerning the three-parameter function are made to make this study independent from the measurement site. The cumulative distribution function is thus determined and calculated. With this function, the interquartile range (IQR) and the kurtosis are obtained in the range of the distribution function parameters. Finally, the usefulness of the distribution function proposed in this study to compact data is established, suggesting that this function might prove satisfactory if used in further wind speed analyses. 2. Experimental description The sodar (DSDPA.90-24, built by METEK GmbH) emits acoustic pulses of a certain frequency (2200 Hz) into the atmosphere and receives the backscattered signal, whose frequency is shifted according to the wind component parallel to the propagation of the acoustic waves (Doppler effect). Wind speed is obtained from this frequency shift. The equipment is installed in the Low Atmosphere Research Centre (CIBA), 339,464 m W, 4,630,866 m N UTM Zone 30 (WGS84), some 30 km NW of Valladolid (Spain). The location is a very extensive plateau 840 m above MSL, with no relief elements, thus ensuring horizontal homogeneity. Non-irrigated crops and grass make up the surrounding vegetation, the roughness length thus being only a few centimetres. The 3-year measuring period commenced on 1 August 2002. Data were continuously acquired over this period, the only noticeable interruptions occurring for about 25 days (10 days in May 2003 and 15 days in June 2003). The minimum height proposed for our measurements was 40 m and the maximum 500 m, measurements being limited to 20 m levels. Two kinds of files were generated by the device: the first comprising raw data, with values each 15 s, and the second considering 10-min averages provided simultaneously for the levels investigated. Only the latter files were used in this paper. The files yielded abundant information for each level measured, although only wind direction and speed were considered in this paper. Each value was accompanied by its plausibility code, which represents
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
the results of the plausibility test performed on the averaged power spectra, and was used to control data quality. Data were rejected according to standard criteria considered by the manufacturer. 3. Results 3.1. The Ekman spiral The first level (40 m) was taken as a reference due to its having the greatest amount of data available. The d– d40 and v/v40 relationships were calculated for each level every 10 min. d– d40 was chosen to make the change independent with the height of the wind direction from the air flow. v/v40 histograms showed a better behaviour for nocturnal data (only one distribution) than for v, which was a mixture of distributions. In order to investigate yearly and daily patterns, seasons were taken into account and, following atmospheric stability analyses, the day was considered as the period from 1 h after sunrise to 1 h before sunset. Moreover, three wind speed intervals were investigated, corresponding to v40, v45 and v410 m s1. Medians of both quantities were then calculated and are represented in Fig. 1. This method has the advantage that it provides the behaviour of variables against 10-min averages where the Ekman spiral is frequently obscured by direction shear. Several results were obtained:
v/v40 was higher during the night, reaching the highest value in winter, 2.5, with v40 m s1, whereas with higher wind speeds v/v40 was lower. d– d40 did not evidence as intense a contrast between day and night as v/v40. During autumn nights, it reached nearly 301. As height increased, a clockwise turning, which was more relevant in winter and autumn, appeared.
Diurnal profiles were virtually flat during spring and summer. However, a loop was observed during summer nights, which was attributed to low-level jets. The opposite case was found on winter nights, where the curve was less pronounced, depending on higher wind speed ratios. According to Fig. 1, two groups of seasons could be considered: winter–autumn and spring–summer. Finally, behaviour with height proved robust, since only a few outliers were present at the highest level with high wind speeds due to the low number of observations available.
91
3.2. Daily behaviour In Fig. 2, medians were calculated every 2 h at 100, 200 and 300 m, chosen due to the high number of observations available. The black point corresponds to midnight. The clockwise rotation of values revealed a loop, which was more prominent at 300 m than at 100 m, where wind speed and directional changes were lower. In winter and spring at 300 m, a noticeable contrast was observed between day, with low wind speed ratios, and night, with high wind speed ratios. In autumn, this contrast was lower. However, at 100 m, the loop was more pronounced in spring and summer. The transition between night and day was sharp, since intermediate values were not obtained. However, a soft evolution was observed from day to night. As in the analysis presented previously, summer exhibited a discrepant behaviour compared to the other seasons due to the lower contrast of wind speed between day and night. 3.3. Distribution fitting 3.3.1. Distribution law Histograms of d– d40 and v/v40 were calculated. Class width was 11 for d– d40 and 0.021 for v/v40. Fig. 3 represents only three levels in winter. d– d40, the direction shear of wind between levels, was close to zero for most observations. Hence, an additional advantage is that its angular nature may be ignored. Histograms were nearly symmetric, meaning that d– d40 treatment is simpler than for wind direction. Intensity of wind direction deviation when height increases from 40 m to the level considered may be established by the maximum of this distribution. Additionally, data spread increased with height. For v/v40, only one maximum was obtained during the night. The power law, z a v ¼ , (1) v40 40 was used to calculate a by considering wind speed measurements at only two levels, 40 m and any other level above it. a histograms are presented in Fig. 4 for winter at three levels where the width class is 0.02. An extremely sharp maximum was obtained since most observations lay in a narrow interval, mainly during the night. Figs. 3 and 4 reflect how the variables presented showed similar behaviour, the smooth shape suggesting that they might be described by a distribution function.
ARTICLE IN PRESS 92
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
Fig. 1. Seasonal, daily and wind speed variation for the Ekman spiral obtained by medians of d– d40 and v/v40.
In order to obtain a good fit for these histograms, functions such as f ðxÞ ¼
1 a þ bðx xmax Þ2
(2)
were considered, where f(x) is the number of observations in each histogram class and x is the
variable (d– d40, v/v40 or a). The shape of the distribution was selected to easily calculate the fit parameters by means of linear regressions by inverting Eq. (2). The second-order denominator was chosen since it is the lowest order that provides a finite maximum with a soft change at xmax. a controls the maximum value corresponding to xmax
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
93
Fig. 2. Seasonal evolution and height change of the daily cycle calculated with 2-h medians of d– d40 and v/v40. Black points correspond to midnight.
and b the fall from the maximum. The main emphasis of this paper lies in a and b, since xmax is merely an origin displacement, which was small for d– d40 and a. According to the analysis of directional variables from Anderson-Cook and Noble (2001),
xmax is used to centre the distribution function and to place the discontinuity away from the maximum. However, d– d40 was retained as the x-axis in the figures, since it is the variable directly related to observations.
ARTICLE IN PRESS 94
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
Fig. 3. d– d40 and v/v40 histograms in winter at three selected levels in grey, and their distribution function fitted from the least RMSE, as a black line (class width: 11 for d– d40 and 0.021 for v/v40).
a is a useful coefficient to transform wind speed between different heights, a fixed value normally being considered. This is assumed for practical reasons, although an interval for a values was observed. Negative values were noticeable during the day, which is a consequence of nearly flat profiles. By contrast, positive values prevailed during the night and a change with height was also observed. The following fitting method was used:
As a previous step, a variable number of extreme observations was discarded, since they might distort further linear fittings at the inversion of Eq. (2). Therefore, an xmax variable was proposed, considering the whole range of x values.
a and b parameters were calculated by linear regression for each xmax proposed. Root mean square errors (RMSE) between experimental and theoretical distributions for the whole range of x values were calculated. The xmax, a and b values, which provide the lowest RMSE, were chosen as the result and the theoretical distribution functions were drawn as black lines in Figs. 3 and 4. s(a)/a and s(b)/b, which were, in general, about 101 or even lower, were calculated for d– d40, v/v40 and a in order to quantify the relative dispersion of the two parameters.
The goodness of the fit, which may be described by the RMSE or the correlation coefficient (r) proved quite satisfactory. Although plots for RMSE and r are not presented, RMSE values were about
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
95
between day and night and among seasons were observed. The b range was wider during the day. This behaviour was less intense for a. The trend was nearly the same for all levels, the greater values corresponding to higher levels.
Parameters for a histogram fits were also calculated satisfactorily and are presented in Fig. 6, where lower values correspond to lower levels and higher values to higher levels. Diurnal values were higher than the nocturnal values, mainly at higher levels and seasonal couples were observed with the a parameter: winter–spring and summer–autumn. 3.3.2. Some properties of the 1/(a+b(xxmax)2) function The distribution function may be easily integrated, yielding the corresponding cumulative distribution function ! rffiffiffi 1 b 1 ðx xmax Þ þ F ðx0 Þ, F ðxÞ ¼ pffiffiffiffiffi tg (3) a ab where F(x0) is a fixed term to ensure a zero value for F(x) at x ¼ x0, the lowest x in the range used. The distribution function provided by Eq. (2) is symmetrical from xmax for d– d40 and a. However, it is not symmetrical for v/v40. This fact determines the shape of the cumulative distribution function, which is presented in Fig. 7 for winter nights at 100 m. Fittings were quite satisfactory. Having obtained the cumulative functions, some useful statistics were investigated: the IQR and the robust kurtosis, RK. IQR is a robust measure of spread that corresponds to the interval where half of the central data lie (Wilks, 2006), and the robust kurtosis Fig. 4. a Histograms in winter at three selected levels in grey, and their distribution function fitted from the least RMSE, as a black line (class width: 0.021).
15–20 in the lower levels and 5 in the higher ones and r was, in general, greater than 0.6 and frequently above 0.8. Fig. 5 shows the fitting parameters for d– d40 and v/v40. Some additional information was obtained from this scatterplot:
A linear trend for both parameters with more scatter for b. Additionally, some differences
RK ¼
IQR 2ðD9 D1Þ
(4)
is a measure of the distribution flatness (Sachs, 1982). D1 and D9 correspond to the 10 and 90 percentiles, respectively. IQR was calculated by means of proposed a and b values distributed in the interval of their calculated values, in accordance with Table 1. For v/v40, the distribution is skewed with a lower limit. For this reason 20 xmax values were tried from 1 to 3. The IQR independence against xmax was the first important result.
ARTICLE IN PRESS 96
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
Fig. 5. Fitting parameters for d– d40 and v/v40.
Fig. 6. Fitting parameters for a.
Fig. 7. Cumulative distribution function for winter nights, at 100 m. Experimental in black and theoretical, Eq. (3), in grey.
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
IQR behaviour, presented in Fig. 8, was nearly the same, and a linear relationship between both parameters, for the variables was considered. Lower a corresponded to lower IQR, whereas lower b corresponded to higher IQR. These linear relationships were fitted by means of new linear expressions b ¼ a0 þ b0 a.
(5) 0
The new a parameter was close to zero, except for lower IQR of v/v40 and a, corresponding to higher slopes. For this reason, this parameter was not considered. No linear relationship was apparent between a0 and b0 . However, the b0 parameter did show a power relationship with IQR through 00
b0 ¼ a00 IQRb .
(6)
IQR may then be calculated from Eqs. (5) and (6) neglecting a0 by means of 1=b00 b . (7) IQR ¼ aa00 This expression is no doubt useful since it allows IQR to be calculated in the range considered for a and b, once a00 and b00 are known. These parameters are presented in Table 2. When IQRs were recalculated with Eq. (7), some discrepancies were
97
observed at higher values, hence the limits presented in Table 2. Linear regressions between ordered experimental and theoretical IQRs were calculated with a variable number of points and the limit was the experimental IQR whose corresponding theoretical IQR was lower than the value of the regression line increased by 20%. The goodness of the fit was determined by the correlation coefficient between both IQRs. RK was calculated, although Fig. 9 only presents the results for d– d40 and a, since the behaviour for v/v40 was not so linear due to the relatively even values of this variable in the domain considered. Linear fits according to Eq. (5) were considered. However, a0 and b0 showed a strong linear relationship, hence the method followed for the IQR was not used for the kurtosis.
3.3.3. Bivariate distribution d– d40 and v/v40 distributions may be multiplied to obtain a bivariate distribution. Additionally, the observed histogram was determined and compared with the bivariate distribution in Fig. 10 for winter nights, at 100 m. The number of observations was
Table 2 Parameters of Eq. (6) Table 1 Limits of Eq. (2) parameters used to calculate IQR and kurtosis from the cumulative distribution function Parameter
d– d40 (deg)
v/v40
a
a b
0.001–0.1 0.000001–0.001
0.001–0.1 0.01–1
0.001–0.1 0.01–1
Variable
a00
b00
Upper IQR
r
d-d40 (deg) v/v40 a
5.133 1.846 2.223
2.137 1.932 1.904
109.887 2.76 2.63
0.995 0.992 0.992
Range of application and correlation coefficient between experimental and theoretical IQRs.
Fig. 8. Linear variation of the IQR calculated from the theoretical cumulative distribution function, Eq. (3), by proposing a and b in the ranges presented in Table 1. Numbers on the lines correspond to IQR for d– d40 (deg), v/v40 and a.
ARTICLE IN PRESS 98
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
Fig. 9. Linear variation of the kurtosis calculated from the theoretical cumulative distribution function, Eq. (3), by proposing a and b in the ranges presented in Table 1. Numbers on the lines correspond to kurtosis for d– d40 and a.
low due to the narrow classes considered (11 for d– d40 and 0.021 for v/v40). The goodness of the fit was determined by means of the correlation coefficient, which was calculated for all levels and seasons. This was around 0.9 for the lower levels, where agreement proved satisfactory and decreased gradually as height increased, until reaching the lower values, around 0.25, at the higher levels, due to the lower data availability.
3.3.4. Influence of wind speed The previous study was performed for wind speed above zero. However, the influence of wind speed was analysed through the selection of observations with a variable wind speed threshold in 1 m s1 steps. In other words, only wind speed observations above the thresholds were considered. As just a slight daily and seasonal evolution was observed, only the results for winter nights are presented in Fig. 11 where the ranges of variables were chosen for better observation of the behaviour of variables due to their wide interval of variation. a and b behaviour for d– d40 and v/v40 was quite similar and a presented small differences. In general, changes with wind speed were observed above 5 m s1. The height influence was appreciated for a with d– d40 and v/v40 below 5 m s1.
Fig. 10. Theoretical (dashed line) and experimental (continuous line) bivariate distribution obtained for winter nights, at 100 m (class width: 11 for d– d40 and 0.021 for v/v40). Numbers on the curves correspond to the frequency of observations.
4. Conclusions Three years of sodar data were analysed. Differences between wind direction at one level and the
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
99
Fig. 11. Height variation of a and b parameters depending on wind speed. Numbers on the curves correspond to a (upper graphs) and b (lower graphs) for d– d40 (left), v/v40 (centre) and a (right).
lowest level (40 m) were established as satisfactory variables to describe wind direction. Ratios between wind speed at one level and wind speed at the lowest level proved sufficiently satisfactory to characterise wind speed. The Ekman spiral was obtained as a median of d– d40 and v/v40 and its seasonal evolution presented. The angular deviation reached around 301 and the wind speed ratio 2.5. Low-level jets were well defined during summer nights. Two-hour medians of d– d40 and v/v40 were also calculated and represented to obtain daily cycles where a pronounced contrast between day and night was appreciated with a soft transition from day to night and a sharp change from night to day. d– d40, v/v40 and a, the power law exponent, histograms were obtained. d– d40 and a histograms were similar, symmetrical, with a central maximum and a fast decrease from it. This behaviour allowed us to take d– d40 as a linear variable. v/v40 showed only one maximum during the night that considerably simplified analysis of this variable. These
histograms were originally fitted to a three-parameter distribution law with simple linear regressions by considering the least RMSE. The a and b parameters of this distribution were calculated and increased with height. Once the theoretical distribution law was proposed, a theoretical study of the IQR and the kurtosis by means of the cumulative distribution expression was performed. Values of the a and b distribution law parameters were proposed in their ranges observed and a linear behaviour of both parameters according to IQR was obtained, enabling an empirical relationship between a, b and IQR to be established. For the kurtosis, only the linear behaviour between both parameters for d– d40 and a was observed. A bivariate distribution was obtained by a product of d– d40 and v/v40 distributions with satisfactory agreement with the experimental distribution. Finally, the height variation of the a and b distribution parameters with wind speed has been
ARTICLE IN PRESS 100
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
determined, the main changes being observed above 5 m s1. Acknowledgement The authors wish to acknowledge the financial support of the Interministerial Commission of Science and Technology and the Regional Government of Castile and Leon. References Al-Nassar, W., Alhajraf, S., Al-Enizi, A., Al-Awadhi, L., 2005. Potential wind power generation in the State of Kuwait. Renewable Energy 30, 2149–2161. Anderson-Cook, C.M., Noble, R.B., 2001. An alternate model for cylindrical data. Nonlinear Analysis 47, 2011–2022. Bivona, S., Burlon, R., Leone, C., 2003. Hourly wind speed analysis in Sicily. Renewable Energy 28, 1371–1385. Carta, J.A., Mentado, D., 2007. A continuous bivariate model for wind power density and wind turbine energy output estimations. Energy Conversion and Management 48, 420–432. Farrugia, P.S., Micallef, A., 2006. Comparative analysis of estimators for wind direction standard deviation. Meteorological Applications 13, 29–41. Gariazzo, C., Pelliccioni, A., di Filippo, P., Sallusti, F., Cecinato, A., 2005. Monitoring and analysis of volatile organic
compounds around an oil refinery. Water, Air and Soil Pollution 167, 17–38. Jaramillo, O.A., Borja, M.A., 2004. Wind speed analysis in LaVentosa, Mexico: a bimodal probability distribution case. Renewable Energy 29, 1613–1630. Martano, P., Cava, D., Mastrantonio, G., Argentini, S., Viola, A., 2005. Sodar detected top-down convection in a nocturnal cloud-topped boundary layer: a case study. Boundary-Layer Meteorology 115, 85–103. Pe´rez, I.A., Garcı´ a, M.A., Sa´nchez, M.L., de Torre, B., 2006. Fit of wind speed and temperature profiles in the low atmosphere from RASS sodar data. Journal of Atmospheric and SolarTerrestrial Physics 68, 1125–1135. Puygrenier, V., Lohou, F., Campistron, B., Saı¨ d, F., Pigeon, G., Be´nech, B., Serc- a, D., 2005. Investigation on the fine structure of sea-breeze during ESCOMPTE experiment. Atmospheric Research 74, 329–353. Rehman, S., Al-Abbadi, N.M., 2005. Wind shear coefficients and their effect on energy production. Energy Conversion and Management 46, 2578–2591. Sachs, L., 1982. Applied Statistics, A Handbook of Techniques. Springer, New York. Venkatram, A., Isakov, V., Pankratz, D., Yuan, J., 2005. Relating plume spread to meteorology in urban areas. Atmospheric Environment 39, 371–380. Wilks, D., 2006. Statistical Methods in the Atmospheric Sciences, second ed. Academic Press, Boston. Zhou, W., Yang, H., Fang, Z., 2006. Wind power potential and characteristic analysis of the Pearl River Delta region, China. Renewable Energy 31, 739–753.
Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100 www.elsevier.com/locate/jastp
Description and distribution fitting of transformed sodar wind observations Isidro A. Pe´rez, M. Luisa Sa´nchez, M. A´ngeles Garcı´ a, Beatriz de Torre Department of Applied Physics, Faculty of Sciences, University of Valladolid, c/Prado de la Magdalena, s/n 47071 Valladolid, Spain Received 11 July 2007; received in revised form 2 October 2007; accepted 17 October 2007 Available online 30 October 2007
Abstract Sodars are useful devices in atmospheric research beyond the first tenths of metres, the management of their observations proving equally important as the data itself. In this paper, a 3-year database of wind direction and speed was used, although in this case both variables have been transformed. The lowest measurement level was considered as a reference and differences in wind direction, ratios in wind speed and the power law exponent were calculated up to the 500 m level. First, the Ekman spiral was considered in order to fix its range as well as daily and seasonal evolutions. Diurnal cycles were then determined to establish the sharp contrast between day and night, and the difference in transitions in morning and afternoon. Histograms of wind direction differences, ratios in wind speed and the power law exponent were also calculated. Their shape was similar, evidencing only one mode with a rapid decrease from it, though the histogram was right skewed for transformed wind speed whereas it was symmetrical for the other two variables. These histograms were fitted by an original three-parameter distribution function using simple linear regressions, considering the least root mean square error as fitting criterion. The cumulative distribution function was also theoretically obtained. A close-up expression for the interquartile range as a function of the distribution function parameters was proposed with satisfactory results. Robust kurtosis behaviour was also observed. A bivariate distribution was calculated and satisfactorily contrasted with experimental values. Finally, the influence of wind speed on the distribution parameters was considered. r 2007 Elsevier Ltd. All rights reserved. Keywords: Sodar data; Wind statistics; Wind analysis
1. Introduction To date much is known about the lower atmosphere both in experimental as well as theoretical terms. However, experimental devices such as masts, and theoretical developments such as the similarity theory, are only valid up to a few tenths of Corresponding author. Tel.: +34 983 42 30 00x4187; fax: +34 983 42 30 13. E-mail address: [email protected] (I.A. Pe´rez).
1364-6826/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2007.10.004
metres, in the so-called surface layer. Beyond this layer, availability of observations is low and knowledge of its behaviour even scarcer. Sodars are devices that have been used in meteorological research over the last 30 years. At present, sodar observations are used in air pollution applications (Gariazzo et al., 2005; Venkatram et al., 2005) or meteorological studies (Martano et al., 2005; Puygrenier et al., 2005). Intense experimentation has yielded much robust evidence and the observations to emerge have proven useful
ARTICLE IN PRESS 90
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
to investigate the lower atmosphere beyond the limits of other devices. However, obtaining good observations is as important as the analysis of this information if satisfactory results are to be acquired (Pe´rez et al., 2006). The first step towards gathering information on a variable is by means of its histogram, the distribution function being its natural extension. This is an easy way of describing all the data. In this paper, a broad database was used to gain an insight into wind in the low atmosphere, which might prove useful in further wind energy and air pollution studies. Two variables are considered: wind direction, d, and wind speed, v. Wind direction analysis is infrequent due to its angular nature. Hence the interest in accurately describing it, since treatment of it must differ from that of linear variables (Farrugia and Micallef, 2006). By contrast, wind speed studies are required in wind energy applications. The behaviour of this variable is widely described by the Weibull and Rayleigh distributions (Bivona et al., 2003; Zhou et al., 2006). However, discrepancies between observations and the frequently used theoretical distributions have recently led to the use of different functions (Carta and Mentado, 2007; Jaramillo and Borja, 2004). Sodar data may be difficult to understand due to the random component of atmospheric variables, which is particularly important in wind. For this reason, the transformation of variables is a procedure that should be investigated. The main emphasis of this paper thus lies not in the analysis of wind direction and speed, but in related variables. The lowest level measured (40 m) has therefore been considered as a reference, due to the higher data availability. Wind direction has been used by differences, d– d40, and wind speed by ratios, v/v40. The power law exponent, a, has also been studied. The power law expression is common in wind speed applications (Al-Nassar et al., 2005), where the a value frequently used, 1/7, is a simple but serious restriction, requiring further detailed analysis since daily evolution has been observed (Rehman and AlAbbadi, 2005). Firstly, the Ekman spiral is obtained. This spiral is occasionally analysed in bodies of water, yet infrequently in the atmosphere. Daily cycles are then revisited with the new variables used. Histograms are then also considered. This paper shows that d–d40, v/v40 and a have similar distributions, which are described by a central maximum and frequencies decaying extremely quickly from it.
Consequently, treatment of d– d40, v/v40 and a may be simplified greatly, since their distributions may be considered by the same expression. An original three-parameter function is used, whose parameters are calculated and whose height variation is also explored. Some theoretical considerations concerning the three-parameter function are made to make this study independent from the measurement site. The cumulative distribution function is thus determined and calculated. With this function, the interquartile range (IQR) and the kurtosis are obtained in the range of the distribution function parameters. Finally, the usefulness of the distribution function proposed in this study to compact data is established, suggesting that this function might prove satisfactory if used in further wind speed analyses. 2. Experimental description The sodar (DSDPA.90-24, built by METEK GmbH) emits acoustic pulses of a certain frequency (2200 Hz) into the atmosphere and receives the backscattered signal, whose frequency is shifted according to the wind component parallel to the propagation of the acoustic waves (Doppler effect). Wind speed is obtained from this frequency shift. The equipment is installed in the Low Atmosphere Research Centre (CIBA), 339,464 m W, 4,630,866 m N UTM Zone 30 (WGS84), some 30 km NW of Valladolid (Spain). The location is a very extensive plateau 840 m above MSL, with no relief elements, thus ensuring horizontal homogeneity. Non-irrigated crops and grass make up the surrounding vegetation, the roughness length thus being only a few centimetres. The 3-year measuring period commenced on 1 August 2002. Data were continuously acquired over this period, the only noticeable interruptions occurring for about 25 days (10 days in May 2003 and 15 days in June 2003). The minimum height proposed for our measurements was 40 m and the maximum 500 m, measurements being limited to 20 m levels. Two kinds of files were generated by the device: the first comprising raw data, with values each 15 s, and the second considering 10-min averages provided simultaneously for the levels investigated. Only the latter files were used in this paper. The files yielded abundant information for each level measured, although only wind direction and speed were considered in this paper. Each value was accompanied by its plausibility code, which represents
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
the results of the plausibility test performed on the averaged power spectra, and was used to control data quality. Data were rejected according to standard criteria considered by the manufacturer. 3. Results 3.1. The Ekman spiral The first level (40 m) was taken as a reference due to its having the greatest amount of data available. The d– d40 and v/v40 relationships were calculated for each level every 10 min. d– d40 was chosen to make the change independent with the height of the wind direction from the air flow. v/v40 histograms showed a better behaviour for nocturnal data (only one distribution) than for v, which was a mixture of distributions. In order to investigate yearly and daily patterns, seasons were taken into account and, following atmospheric stability analyses, the day was considered as the period from 1 h after sunrise to 1 h before sunset. Moreover, three wind speed intervals were investigated, corresponding to v40, v45 and v410 m s1. Medians of both quantities were then calculated and are represented in Fig. 1. This method has the advantage that it provides the behaviour of variables against 10-min averages where the Ekman spiral is frequently obscured by direction shear. Several results were obtained:
v/v40 was higher during the night, reaching the highest value in winter, 2.5, with v40 m s1, whereas with higher wind speeds v/v40 was lower. d– d40 did not evidence as intense a contrast between day and night as v/v40. During autumn nights, it reached nearly 301. As height increased, a clockwise turning, which was more relevant in winter and autumn, appeared.
Diurnal profiles were virtually flat during spring and summer. However, a loop was observed during summer nights, which was attributed to low-level jets. The opposite case was found on winter nights, where the curve was less pronounced, depending on higher wind speed ratios. According to Fig. 1, two groups of seasons could be considered: winter–autumn and spring–summer. Finally, behaviour with height proved robust, since only a few outliers were present at the highest level with high wind speeds due to the low number of observations available.
91
3.2. Daily behaviour In Fig. 2, medians were calculated every 2 h at 100, 200 and 300 m, chosen due to the high number of observations available. The black point corresponds to midnight. The clockwise rotation of values revealed a loop, which was more prominent at 300 m than at 100 m, where wind speed and directional changes were lower. In winter and spring at 300 m, a noticeable contrast was observed between day, with low wind speed ratios, and night, with high wind speed ratios. In autumn, this contrast was lower. However, at 100 m, the loop was more pronounced in spring and summer. The transition between night and day was sharp, since intermediate values were not obtained. However, a soft evolution was observed from day to night. As in the analysis presented previously, summer exhibited a discrepant behaviour compared to the other seasons due to the lower contrast of wind speed between day and night. 3.3. Distribution fitting 3.3.1. Distribution law Histograms of d– d40 and v/v40 were calculated. Class width was 11 for d– d40 and 0.021 for v/v40. Fig. 3 represents only three levels in winter. d– d40, the direction shear of wind between levels, was close to zero for most observations. Hence, an additional advantage is that its angular nature may be ignored. Histograms were nearly symmetric, meaning that d– d40 treatment is simpler than for wind direction. Intensity of wind direction deviation when height increases from 40 m to the level considered may be established by the maximum of this distribution. Additionally, data spread increased with height. For v/v40, only one maximum was obtained during the night. The power law, z a v ¼ , (1) v40 40 was used to calculate a by considering wind speed measurements at only two levels, 40 m and any other level above it. a histograms are presented in Fig. 4 for winter at three levels where the width class is 0.02. An extremely sharp maximum was obtained since most observations lay in a narrow interval, mainly during the night. Figs. 3 and 4 reflect how the variables presented showed similar behaviour, the smooth shape suggesting that they might be described by a distribution function.
ARTICLE IN PRESS 92
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
Fig. 1. Seasonal, daily and wind speed variation for the Ekman spiral obtained by medians of d– d40 and v/v40.
In order to obtain a good fit for these histograms, functions such as f ðxÞ ¼
1 a þ bðx xmax Þ2
(2)
were considered, where f(x) is the number of observations in each histogram class and x is the
variable (d– d40, v/v40 or a). The shape of the distribution was selected to easily calculate the fit parameters by means of linear regressions by inverting Eq. (2). The second-order denominator was chosen since it is the lowest order that provides a finite maximum with a soft change at xmax. a controls the maximum value corresponding to xmax
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
93
Fig. 2. Seasonal evolution and height change of the daily cycle calculated with 2-h medians of d– d40 and v/v40. Black points correspond to midnight.
and b the fall from the maximum. The main emphasis of this paper lies in a and b, since xmax is merely an origin displacement, which was small for d– d40 and a. According to the analysis of directional variables from Anderson-Cook and Noble (2001),
xmax is used to centre the distribution function and to place the discontinuity away from the maximum. However, d– d40 was retained as the x-axis in the figures, since it is the variable directly related to observations.
ARTICLE IN PRESS 94
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
Fig. 3. d– d40 and v/v40 histograms in winter at three selected levels in grey, and their distribution function fitted from the least RMSE, as a black line (class width: 11 for d– d40 and 0.021 for v/v40).
a is a useful coefficient to transform wind speed between different heights, a fixed value normally being considered. This is assumed for practical reasons, although an interval for a values was observed. Negative values were noticeable during the day, which is a consequence of nearly flat profiles. By contrast, positive values prevailed during the night and a change with height was also observed. The following fitting method was used:
As a previous step, a variable number of extreme observations was discarded, since they might distort further linear fittings at the inversion of Eq. (2). Therefore, an xmax variable was proposed, considering the whole range of x values.
a and b parameters were calculated by linear regression for each xmax proposed. Root mean square errors (RMSE) between experimental and theoretical distributions for the whole range of x values were calculated. The xmax, a and b values, which provide the lowest RMSE, were chosen as the result and the theoretical distribution functions were drawn as black lines in Figs. 3 and 4. s(a)/a and s(b)/b, which were, in general, about 101 or even lower, were calculated for d– d40, v/v40 and a in order to quantify the relative dispersion of the two parameters.
The goodness of the fit, which may be described by the RMSE or the correlation coefficient (r) proved quite satisfactory. Although plots for RMSE and r are not presented, RMSE values were about
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
95
between day and night and among seasons were observed. The b range was wider during the day. This behaviour was less intense for a. The trend was nearly the same for all levels, the greater values corresponding to higher levels.
Parameters for a histogram fits were also calculated satisfactorily and are presented in Fig. 6, where lower values correspond to lower levels and higher values to higher levels. Diurnal values were higher than the nocturnal values, mainly at higher levels and seasonal couples were observed with the a parameter: winter–spring and summer–autumn. 3.3.2. Some properties of the 1/(a+b(xxmax)2) function The distribution function may be easily integrated, yielding the corresponding cumulative distribution function ! rffiffiffi 1 b 1 ðx xmax Þ þ F ðx0 Þ, F ðxÞ ¼ pffiffiffiffiffi tg (3) a ab where F(x0) is a fixed term to ensure a zero value for F(x) at x ¼ x0, the lowest x in the range used. The distribution function provided by Eq. (2) is symmetrical from xmax for d– d40 and a. However, it is not symmetrical for v/v40. This fact determines the shape of the cumulative distribution function, which is presented in Fig. 7 for winter nights at 100 m. Fittings were quite satisfactory. Having obtained the cumulative functions, some useful statistics were investigated: the IQR and the robust kurtosis, RK. IQR is a robust measure of spread that corresponds to the interval where half of the central data lie (Wilks, 2006), and the robust kurtosis Fig. 4. a Histograms in winter at three selected levels in grey, and their distribution function fitted from the least RMSE, as a black line (class width: 0.021).
15–20 in the lower levels and 5 in the higher ones and r was, in general, greater than 0.6 and frequently above 0.8. Fig. 5 shows the fitting parameters for d– d40 and v/v40. Some additional information was obtained from this scatterplot:
A linear trend for both parameters with more scatter for b. Additionally, some differences
RK ¼
IQR 2ðD9 D1Þ
(4)
is a measure of the distribution flatness (Sachs, 1982). D1 and D9 correspond to the 10 and 90 percentiles, respectively. IQR was calculated by means of proposed a and b values distributed in the interval of their calculated values, in accordance with Table 1. For v/v40, the distribution is skewed with a lower limit. For this reason 20 xmax values were tried from 1 to 3. The IQR independence against xmax was the first important result.
ARTICLE IN PRESS 96
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
Fig. 5. Fitting parameters for d– d40 and v/v40.
Fig. 6. Fitting parameters for a.
Fig. 7. Cumulative distribution function for winter nights, at 100 m. Experimental in black and theoretical, Eq. (3), in grey.
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
IQR behaviour, presented in Fig. 8, was nearly the same, and a linear relationship between both parameters, for the variables was considered. Lower a corresponded to lower IQR, whereas lower b corresponded to higher IQR. These linear relationships were fitted by means of new linear expressions b ¼ a0 þ b0 a.
(5) 0
The new a parameter was close to zero, except for lower IQR of v/v40 and a, corresponding to higher slopes. For this reason, this parameter was not considered. No linear relationship was apparent between a0 and b0 . However, the b0 parameter did show a power relationship with IQR through 00
b0 ¼ a00 IQRb .
(6)
IQR may then be calculated from Eqs. (5) and (6) neglecting a0 by means of 1=b00 b . (7) IQR ¼ aa00 This expression is no doubt useful since it allows IQR to be calculated in the range considered for a and b, once a00 and b00 are known. These parameters are presented in Table 2. When IQRs were recalculated with Eq. (7), some discrepancies were
97
observed at higher values, hence the limits presented in Table 2. Linear regressions between ordered experimental and theoretical IQRs were calculated with a variable number of points and the limit was the experimental IQR whose corresponding theoretical IQR was lower than the value of the regression line increased by 20%. The goodness of the fit was determined by the correlation coefficient between both IQRs. RK was calculated, although Fig. 9 only presents the results for d– d40 and a, since the behaviour for v/v40 was not so linear due to the relatively even values of this variable in the domain considered. Linear fits according to Eq. (5) were considered. However, a0 and b0 showed a strong linear relationship, hence the method followed for the IQR was not used for the kurtosis.
3.3.3. Bivariate distribution d– d40 and v/v40 distributions may be multiplied to obtain a bivariate distribution. Additionally, the observed histogram was determined and compared with the bivariate distribution in Fig. 10 for winter nights, at 100 m. The number of observations was
Table 2 Parameters of Eq. (6) Table 1 Limits of Eq. (2) parameters used to calculate IQR and kurtosis from the cumulative distribution function Parameter
d– d40 (deg)
v/v40
a
a b
0.001–0.1 0.000001–0.001
0.001–0.1 0.01–1
0.001–0.1 0.01–1
Variable
a00
b00
Upper IQR
r
d-d40 (deg) v/v40 a
5.133 1.846 2.223
2.137 1.932 1.904
109.887 2.76 2.63
0.995 0.992 0.992
Range of application and correlation coefficient between experimental and theoretical IQRs.
Fig. 8. Linear variation of the IQR calculated from the theoretical cumulative distribution function, Eq. (3), by proposing a and b in the ranges presented in Table 1. Numbers on the lines correspond to IQR for d– d40 (deg), v/v40 and a.
ARTICLE IN PRESS 98
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
Fig. 9. Linear variation of the kurtosis calculated from the theoretical cumulative distribution function, Eq. (3), by proposing a and b in the ranges presented in Table 1. Numbers on the lines correspond to kurtosis for d– d40 and a.
low due to the narrow classes considered (11 for d– d40 and 0.021 for v/v40). The goodness of the fit was determined by means of the correlation coefficient, which was calculated for all levels and seasons. This was around 0.9 for the lower levels, where agreement proved satisfactory and decreased gradually as height increased, until reaching the lower values, around 0.25, at the higher levels, due to the lower data availability.
3.3.4. Influence of wind speed The previous study was performed for wind speed above zero. However, the influence of wind speed was analysed through the selection of observations with a variable wind speed threshold in 1 m s1 steps. In other words, only wind speed observations above the thresholds were considered. As just a slight daily and seasonal evolution was observed, only the results for winter nights are presented in Fig. 11 where the ranges of variables were chosen for better observation of the behaviour of variables due to their wide interval of variation. a and b behaviour for d– d40 and v/v40 was quite similar and a presented small differences. In general, changes with wind speed were observed above 5 m s1. The height influence was appreciated for a with d– d40 and v/v40 below 5 m s1.
Fig. 10. Theoretical (dashed line) and experimental (continuous line) bivariate distribution obtained for winter nights, at 100 m (class width: 11 for d– d40 and 0.021 for v/v40). Numbers on the curves correspond to the frequency of observations.
4. Conclusions Three years of sodar data were analysed. Differences between wind direction at one level and the
ARTICLE IN PRESS I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
99
Fig. 11. Height variation of a and b parameters depending on wind speed. Numbers on the curves correspond to a (upper graphs) and b (lower graphs) for d– d40 (left), v/v40 (centre) and a (right).
lowest level (40 m) were established as satisfactory variables to describe wind direction. Ratios between wind speed at one level and wind speed at the lowest level proved sufficiently satisfactory to characterise wind speed. The Ekman spiral was obtained as a median of d– d40 and v/v40 and its seasonal evolution presented. The angular deviation reached around 301 and the wind speed ratio 2.5. Low-level jets were well defined during summer nights. Two-hour medians of d– d40 and v/v40 were also calculated and represented to obtain daily cycles where a pronounced contrast between day and night was appreciated with a soft transition from day to night and a sharp change from night to day. d– d40, v/v40 and a, the power law exponent, histograms were obtained. d– d40 and a histograms were similar, symmetrical, with a central maximum and a fast decrease from it. This behaviour allowed us to take d– d40 as a linear variable. v/v40 showed only one maximum during the night that considerably simplified analysis of this variable. These
histograms were originally fitted to a three-parameter distribution law with simple linear regressions by considering the least RMSE. The a and b parameters of this distribution were calculated and increased with height. Once the theoretical distribution law was proposed, a theoretical study of the IQR and the kurtosis by means of the cumulative distribution expression was performed. Values of the a and b distribution law parameters were proposed in their ranges observed and a linear behaviour of both parameters according to IQR was obtained, enabling an empirical relationship between a, b and IQR to be established. For the kurtosis, only the linear behaviour between both parameters for d– d40 and a was observed. A bivariate distribution was obtained by a product of d– d40 and v/v40 distributions with satisfactory agreement with the experimental distribution. Finally, the height variation of the a and b distribution parameters with wind speed has been
ARTICLE IN PRESS 100
I.A. Pe´rez et al. / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 89–100
determined, the main changes being observed above 5 m s1. Acknowledgement The authors wish to acknowledge the financial support of the Interministerial Commission of Science and Technology and the Regional Government of Castile and Leon. References Al-Nassar, W., Alhajraf, S., Al-Enizi, A., Al-Awadhi, L., 2005. Potential wind power generation in the State of Kuwait. Renewable Energy 30, 2149–2161. Anderson-Cook, C.M., Noble, R.B., 2001. An alternate model for cylindrical data. Nonlinear Analysis 47, 2011–2022. Bivona, S., Burlon, R., Leone, C., 2003. Hourly wind speed analysis in Sicily. Renewable Energy 28, 1371–1385. Carta, J.A., Mentado, D., 2007. A continuous bivariate model for wind power density and wind turbine energy output estimations. Energy Conversion and Management 48, 420–432. Farrugia, P.S., Micallef, A., 2006. Comparative analysis of estimators for wind direction standard deviation. Meteorological Applications 13, 29–41. Gariazzo, C., Pelliccioni, A., di Filippo, P., Sallusti, F., Cecinato, A., 2005. Monitoring and analysis of volatile organic
compounds around an oil refinery. Water, Air and Soil Pollution 167, 17–38. Jaramillo, O.A., Borja, M.A., 2004. Wind speed analysis in LaVentosa, Mexico: a bimodal probability distribution case. Renewable Energy 29, 1613–1630. Martano, P., Cava, D., Mastrantonio, G., Argentini, S., Viola, A., 2005. Sodar detected top-down convection in a nocturnal cloud-topped boundary layer: a case study. Boundary-Layer Meteorology 115, 85–103. Pe´rez, I.A., Garcı´ a, M.A., Sa´nchez, M.L., de Torre, B., 2006. Fit of wind speed and temperature profiles in the low atmosphere from RASS sodar data. Journal of Atmospheric and SolarTerrestrial Physics 68, 1125–1135. Puygrenier, V., Lohou, F., Campistron, B., Saı¨ d, F., Pigeon, G., Be´nech, B., Serc- a, D., 2005. Investigation on the fine structure of sea-breeze during ESCOMPTE experiment. Atmospheric Research 74, 329–353. Rehman, S., Al-Abbadi, N.M., 2005. Wind shear coefficients and their effect on energy production. Energy Conversion and Management 46, 2578–2591. Sachs, L., 1982. Applied Statistics, A Handbook of Techniques. Springer, New York. Venkatram, A., Isakov, V., Pankratz, D., Yuan, J., 2005. Relating plume spread to meteorology in urban areas. Atmospheric Environment 39, 371–380. Wilks, D., 2006. Statistical Methods in the Atmospheric Sciences, second ed. Academic Press, Boston. Zhou, W., Yang, H., Fang, Z., 2006. Wind power potential and characteristic analysis of the Pearl River Delta region, China. Renewable Energy 31, 739–753.