2.1.1 Sensor response time and resolved variability

The pyranometer and pyrheliometer instruments are thermopiles, meaning that there is a nonzero response time to variations in incoming radiation, i.e., the time it takes for the thermopile to adjust to changes in irradiance signal. Thus, the true resolved resolution is not 1 Hz, although whether or not this impact is noticeable depends on the magnitude and rate of change in irradiance variability. According to the manufacturer's specifications, the CH1 pyrheliometer has a 7 s (95 %) or 10 s (99 %) response time (Kipp & Zonen, 2001), whereas these values are 1.66 s (66 %) or 5 s (95 %) for the CM22 pyranometers (Kipp & Zonen, 2004). This results in a likely underestimation of variability at 1 Hz, but the exact extent of this underestimation is hard to quantify given the absence of solar irradiance measurements with fast-responding sensors at a similar location as well as the lack of measurements that are long enough to sample diverse weather conditions. Figure 2 illustrates the power spectral density (PSD) of a year of BSRN 1 Hz data. Tabar et al. (2014) present spectra (their Fig. 2) for 1 year of data from a semiconductor pyranometer (with a >1 Hz response time) that show an order of magnitude higher PSD at 1 Hz compared with 0.1 Hz, as in our Fig. 2. However, these data were collected in Hawaii, which is a very different geographical location and climate compared with Cabauw. Alternatively, 2 weeks of summer time 10 Hz irradiance observations from fast radiometers (Mol and Heusinkveld, 2022) show a steeper decline between 0.1 and 1 Hz than the BSRN data. This supports the idea that the true cloud-driven irradiance variability does not always extend towards 1 Hz and higher frequencies and that the slow response time often has no noticeable effect.

Figure 2Power spectral density of 1 year (2016) of 1 Hz BSRN data (GHI component). The 1 min spectrum is based on resampled 1 Hz data. As a reference, $f^{-\mathrm{5}/\mathrm{3}}$ scaling is added to emphasize the steep decline after 0.1 Hz towards 1 Hz. An additional comparison is shown for the spectrum from a semiconductor radiometer deployed during the FESSTVaL (Field Experiment on submesoscale spatio-temporal variability in Lindenberg) campaign from 14 to 30 June 2021 near Berlin, Germany.

Download

Based on the technical specifications of the pyranometer and pyrheliometer, we are at least confident that variability is resolved up to 0.1 Hz (10 s). This is also supported by van Stratum et al. (2023), who showed agreement between the BSRN dataset presented in this paper and irradiance spectra from semi-realistic large-eddy simulation up to 0.1 Hz (their Fig. 6). Given the uncertainty regarding how much of the true variability is resolved between 0.1 and 1 Hz, we advise analyses at 1 Hz only in combination with additional constraints, such as knowledge of cloud properties and their velocity, to estimate the fastest possible changes in irradiance. For example, in Mol et al. (2023), knowledge of wind speed, cloud size distributions, and cloud edge transparency were combined to utilize the data down to 1 Hz.

One last implication is that the slow response time slightly reduces the contrast in data. Ehrlich and Wendisch (2015) demonstrated a reconstruction technique of the true 1 Hz signal through deconvolution, essentially a sharpening technique, which is not a trivial exercise. We have not applied this here, as we cannot validate if it works reliably for our dataset, but we mention it as an option to anyone who might want to attempt to apply the method despite its challenges. The only preprocessing that we apply to the data, namely gap filling and quality control, is discussed in Sect. 3.1.

2.2.1 Solar position and direct horizontal irradiance

