Østerild: A natural laboratory for atmospheric turbulence

Understanding atmospheric turbulence is without a doubt one of the most complex subjects in meteorology. However, its behavior can be more easily investigated by analysis of high-quality measurements of the velocity and temperature ﬂuctuations. Here, we show that at the large test station of wind turbines in Østerild, northern Denmark, measurements from a 250-m lightning mast provide unique insights into the behavior of atmospheric turbulence within the range of heights where modern large wind turbines operate. We illustrate that for the predominant wes-terly winds at the site, the ﬂow can be assumed to be close to homogeneous. This allows the analysis of the behavior with atmospheric stability and height of the wind speed and direction and of different turbulence measures, i.a., velocity variances and covariances, as well as turbulence spectra and their characteristics. It is shown, e.g., that for a wide range of atmospheric stability conditions, turbulence parameters up to 241 m are directly related to mean ﬂow characteristics, which will aid in improving the simulation of turbulence in, i.a., aeroelastic simulations. These measurements can be used for the evaluation of a wide range of models describing the meteorological conditions and the atmospheric ﬂow and its turbulent behavior.


I. INTRODUCTION
Boundary-layer meteorology and, for our particular interest, wind-power meteorology rely on flow and meteorological models that can describe the behavior of atmospheric turbulence as accurate as possible. However, the complexity of the description of atmospheric turbulence hinders the development of models of different fidelities. The development and improvement of such models thus rely heavily on our understanding of turbulence from measurements.
However, high-quality, long-term, and accurate turbulence measurements extending beyond the first tens of meters from the ground are rarely found. Cup anemometers are robust instruments, but they can only measure the horizontal velocity accurately and, thus, only the variance of the horizontal velocity. Sonic anemometers, nowadays almost as robust as cup anemometers, are able to measure all velocity and temperature fluctuations, but due to their measurement configuration, they are affected by self-induced flow distortion. 1,2 Further, when deployed on masts, both cup and sonic anemometers suffer from structure-induced flow distortion. 3 Alternatives to these forms of in situ anemometry, which avoid issues of flow distortion, are remotesensing techniques. Generally, remote sensors are able to probe the atmosphere in much more detail (in terms of vertical spacing) and cover higher ranges than conventional meteorological masts, but their measurement availability is generally lower, and due to the characteristics of the probe volumes and scanning configuration, their turbulence estimates are normally biased when compared to those from sonic-and cup-anemometer measurements. 4,5 Although relatively recently established, the measurements performed at the large test station of wind turbines in Østerild, northern Denmark, have already been used for different purposes. Lars en et al. 6 used cup and sonic anemometer measurements to analyze the height dependence of the wind speed power spectrum from the synoptic scale to the spectral gap. Karagali et al. 7 examined the imprint of surface heterogeneity on the mean wind flow around Østerild by analyzing measurements from two scanning lidars performing horizontal wide scans at 50 and 200 m. The latter campaign is part of the experiments of the New European Wind Atlas. 8 Kelly 9 used sonic anemometer and lidar measurements performed in Østerild, prior to the construction of the test station, to investigate the relation between the turbulence length scale and the mean vertical wind shear. Hannesd ottir and Kelly 10 characterized ramplike wind-speed fluctuations using the sonic anemometer measurements.
Here, the purpose is to introduce the Østerild measurements to boundary-layer and wind-power meteorologists, in particular, those from the sonic anemometers that cover the first 241 m. This is mainly done by analysis of the behavior with the atmospheric stability and height of wind speed, direction, velocity variances, covariances, and turbulence spectra. Further, we show how turbulence spectral tensor parameters behave with the height and atmospheric stability and the relation of the turbulence length scale with standard wind-power measurements, which is useful to improve the description of the threedimensional turbulence structure of a given site.
Section II introduces the Østerild facility and the measurements deployed on the lightning meteorological mast. Data handling (selection, filtering, and flow correction) is described in Sec. III. The general wind conditions at the site are described in Sec. IV including yearly and diurnal distributions of atmospheric stability and the behavior of wind speed distributions with the height. Section V describes the land use characteristics of the area and provides estimates of the roughness length for different wind directions. Turbulence intensity and atmospheric stability conditions are examined in Secs. VI A and VI B. Section VI C focuses on the close-to-homogeneous sector at Østerild: the behavior of the vertical profiles of the mean wind and direction is illustrated together with that of the velocity variances and covariances. For the latter sector, the turbulence spectra and their characteristics with height and atmospheric stability are examined in Sec. VI C 3. The last section provides some perspectives for future work and conclusions.

II. THE ØSTERILD FACILITY A. Østerild test station
The Østerild test station is located in northern Jutland, Denmark, on a land portion between the waters from Limfjorden and the North Sea [see Fig. 1(a)]. The station mainly consists of seven turbine stands (two more are being added), each with a mast deployed about 2.5-4 rotor diameters west of the turbine reaching the correspondent hub height (see https://www.vindenergi.dtu.dk/english/test-centers/oesterild for details of the turbines). The turbines are aligned in the northsouth direction extending almost 6 km. Two lightning masts are deployed north and south of the row of turbines. The terrain is flat, as shown in Fig. 1(b), and heterogeneous, with a mix between crop and agricultural land, urban areas, and forest.

B. The Østerild lightning masts
The Østerild lightning masts are 250-m tall towers with a triangular lattice structure (see Fig. 2). They are equally equipped with instruments at different levels. Here, we are mostly interested in the measurements from the Risø P2546A cup anemometers at 10, 40, 70, 106, 140, 178, 210, and 244 m above the ground level (agl-all heights are here referred to agl unless otherwise stated), a reference wind vane at 40 m, and Metek USA-1 sonic anemometers at 7, 37, 103, 175, and 241 m. The difference between the south and the north mast is that for the first, sonic and cup anemometers are mounted on booms at 180 from the north and for the second at 0 . There are redundant cup anemometers at 106 and 178 m.

III. DATA HANDLING A. Data selection and filtering
The analysis of the velocity and turbulence profiles is based on the sonic anemometer measurements mainly and is performed on 10-min statistics. Therefore, sonic and cup anemometer measurements were extracted when all sonic anemometers deployed on the south lightning mast were available. We focus our analysis on this mast because (1) it has more availability than the north mast (99% vs 71%), (2) the frequency of south winds is higher than that of the north winds (see Sec. IV), and so the wakes from the wind turbines and direct mast shadow influence less often the measurements on the south than those on the north mast, and (3) there is a forest clearance west of the south mast; since the predominant winds at Østerild are westerlies, the flow, to a reasonable extent, can be considered homogeneous at this position. Therefore, we can eventually compare the observations with canonical flow solutions and simulations (see Sec. VI C for details). Our analysis concentrates on observations covering the period starting on April 1, 2015 to March 5, 2019. Journal of Renewable and Sustainable Energy ARTICLE scitation.org/journal/rse We first filtered data with regard to the minimum observed mean wind speed (hereafter, the term mean refers to the average over a 10-min period); the cup anemometer at 10 m must have recorded wind speeds higher than 3 m s À1 . 103 936 10-min (%24 months) were left after applying this criterion.
Second, we filtered data with regard to the quality of the sonic anemometer 20 Hz records within a 10-min period; for all records, the sonic anemometer status must have been either 0 or 1, i.e., that new data were acquired by the system without errors. Also, we only allowed complete 10-min time series for the analysis. This means that all sonic-derived statistics were based on 10-min time series containing 12 000 records. A total of 63 427 10-min (%15 months) were left for the analysis.

B. Correction and usage of sonic anemometer measurements
For each 10-min period, we removed the 2D correction and applied the 3D correction suggested by METEK GmbH, 11 which was shown to correct the velocity components in such a way that the ratio of the w-to the u-velocity spectrum and that of the v-to the u-velocity spectra within the inertial subrange become very close to 4/3 as expected from Kolmogorov's hypothesis. 2 We also applied the corrections suggested by Liu et al. 12 for the computation of the sonic heat flux, which we assume to be equal to the virtual potential kinematic heat flux. We applied yaw and pitch rotations to the sonic anemometer 10-min time series, so that the u-velocity component becomes aligned with the mean wind with the rotation angles determined each 10-min period.
We compute fluctuations of velocities, their variances and covariances, and other turbulent fluxes from the sonic anemometer measurements. The friction velocity is computed as where the prime ( 0 ) denotes fluctuations and the overbar a mean. The Obukhov length is computed as where T s is the sonic anemometer temperature, j is the von K arm an constant (we use a value of 0.4), and g is the earth's gravitational acceleration. We also compute the three velocity component autospectra and uw-cospectra. We classified the sonic anemometer observations within atmospheric stability classes by estimating the dimensionless stability parameter z=L O , where z is the height above ground, and using the ranges in Table I.

C. Mast and wind turbine shading
Since we concentrate the analysis on the sonic anemometer measurements, we want to avoid the direct mast shadow on the observations and the influence of the wind turbines' wakes in the stands. Turbine 7, which is the southern in the row, directly shadows the south mast for 32 winds, and thus, assuming a wide wake expansion of 15 , we guarantee a wake-free sector from 47 . Turbine 1 direct "shadow" (if any) is at 3 from the south mast, and so the wake-free sector is within the range of 47 -348 . Figure 3 shows the ratio of the horizontal wind speed measured by the sonic to the twin cup anemometer at two levels as a function of wind direction (similar results are found for the levels at 103 and 175 m). Three sets of data are shown: a group of 10-min periods within an  , where the sonic anemometer measurements show a systematic misbehavior, which is clearly illustrated as cup to sonic wind speed ratios are far from one, 10-min periods avoiding the direct influence of turbine wakes and mast shadow, and 10-min periods within the 270 6 15 sector with reference to the direction measured by the 37-m sonic anemometer. It is clearly illustrated that most outliers are filtered out by removing data from the eight-month period. For the wake-and shadow-free directions in Fig. 3(a), a few outliers are shown; these correspond to groups of consecutive 10-min periods, and so they might be caused by temporal nonreported data acquisition issues. Data left after filtering the eight-month period and the turbine/mast wake sectors amount for 51 209 10-min (%12 months), and those within the 270 6 15 sector are 14 355 10-min (%3 months). Figure 4 shows the wind roses from the measurements of two cup anemometers separated 204 m vertically and illustrates the predominance of westerlies at the Østerild site, being the winds from the northern directions, the less frequent. In this section, we analyzed the cup anemometer measurements only, and we do not follow the data filtering criteria described in Sec. III A, since we want to show all possible recorded winds for wind climatology. Figure 5(a) shows the histogram of wind speed from measurements of a number of cup anemometers together with the distribution of wind speed derived from computing the Weibull scale and shape parameters A w and k w using a maximum likelihood estimator on the time series. As illustrated, these all-sectors histograms closely follow Weibull distributions. A similar analysis was performed for the sector 270 6 15 (not shown), and the Weibull distributions follow even closer the histograms compared to the all-sector analysis. As expected A w increases with the height as it is directly related to the mean wind speed, whereas k w increases up to a maximum, the reversal height, and decreases slowly above that level. This behavior combines the effect of the diurnal and the synoptic variation. 13 Interestingly, for the all-sector dataset, k w values are slightly lower, and thus, the variability of 10-min mean winds is higher than that for the 270 6 15 sector, which is expected due to the heterogeneity of the area. Figure 6 illustrates the frequency of occurrence of the atmospheric stability classes in Table I

V. ROUGHNESS AND LAND-USE CHARACTERISTICS
The roughness length can be estimated using different methodologies (e.g., from surface characterization or by analyzing wind measurements). One can derive it also in several ways from observations. Figure 7(a) shows for 30 wind sectors the roughness rose derived using the 7-m sonic anemometer observations and that using the 37-m sonic anemometer observations within neutral atmospheric conditions assuming that the logarithmic wind profile can be applied locally, i.e., where j is the von K arm an constant (¼ 0.4) and U ¼ ðu 2 þ v 2 Þ 1=2 is the magnitude of the mean horizontal wind. For each sector, the value shown in Fig. 7(a) is exp h ln z 0 i, which represents an ensembleaverage roughness length.
The results for each of the sonic-anemometer observations show different patterns, except that in both cases, the 270 6 15 sector shows the lowest values. By using the measurements from the 7-m level, lower z 0 -values are generally found, except for the 120 and 150 sectors. The trees around the south lightning mast are 15-20-m tall, and so the roughness estimation for all but 270 6 15 is biased by the displacement height, which cannot be accurately estimated due to the surface heterogeneity. Figure 7(b) shows the distributions of z 0 computed from Eq. (3) for the 270 6 15 sector under neutral stability conditions. The values estimated from the 37-m level are more scattered than those at 7 m, as expected, because more nonlocal effects (e.g., from the surrounding forest and villages) affect the observations. They also show a higher Journal of Renewable and Sustainable Energy ARTICLE scitation.org/journal/rse ensemble-average value than those of the 7-m level, illustrating the influence of the high roughness elements on the 37-m level. Meteorological models, however, rely on land use datasets, and it is up to the model and the "user experience" to assign roughnesslength values to given land-use classes. Floors et al. 14 discussed a technique, which objectively produces roughness-length maps from treeheight maps that are created using airborne laser scans. The technique basically estimates the elevations of the point cloud within a spatial grid cell of typically 20 m Â 20 m. The tree height is estimated as the maximum of the difference between the elevation of all points within the grid cell and the median of those points. The roughness map around Østerild in Fig. 8 assumes that the roughness length is 10% the estimated tree height. 15 Figure 9 shows the behavior of the turbulence intensity (TI ¼ r u =U, where r u is the standard deviation of the u-velocity component) with wind speed for the different heights and for both the wake-and shadow-free and the 270 6 15 sectors. The range of TIs is large (%0.05-0.3), which shows the variety of turbulence conditions of the site. For both cases, the TI behavior with wind speed is very similar for the two highest levels and the main differences appear at the lowest three levels. For the latter height ranges, the TI is higher for the broader than for the narrow sector as the forest nearby strongly influences these heights.

VI. TURBULENCE CONDITIONS A. Turbulence intensity
Further, for both cases, TI behaves with wind speed as expected; it decreases with increasing wind and decreases with the height. For the three highest levels, the TI tends to stabilize (or even increase) for winds !12 m s À1 , which is a behavior typical of near-neutral conditions. 16 B. Atmospheric stability Figure 10 illustrates frequency distributions of the dimensionless stability for both the wake-and the shadow-free and the 270 6 15 sectors for the different heights. Although for the narrower sector, the frequencies are more evenly distributed than for the broader sector, both show the same pattern with the height; at 7 m, most conditions are close to neutral (z=L O % 0) and the distributions systematically spread, covering further stable and unstable conditions, the higher the level of the sonic anemometer. For the narrow sector, the ranges of dimensionless stabilities in Table I are also shown. Figure 11 shows, based on either the 37-or the 241-m sonic anemometer measurements, the occurrence of atmospheric stability classes per wind speed bin for the wake-and shadow-free sectors. Based on the 37-m measurements, the stability conditions are predominantly near-neutral for all wind speed ranges with the very unstable and very stable conditions within the lowest wind speed ranges. Based on the 241-m measurements, however, neutral conditions are generally less predominant within the whole wind speed range, although always present. Unstable and very unstable conditions are predominant within the low wind speed range, and stable and very stable conditions prevail within the high wind speed range.   As explained earlier, due to the surface characteristics within the 270 6 15 sector, we might have the possibility to study canonical flows at Østerild. To understand how homogeneous the atmospheric flow is within this range of directions, we can start by studying the behavior of the surface-layer dimensionless wind shear, with dimensionless stability z=L O . The dimensionless wind shear should follow surface-layer similarity 17 under homogeneous flow conditions, i.e., it should not deviate largely from the empirical expressions relating both dimensionless parameters. 18 Here, dU/dz was computed at 37 m by first fitting the vertical profile of wind speed U to the polynomial form of H€ ogstr€ om 19 from measurements of the three lowest sonic anemometers, and / m is thus referred to 37 m as we use the u Ã value at this level. Figure 12 shows the behavior of / m with z=L O from the sonic anemometer measurements together with the "modified" expressions of the Kansas experiment. 19 As illustrated, the observed dimensionless wind shear is in good agreement with the modified Kansas expressions for all stability ranges except that the latter overpredict the dimensionless wind shear for z=L O ! 0.5 as observed in other sites and experimental campaigns. 20,21 Since the behavior with the stability of the dimensionless wind shear closely follows surface-layer similarity, we can expect that by taking the ratio of the wind to the surfacelayer velocity scale (u Ã ), the vertical wind profile under neutral conditions follows the logarithmic wind profile within the surface layer and that both unstable and stable conditions deviate from it in accordance with surface-layer similarity. This is illustrated in Fig. 13 for all stability classes. At 7 m, the differences in the   Table I.

ARTICLE
scitation.org/journal/rse computed hU=u Ã37 i values between all stability classes become small, which reflects that the roughness length is close in all stability cases; the largest discrepancies are found in very unstable and very stable conditions, as expected, since for these two ranges, atmospheric stability has a large effect on the normalized wind speed already at 7 m, and these stability ranges occur at particular months during the year, where the roughness length might be different from the all-season estimation. In the plot, the logarithmic wind profile prediction is also illustrated using the roughness length value estimated in Sec. V for the 37-m height. The prediction fits well the observations, but the latter deviates the highest we observe above the ground. This is partly because the turbulence length scale does not increase linearly with the height but has a limiting value. 22 It can also be seen that the standard error is the highest within the very stable case, which is due to the higher wind variability under these atmospheric conditions. Similar to Fig. 13, we computed the turning of wind with the height for the number of stability conditions observed at the 37-m level and computed the relative direction change using the 7-m level as the reference (see Fig. 14). The behavior of the turning of the wind with the height seems to be consistent with what is expected; the relative direction change mostly increases with the height (with some exceptions for the unstable classes within the two first height levels),   Fig. 7(a).

ARTICLE
scitation.org/journal/rse with the change being the highest for very stable conditions (up to more than 18 at 241 m) and lowest for very unstable conditions (close to 3 at 241 m). Figure 15 shows the behavior with the height and stability of the three velocity variances and the uw-covariance normalized by the squared of the friction velocity at 37 m, u 2 Ã37 . As expected, the highest to lowest velocity variances are those of the u-, v-, and w-velocity components, respectively, followed by the uw-covariance. The more unstable the atmosphere, the closer the u-, v-, and w-variances become; for very unstable conditions, r 2 u % r 2 v % r 2 w , and for the two highest sonic levels as the atmospheric flow becomes close to isotropic. The normalized velocity variances under neutral stability conditions show the lowest values at several sonic levels.

Velocity variances and covariances
Generally, the behavior with the height of the variances follows the same pattern; without accounting for the 7-m level, whose values can be higher or lower than those at 37 m, suggesting that local terrain features highly influence the measurements at the lowest sonic anemometer level, the velocity variances and uw-covariance decrease with the height. The only exceptions are found for very stable conditions, in which all variances and uw-covariance peak at 103 m, and for unstable and very unstable conditions for the w-variance, where there is an increasing behavior up to 175 m, as it is also normally expected for such stability conditions. 23 When concentrating the analysis to neutral atmospheric conditions only and to the 37-m level, the relations r u =u Ã ¼ 2:26; r u =r v ¼ 1:29, and r u =r w ¼ 1:68 are found, which are close to those estimated over a less rough site within an homogeneous sector. 16 For all close-to-neutral conditions, the ratio hu 0 w 0 =u 2 Ã37 i becomes very close to À1 close to the surface, as expected from the very low contribution of the vw-covariance in u Ã [see Eq. (1)].

Spectral analysis
The three velocity component autospectra and uw-cospectra were also classified into the stability classes in Table I based on z=L O values from the 37-m sonic anemometer and then ensemble-averaged for each height. We fit the three-dimensional turbulence spectral model of Mann 24 (hereafter, the Mann model) to those spectra. This model is commonly used for describing atmospheric turbulence in a wide range of turbulence conditions 22,25 and for generating turbulence fields useful for, i.a., aeroelastic calculations. 26 The fitting procedure outputs the Mann model parameters, which, besides the wave vector k, describe the velocity spectral tensor: ae 2=3 , L, and C, with a being the spectral Kolmogorov constant, e the specific rate of destruction of turbulent kinetic energy, L a length scale related to the size of the turbulent eddies, and C a parameter describing the anisotropy of the turbulence. Figure 16 illustrates both the ensemble-averaged measured spectra and the Mann model fit for neutral conditions and for all sonic-anemometer levels. The same procedure was performed to the spectra within the other stability classes (not shown). The measured spectra behave as expected for atmospheric flow under flat and homogeneous land; the u-spectra peak the highest followed by the v-and the w-spectra, and the uw-cospectra are negative due to vertical wind shear. When moving away from the surface, the spectra peak generally lower (except for the 7-m level), indicating a decrease in the turbulence level, and the spectra peak at lower wavenumbers, indicating that the length scales are larger than those close to the surface. 22 As shown, the Mann model fits well the measured spectra for all heights, with the largest differences for the uw-cospectra for which the Mann model overpredicts its energy content. 27 Although not shown, similar results are found under the other stability classes at all levels. Figure 17 illustrates the behavior with the height and atmospheric stability class of the Mann parameters, which were fitted from the observed spectra. In agreement with the results by Peña et al. 22 for flow over a more homogeneous and flat area, we find that C generally decreases with the height and that the relative differences with the height are larger, the more unstable the atmosphere is, reaching values about unity at the highest levels. For all stable and near-neutral stabilities, we find C % 3.
Also, we find that ae 2=3 decreases with the height, as expected due to the decrease in turbulence intensity with the height being the highest the closest the atmospheric stability classes are to neutral, also as expected, as the close to neutral conditions show the highest velocity variances. Finally, L, the turbulence length scale, shows an increasing behavior with the height, as expected since the turbulent eddies also grow with the height, and the size of L is the highest the more unstable and the lowest the more stable the atmosphere is, also in very good agreement with the results obtained by Peña et al. 22 For completeness, the suggested values by IEC 28 for C and L, when designing wind turbines, are also shown. The measurements at Østerild show that turbulence is more isotropic and that the length scales are generally larger than the suggested values by the standard particularly the highest we observe above the ground.
An interesting question is whether or not we can estimate the Mann parameters from "traditional" measured flow parameters such as the velocity variances. Peña et al. 22 showed that the Mann-model length scale is related to the length scale of the wind profile, i.e., the Journal of Renewable and Sustainable Energy ARTICLE scitation.org/journal/rse mixing-length scale, which from mixing-length theory is the ratio of the friction velocity to the vertical wind shear. They showed that for a range of atmospheric stability conditions, wind speed ranges, and heights, L % 1:70u Ã =ðdU=dzÞ. Kelly 9 proposed the relation L % r u = ðdU=dzÞ, which can be computed directly from cup anemometer measurements. Figure 18 shows the behavior of L with the latter relation for a number of atmospheric stability conditions and heights. L behaves linearly with both ratios within all stable and close to neutral ranges. Under unstable and very unstable regimes at the two highest measurement levels, the behavior is not linear mainly due to the very low vertical wind shear.

VII. CONCLUSIONS AND PERSPECTIVES
The Østerild test station provides us with a high-quality, long-term, and accurate dataset of meteorological measurements. Such measurements are needed for the evaluation and improvement of meteorological and turbulence models. We showed that, in particular, the dataset from the southernmost lighting mast allows us to study the long-term wind climatology and its variation, which is required, e.g., in wind resource assessment. It was shown that the winds at Østerild experience all types of atmospheric stability conditions during the year, which is useful to understand the behavior of the mean wind and direction and turbulence within a wide range of flow conditions. Particularly, for the predominant winds, we argued that the flow can be considered homogeneous and showed that close to the surface, it follows Monin-Obukhov similarity. The velocity variances and uwcovariance showed the expected behavior with both stability and height. Finally, spectral analysis revealed the behavior of the velocity spectral tensor characteristics with both the height and stability, and it was further demonstrated that the turbulence length scale is linearly related to the variance of the wind and the vertical wind shear.
Evaluation of meteorological models, ranging from meso-to microscales, and microscale flow models of different fidelities, ranging from linearized Navier-Stokes equations to large-eddy simulations, is currently undergoing using the Østerild measurements. The particular predominant westerlies at the south lightning mast allow us the Journal of Renewable and Sustainable Energy ARTICLE scitation.org/journal/rse comparison of idealized simulations of canonical flows and help understanding the behavior of the models over close to homogeneous conditions. The heterogeneity of the area challenges the microscale model abilities. Also of particular interest is to evaluate the accuracy of meso-micro-coupling techniques at Østerild. Inputs to mesoscale models, such as land use information, face challenges for representing the heterogeneity of the area and show the need to properly downscale from the meso-to the microscales. Work in progress is related to determining what kind of inputs (e.g., land use, boundary, and initial conditions) both types of models should use to get the most accurate  coupling result. Remote-sensing techniques, i.a., wind and aerial lidars, and satellite information, are constantly being used at Østerild to characterize both the flow beyond the range of the mast and the surface heterogeneity. Journal of Renewable and Sustainable Energy ARTICLE scitation.org/journal/rse