Information about the Sun's position is important for quality control, data analysis, and interpretation of results. Calculations of the Sun's position (elevation and azimuth angle) are done using the Pysolar (https://github.com/pingswept/pysolar, last access: 16 May 2023) Python package at a 1 min resolution, linearly interpolated to 1 s. For the purposes of this research, the calculations are indistinguishable from highly accurate peer-reviewed code such as the Solar Position Algorithm (SPA; https://midcdmz.nrel.gov/spa/, last access: 16 May 2023). Using the solar elevation angle α (degrees above horizon, or $\mathit{\theta }=\mathrm{90}-\mathit{\alpha }$ as a zenith angle), we calculate the horizontal component of direct irradiance: DHI = DNI ⋅ sin(α). An alternative calculation is DHI = GHI − DIF, which (in the case of a good measurement setup) should be equal to DNI ⋅ sin(α) and is the basis for one of the checks in data quality control (discussed in Sect. 3.1).

2.2.2 Clear-sky irradiance and atmospheric composition

Clear-sky global horizontal irradiance (GHI_cs) is the total downwelling horizontal solar irradiance in the absence of clouds. We use Copernicus Atmosphere Monitoring Service (CAMS) McClear version 3.5 (released in September 2022; Gschwind et al., 2019) as the GHI_cs reference for our dataset. CAMS McClear includes corrections based on atmospheric composition reanalysis, such as aerosols and total column atmospheric water vapor. This allows us to define times when GHI exceeds GHI_cs only through the effect of clouds as best we can, as opposed to simpler techniques that do not correct for aerosols and chemical composition. Atmospheric composition input for McClear is included in their publicly available dataset (https://www.soda-pro.com/web-services/radiation/cams-mcclear, last access: 16 May 2023), which we add to our dataset for context. The only further processing applied to these data is linear interpolation from 1 min to 1 s to match the irradiance observations.

2.3.1 Wind profiles from Cabauw

A 213 m high tower provides wind speed and direction measurements at 2, 10, 20, 40, 80, 140, and 200 m above ground level at a 10 min intervals (Wauben et al., 2010). Figure 1b shows the tower with respect to the BSRN site, which is a few hundred meters to the south. We apply no further processing to these data, apart from creating daily files from the monthly files. The original tower data (including temperature, visibility, and humidity) are shared publicly by the KNMI on their open-data platform: https://dataplatform.knmi.nl/dataset/cesar-tower-meteo-lb1-t10-v1-2 (last access: 16 May 2023​​​​​​​).

2.3.2 NubiScope

To obtain detailed cloud cover observations for analysis, validation of satellite observations (Sect. 2.4), and irradiance-based sky-type classification (Sect. 3.3.2), we use the NubiScope located within a few meters of the BSRN instrumentation (Wauben et al., 2010). In 10 min, this instruments makes a hemispherical scan of the sky using infrared sensors to determine the cloud fraction and various categories of sky type. For validation and analyses purposes, we focused on a 3-year subset (2014–2016) of NubiScope data. If necessary, additional data (May 2008 to April 2017) are publicly available at https://dataplatform.knmi.nl/dataset/cesar-nubiscope-cldcov-la1-t10-v1-0 (last access: 16 May 2023​​​​​​​). Again, we applied no further processing, apart from turning monthly files into daily files.

3.2.1 Satellite processing

We classify cloud types using a simple cloud-top pressure (CTP) and cloud optical thickness (COT) categorization (see the International Satellite Cloud Climatology Project, ISCCP, algorithm description in NOAA, 2022, their Fig. 20). The nine cloud types are cumulus (Cu), stratocumulus (Sc), and stratus (St) for low-level clouds; altocumulus (Ac), altostratus (As), and nimbostratus (Ns) for mid-level clouds; and cirrus (Ci), cirrostratus (Cs), and cumulonimbus (Cb) for high-level clouds. Cu, Ac, and Ci are the optically thinnest clouds for each altitude, and St, Ns, and Cb the thickest; as the only exception in this list, Cb spans from low to high altitude. In this study, we used the abovementioned classification to group cloud conditions of various altitudes and optical thicknesses together in a more intuitive way, although analyses can be done on the input COT and CTP data rather than the derived cloud types. The main limitations are that both the cloud fraction and actual cloud optical thickness contribute to a higher reported optical thickness in a satellite pixel, which is a result of limited spatial resolution, and higher clouds can obscure lower clouds. The spatial satellite product is converted to a time series representative of the BSRN station by determining the most common cloud class within a 5, 10, or 15 km radius of Cabauw, as illustrated by the circles in Fig. 1a. A smaller radius is not possible due to satellite resolution (pixel area of ≈20 km²), and larger radii become unrepresentative for Cabauw. Cloud cover is derived by calculating the fraction of pixels with clouds within a given radius, which is likely an overestimation due to sub-pixel cloud fractions not always being 1. However, overall agreement with the NubiScope is not bad, as illustrated by the similar probability densities in Fig. 4. Correlation coefficients between the two only show marginal improvement between 10 and 15 km. The satellite-derived cloud cover overestimates the extremes at [0.0–0.1) and [0.9–1.0] compared to values bins of [0.1–0.2) and [0.8–0.9), and it does not have the nuances that the NubiScope can resolve. For r=5 km, there are only four pixels, so cloud cover from this is too coarse for most applications, but the dominant cloud type derived from this narrow area around Cabauw is expected to be most representative. This might change for high-altitude clouds at low solar elevation angles, for example, and may perhaps require more sophisticated approaches; therefore, we have included the original spatial satellite fields in our dataset.

Figure 4Comparison of the cloud fraction derived from satellite to the ground-based NubiScope. The analysis is done for radii of 5, 10, and 15 km and shows the probability density for 10 bins with cloud fractions between 0 and 1. Satellite radii range from 5 km (dark, small) to 15 km (light, large), and are shown using the cross markers. Correlation coefficients are shown in the top left for each radius. Data range from January 2014 to December 2016 and are interpolated (using nearest-neighbor interpolation) to a common 5 min time axis.

Download

3.2.2 Validation dataset

Employing the clear-sky and overcast classifications based on a combination of the NubiScope and satellite data, we derive a validation dataset to be used for statistical verification of sky types based on irradiance observations (Sect. 3.3.2). Because both instruments have their limitations, the idea is to only identify a situation as clear sky or overcast if both datasets are in agreement. The rest of the preprocessing involves interpolation to a common 1 min resolution grid, and we also mask out data if either of the two products is missing for given a given time. The validation dataset is provided in Mol et al. (2022) for all three radii, although we mostly use r=10 km, for the years 2014 to 2016. Disagreement between the two observational datasets is common, with the NubiScope being roughly 3 or 1.5 times more conservative with respect to classifying a sky type as clear or overcast, respectively. This is likely not only to do with the difference in the type of observation (ground-based versus remote sensing) but is also due to the fact that the NubiScope is a more sensitive instrument, as we illustrate in Fig. 4 and discuss in Sect. 3.2.1. In most cases, the NubiScope is more conservative; thus, both the clear-sky and overcast conditions are mostly controlled by what the NubiScope sees. Figure 8 illustrates this best, with the satellite and NubiScope differing with respect to the seasonal cycle for clear-sky conditions and with respect to yearly averages for both clear-sky and overcast conditions.

3.3.1 Cloud shadow and enhancement

Broken cloud fields making patterns of cloud shadows and cloud enhancements are one of the most noticeable drivers of intra-day irradiance variability (e.g., Yordanov et al., 2015; Gueymard, 2017; Veerman et al., 2022). Cloud shadows occur when (most) direct irradiance is blocked, and cloud enhancement occurs when light scattered by clouds locally coincides with direct irradiance to increase global horizontal irradiance above clear-sky values. We define a shadow as a region where direct normal irradiance (DNI) is below 120 W m⁻², which is the inverse of what the World Meteorological Organization (WMO) defines as sunshine (DNI ≥ 120 W m⁻²; WMO, 2014), as this is a straightforward implementation. Cloud enhancement requires a more careful approach. We define the occurrence of cloud enhancement, using a single measurement, as when the global horizontal irradiance (GHI) exceeds the reference clear-sky irradiance (GHI_cs). In reality, observed GHI can still fluctuate noticeably under cloud-free conditions, and the GHI_cs reference may not be perfect; therefore, there is some uncertainty in the detection for weak cases of cloud enhancement. To prevent false positives in the detection algorithm, we first apply an activation threshold defined by the GHI exceeding GHI_cs by 1 % and 10 W m⁻². Both a relative and absolute threshold are used, as clear-sky irradiance ranges from 10⁰ to 10³ W m⁻² as function of the solar elevation angle (i.e., time of day). A value of 10 W m⁻² is based on 1 % of the typical order of magnitude for clear-sky irradiance around noon for Cabauw. When the threshold is reached, adjacent measurements are also marked as “cloud enhancement” as long as they exceed GHI_cs by 0.1 %. Edge cases at low solar elevation angles are removed by requiring the DNI to be at least 10 W m⁻². All of these thresholds are chosen to enable us to capture all but the weakest of cloud enhancements, which are arguably not important. The residual third class is called “sunshine” and is based on the WMO definition of sunshine minus cloud enhancement. Detection criteria can be adjusted in the code and recalculated, or another level of filtering can be applied after classification through the derived event statistics (discussed in Sect. 3.4). Examples of the classifications are shown in the top color-coded bar beneath the time series in Figs. 5, 6, and 7.

3.3.2 Overcast, clear-sky, and variable irradiance

The second group of classifications represents the irradiance “weather” type, or sky type, based on irradiance data only. The weather types are clear sky and overcast for smooth and predictable surface irradiance, and a third class of “variable” irradiance represents pronounced and unpredictable 3D radiative effects to contrast the former two. The way that these classifications are derived is partially based on subjective thresholds and assumes good-quality clear-sky data; thus, we validate against satellite and ground-based cloud cover observations.

Table 1Skill scores for irradiance-based sky-type classifications compared to the validation dataset (satellite + NubiScope). Scores are based on a contingency table approach. POD is the probability of detection and FAR is the false alarm ratio.

Download Print Version | Download XLSX

We classify those points in the time series for which GHI stays within 3 % or 5 W m⁻² of GHI_cs in a 15 min centered moving window, with a maximum standard deviation of the ratio $\mathrm{GHI}/{\mathrm{GHI}}_{\mathrm{cs}}=\mathrm{0.01}$ within that window, as clear sky. This irradiance-based algorithm emphasizes smoothness, and thus predictability, more than exactly matching GHI_cs, so as to not rely too much on CAMS McClear being perfectly accurate. Clear-sky conditions are uncommon, occurring between 5 % and 15 % of the time depending on the observational method (Fig. 8). Skill scores indicate that the irradiance-based classification misses almost half of the cases (probability of detection close to 50 %) and is generally too conservative (bias <1) with respect to the validation dataset (see Table 1). The negative bias against the validation dataset, which is not what Fig. 8 shows, is due to the skill scores only being calculated for cases in which there is agreement between the satellite and NubiScope, rather than for all available data in Fig. 8. The order of magnitude of occurrence is similar to what the validation dataset suggests, but the seasonal cycle is not reproduced, although seasonality between the irradiance classification and satellite alone is similar. For 2014 and 2015, the seasonal correlation to the NubiScope is also much better, but there were many cases in 2016 where the NubiScope saw thin cirrus (cloud cover < 5 %, 17, 18, 23, 24, and 25 August) rather than clear sky. In all of these cases, GHI < GHI_cs, which is consistent with thin cirrus attenuating incoming solar radiation slightly; these instances are not part of the validation set due to disagreement between satellite and NubiScope measurements. It appears (from Table 1 and Fig. 10) that there is poor skill in the irradiance-based classification, although manual inspection of quicklooks gives a different impression and most of the bias shown in Fig. 10 appears to stem from cases with thin cirrus. If one wants to filter out the thin-cirrus cases, the classification threshold can be made more strict, thereby limiting cases to more true clear-sky conditions. Figure 6 shows a case with overall agreement between the NubiScope-, satellite-, and irradiance-based clear-sky classifications. We refer the reader to the public dataset (Mol et al., 2022) with time series quicklooks for many more examples.

Figure 8Comparison of sky-type classifications based on satellite time series (r=10 km), NubiScope, and 1 Hz irradiance observations. The validation period is from 2014 to 2016, in units relative to available data during daylight (solar elevation angle α>0^∘).

Download

We define overcast weather as a period of 45 min during which the sum of DNI is below 1 % of GHI_cs and the average is below 10 W m⁻² (in order to catch rare edge cases). This class is indicative of continuous, optically thick, and persistent cloud cover, which is a common occurrence in the Netherlands (20 %–50 % of the time depending on the season), and it does well against the validation dataset. The probability of detection is high (92 %), with 6 % false alarms and only a slight negative bias of −1.7 %. Scores move slightly toward a positive (+2.6 %) or negative (−5.1 %) bias when shortening or lengthening the moving window to 15 and 60 min, respectively. Although a cloud cover of 100 %, as seen by the NubiScope and satellite, is classified as overcast in the validation dataset, cloud cover does not imply that the optical thickness is high enough to block all direct irradiance. This distinction is illustrated in Fig. 6, which shows good agreement for overcast conditions overall, although there is still some irradiance variability with 100 % cloud cover around 10:00 UTC. This example emphasizes our definition of overcast as smooth and predictable diffuse irradiance weather, as opposed to a sky type with 100 % cloud cover. The seasonal cycles between the irradiance-based sky type and NubiScope correlate well (Fig. 10), whereas the satellite yearly cycle is less pronounced.

Finally, the variable weather class, which is built upon the instantaneous classifications defined in Sect. 3.3.1, is defined as any 60 min window in which 10 transitions from a shadow to cloud enhancement or vice versa occur. It is indicative of weather associated with a characteristic bimodal distribution of irradiance, cloud enhancements, and at least a handful of large fluctuations in a short time frame, all of which current numerical weather prediction models cannot reproduce. This classification does well with respect to locating highly variable irradiance conditions, examples of which are shown in Figs. 5 and 7, and can, for example, be used to find case studies.