Irradiance and cloud optical properties from solar photovoltaic systems

Solar photovoltaic power output is modulated by atmospheric aerosols and clouds and thus contains valuable information on the optical properties of the atmosphere. As a ground-based data source with high spatiotemporal resolution it has great potential to complement other ground-based solar irradiance measurements as well as those of weather models and satellites, thus leading to an improved characterisation of global horizontal irradiance. In this work several algorithms are presented that can retrieve global tilted and horizontal irradiance and atmospheric optical properties from solar photovoltaic data and/or pyranometer measurements. The method is tested on data from two measurement campaigns that took place in the Allgäu region in Germany in autumn 2018 and summer 2019, and the results are compared with local pyranometer measurements as well as satellite and weather model data. Using power data measured at 1 Hz and averaged to 1 min resolution along with a non-linear photovoltaic module temperature model, global horizontal irradiance is extracted with a mean bias error compared to concurrent pyranometer measurements of 5.79 W m⁻² (7.35 W m⁻²) under clear (cloudy) skies, averaged over the two campaigns, whereas for the retrieval using coarser 15 min power data with a linear temperature model the mean bias error is 5.88 and 41.87 W m⁻² under clear and cloudy skies, respectively.

During completely overcast periods the cloud optical depth is extracted from photovoltaic power using a lookup table method based on a 1D radiative transfer simulation, and the results are compared to both satellite retrievals and data from the Consortium for Small-scale Modelling (COSMO) weather model. Potential applications of this approach for extracting cloud optical properties are discussed, as well as certain limitations, such as the representation of 3D radiative effects that occur under broken-cloud conditions. In principle this method could provide an unprecedented amount of ground-based data on both irradiance and optical properties of the atmosphere, as long as the required photovoltaic power data are available and properly pre-screened to remove unwanted artefacts in the signal. Possible solutions to this problem are discussed in the context of future work.

1 Introduction

An accurate determination of incoming solar radiation at the Earth's surface is important not only for both climate and weather research, but also for the stable operation of the electricity grid in the future. In Germany alone there are 2.6 million photovoltaic (PV) systems installed, with a nominal power of 71 GWp (Holm, 2023) so that accurate forecasts of solar PV power generation are indeed becoming indispensable for cost-effective grid operation. In this context the proliferation of PV systems provides a unique opportunity to characterise global irradiance with unprecedented spatiotemporal resolution, which would lead to improvements in both weather and climate models. Solar panels can in this way be seen as a dense network of sensors for atmospheric optical properties. This new information could facilitate the development of highly resolved PV power forecasts and can play a role in improving climate models, in particular since the highly variable nature of cloud cover as well as uncertainties in cloud microphysics result in the greatest uncertainty in our understanding of the radiative forcing of the climate.

It has been shown by several authors (see for instance Urraca et al., 2018; Ohmura et al., 1998; Frank et al., 2018; Zubler et al., 2011) that the estimates of global horizontal irradiance (GHI) from both the global ECMWF (ERA5) and the regional (COSMO-REA6) numerical weather prediction (NWP) model reanalyses deviate from ground-based measurements. In Urraca et al. (2018), comparisons are made with pyranometer measurements from the Baseline Surface Radiation Network (BSRN) (Ohmura et al., 1998) and from a dense network of pyranometers operated by European meteorological services. In general the model reanalyses overestimate the irradiance under cloudy skies, which is largely due to an underestimation of cloud optical depth (COD). The mean positive bias of ERA5 daily mean irradiance is +4.05 W m⁻² (3.47 %) over Europe and +4.54 W m⁻² (2.92 %) worldwide. On the other hand, the regional COSMO-REA6 dataset underestimates GHI on clear-sky days, with a mean bias of −5.29 W m⁻² (−3.22 %), which can be attributed to the use of an aerosol climatology with a too large aerosol optical depth (AOD), as discussed in Frank et al. (2018). Although the COSMO-D2 data use a different aerosol scheme, these negative biases in the GHI are still present, especially in summer (Zubler et al., 2011). Satellite datasets perform a lot better, with data from the Surface solAr RAdiation Heliosat (SARAH) showing a mean bias of only +0.86 W m⁻² in the daily mean GHI (compared to +4.22 W m⁻² from ERA5) over Europe (Urraca et al., 2018). Interestingly SARAH overestimates in most cases, with only some stations showing a negative bias related to snow detection. Overall the satellite measurements display a smaller absolute error than reanalysis products. The positive bias of the GHI from satellite retrievals is confirmed by Yang and Bright (2020): their comprehensive global evaluation of hourly satellite irradiance data reveals a mean bias error¹ of 4.67 W m⁻² for hourly SARAH-2 irradiance compared to the nine BSRN stations over Europe (excluding the Austrian station Sonnblick at 3100 m altitude), compared to 7.93 W m⁻² for the Copernicus Atmospheric Monitoring Service (CAMS) radiation data (see Sect. 3.3).

The idea of using PV systems as radiation sensors has been explored by several authors. In Engerer and Mills (2014), Killinger et al. (2016), and Marion and Smith (2017), methods are developed in order to use the output of one PV system to predict that of another, which is in essence done by inferring GHI from PV power measurements. In all three cases empirical models for the decomposition of irradiance into direct and diffuse components are used, and system parameters such as orientation and PV module efficiency are required inputs. Engerer and Mills (2014) achieve a mean bias error of 1.09 % for the PV power output under clear-sky conditions, but the accuracy diminishes for partly cloudy skies, as expected; Killinger et al. (2016) achieve a mean bias error between −3.9 % and −9.8 % for the GHI, depending on the empirical model used for irradiance transposition, and Marion and Smith (2017) achieve a mean bias error for the GHI of within ±1.5 % using south-facing PV modules at 10, 25, and 40^∘ tilt angles. A similar approach is taken in Elsinga et al. (2017), in this case using a single-diode PV model and a different decomposition model. In Halilovic et al. (2019) the authors replaced the original iterative approach used in Killinger et al. (2016) by an analytical method, to minimise computational cost, and achieved a mean bias error of between 0.1 % and 2.1 % for the resulting GHI, using data from silicon reference cell measurements at different tilt and azimuth angles.

In Nespoli and Medici (2017) a different method is introduced, in this case without the need for system-specific information such as orientation or nominal power. A similar approach is taken in Saint-Drenan (2015) and Saint-Drenan et al. (2015), where system parameters are estimated by statistical methods. In addition, Scolari et al. (2018), Laudani et al. (2016), Carrasco et al. (2014), and Abe et al. (2020) have also described the inference of solar irradiance from PV current and voltage measurements using an equivalent-circuit model. In this case greater accuracy is achievable, provided the module temperature is also measured.

This work builds upon the proof-of-concept study presented in Buchmann (2018) (for clear-sky days only); however it is unique in that empirical models for the separation of radiation components are avoided – rather an explicit simulation of the diffuse radiance distribution is performed using libRadtran (Mayer and Kylling, 2005; Emde et al., 2016). Although this is computationally more intensive, it has several advantages over the usual approach (see for instance Perez et al., 1992): by using a state-of-the-art radiative transfer code one can more accurately model the clear-sky irradiance, especially for larger solar zenith angles, and one can explicitly take into account information on aerosol load or ground albedo from freely available datasets. In addition it is possible to include information on the state of the atmosphere from weather models, which is particularly relevant for including the effects of precipitable water on incoming irradiance. The radiative transfer solvers DISORT (DIScrete Ordinate Radiative Transfer) (Stamnes et al., 1988; Buras et al., 2011) and MYSTIC (Monte Carlo code for the physically correct tracing of photons in cloudy atmospheres) (Mayer, 2009) are used for forward model calibration as well as for inferring atmospheric optical properties and GHI from ground-based irradiance measurements and PV power data.

In order for a PV-based determination of solar irradiance to viably complement the global coverage of state-of-the-art satellite measurements, a mean bias error of the order of 5 W m⁻² would be desirable (see the discussion on CAMS and other satellite-based products above). This level of accuracy also corresponds to the target accuracy for global radiation measurements from the BSRN (McArthur, 2005). However, even if this is not achieved, ground-based irradiance measurements and retrievals can be seen as complementary since they have the added advantage of superior spatiotemporal resolution. The first step to achieve this is to accurately model the generated power as a function of system-specific parameters, such as the array's elevation and azimuth angle, conversion efficiency, and temperature dependence, and then extract those parameters from measured power data using a fitting procedure. In order to remove any biases related to atmospheric conditions, it makes sense to first calibrate the systems under clear skies. Once this has been done to sufficient accuracy one can use measured PV power to infer atmospheric optical parameters such as aerosol or cloud optical depth under different sky conditions, enabling the inference of GHI as well as in some cases of direct and diffuse irradiance components.

The more parameters used to model the PV power, the greater the uncertainty in the retrieved irradiance. For this reason it is of course desirable to obtain as much a priori metadata about the PV systems as possible, such as datasheet parameters and array orientation. However, this information is not always readily available, especially when considering a large number of PV systems over a wide area. In that sense, there will always be a trade-off between quantity and quality of the data, which then plays itself out in the accuracy of the retrieved irradiance. The advantage of PV systems or any ground-based devices is that one can achieve a much higher spatiotemporal resolution compared to satellite data or weather models, which thus allows one to study high-frequency fluctuations in global irradiance.

In the European context, irradiance variability is dominated by the optical properties of clouds and less by those of aerosols. Ground-based COD retrievals using broadband measurements from pyranometers have been carried out in several studies (see for example Leontyeva and Stamnes, 1994; Boers, 1997; Boers et al., 1999; Deneke, 2002). Indeed, the transmission of irradiance through a cloud is the most sensitive to its optical depth and less sensitive to droplet radius, single-scattering albedo, or asymmetry factor (Leontyeva and Stamnes, 1994). In most previous studies the clouds are assumed to be horizontally homogeneous in a plane-parallel atmosphere with 1D radiative transfer, which leads to a bias in the extraction of cloud optical properties, in particular under broken-cloud conditions. By neglecting 3D effects, the horizontal transport of photons is not considered, which however plays an important role in real-life situations. These 3D effects can for example lead to an enhancement of solar irradiance (Schade et al., 2007) so that the GHI exceeds the clear-sky irradiance due to reflected light from the edges of clouds. The inherent 4D variability in clouds also complicates the comparison of ground-based and satellite retrievals of cloud properties, since one compares the time average of a point measurement with a spatially averaged quantity.

The goal of this work is to demonstrate that PV systems can indeed be used as ground-based sensors for GHI as well as to infer the optical properties of the atmosphere, in particular the COD. The first results are presented from two measurement campaigns carried out in autumn 2018 and summer 2019 in the Allgäu region in southern Germany, as part of the research project entitled “Development of innovative satellite-based methods for improved PV yield prediction on different time scales for distribution grid level applications” (MetPVNet) (Meilinger et al., 2021b, a). In Sect. 2 the forward model and its calibration are discussed, and the inversion methods are outlined in detail. Section 3 provides a detailed description of the data from the measurement campaigns. The results are presented in Sect. 4, with a focus on both tilted and horizontal irradiance as well as cloud optical depth under overcast skies, and a summary and conclusions are given in Sect. 5. Further details of the PV modelling aspects and radiative transfer simulation can be found in Appendix A.

2 Photovoltaic power model: calibration and inversion

In order to infer local atmospheric optical properties from measured PV data, accurate modelling of both atmospheric radiative transfer and PV power generation is required. In this section both the PV model and the libRadtran radiative transfer model is described, along with the calibration and inversion procedure.

2.1 Forward model: from atmospheric properties to photovoltaic power

The power generated by a solar PV module depends primarily on incoming short-wave solar irradiance and module temperature, both of which depend on atmospheric conditions. Once this dependence is properly described in a model, PV power and/or current measurements can be used to infer the irradiance in the plane of the array, taking into account the geometry of the system, i.e. the elevation and azimuth angle of the solar panels. After extracting this global tilted irradiance (GTI) from PV data, one can go on to infer atmospheric optical properties such as cloud optical depth and global horizontal irradiance by further inverting the radiative transfer model.

The most physically correct method of modelling the power output of a PV plant is with an equivalent-circuit model that captures the properties of semiconductors, such as the two-diode model (see for instance Mertens, 2014). In this way the temperature dependence of the current and voltage is explicitly defined according to the Shockley equation (Shockley, 1949). A drawback of such models is their computational complexity and reliance on parameters found on module datasheets, which are in the most general case not always available. There are however several parameterised models in the literature that attempt to reduce the power generation equation to a simple relation among incoming plane-of-array irradiance, module area, and temperature-dependent efficiency, with the latter described as a function of ambient conditions. Several such models exist (see Skoplaki and Palyvos, 2009, for a review), with some of the most popular being that of the PV Performance Modeling Collaborative (https://pvpmc.sandia.gov/, last access: 17 October 2023) from Sandia National Laboratories (King et al., 2004, 2007) or the Huld model used in the online PVGIS tool (Huld et al., 2011, https://re.jrc.ec.europa.eu/pvgis.html, last access: 17 October 2023). Since the goal here is an inversion, the choice of model depends on the availability of measured data: in this work and in the context of the AC power data from the MetPVNet campaign, a simplified parametric power model is employed. The model is described here briefly, and more details are given in Appendix A.

In order to correctly capture the effects of the variable solar spectrum one also needs to take into account the spectral response of the PV technology in question (Alonso-Abella et al., 2014), which in the case of an equivalent-circuit model can then be included in the calculation of the photocurrent (see for instance the libRadtran-based spectral PV model described in Herman-Czezuch et al., 2022). In the case of parametrised PV power models, this so-called “spectral mismatch”, i.e. the difference between the entire spectrum of incoming radiation and the range utilised by a certain PV module, is usually simply absorbed into the PV model parameters, leading to a site-specific bias that may not take into account variations in the spectrum from local atmospheric conditions. By using libRadtran for calibration and inversion along with information on the state of the atmosphere from weather models, one can take these variations into account in the radiative transfer (RT) simulation and subsequent inversion, as discussed in Sect. 2.2 below. In particular the water vapour column and aerosol optical depth at each site need to be taken into account (see Sect. A3 in Appendix A).

It can be shown using the diode model (see for instance Sauer, 1994; Abe et al., 2020) that the maximum power point (MPP) current generated by a PV module is linearly dependent on the incident irradiance and only very weakly dependent on temperature. However, the dependence of MPP voltage on temperature (which itself is a function of irradiance) is an order of magnitude greater (roughly −0.4 % K⁻¹) so that this simple linear relationship breaks down when considering the PV power. In this work a parameterised power model is used (see Buchmann, 2018; Skoplaki and Palyvos, 2009; Dows and Gough, 1995), with AC PV power described as²

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathit{P}}_{\mathrm{AC},\mathrm{mod}}& \simeq {\mathit{G}}_{\mathrm{tot},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}(b_{\mathrm{1}}+b_{\mathrm{2}}{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}\\ & +b_{\mathrm{3}}{\mathit{T}}_{\mathrm{ambient}}+b_{\mathrm{4}}{\mathit{v}}_{\mathrm{wind}}+b_{\mathrm{5}}{\mathit{T}}_{\mathit{s}\mathit{k}\mathit{y}}),\end{array}\end{array}$$

in the case of the linear temperature model defined in Eq. (A3) (TamizhMani et al., 2003) or as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{P}}_{\mathrm{AC},\mathrm{mod}}& \simeq {\mathit{G}}_{\mathrm{tot},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}(b_{\mathrm{1}}^{\prime }+{\displaystyle \frac{{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}}{b_{\mathrm{2}}^{\prime }+b_{\mathrm{4}}^{\prime }{\mathit{v}}_{\mathrm{wind}}}}\\ & +b_{\mathrm{3}}^{\prime }{\mathit{T}}_{\mathrm{ambient}}+b_{\mathrm{5}}^{\prime }{\mathit{T}}_{\mathrm{sky}}),\end{array}\end{array}$$

in the case of the non-linear temperature model defined in Eq. (A4) (Faiman, 2008; Barry et al., 2020). This means that the modelled AC power P_AC,mod is a non-linear function of plane-of-array irradiance ${\mathit{G}}_{\mathrm{tot},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}$, together with the effects of ambient temperature T_ambient, wind speed v_wind, and sky temperature T_sky that influence module temperature and thus efficiency. Note that the subscript PV for the tilted irradiance ${\mathit{G}}_{\mathrm{tot},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}$ refers to the fact that only the relevant wavelength (in this case 300 to about 1200 nm for silicon PV modules) is considered, and the subscript τ indicates that transmission through the glass surface of the PV panels has been taken into account with an optical model. Further details of the model employed here are given in Sects. A1 and A2 in Appendix A, and the refractive index n of the glass covering is one of the parameters varied in the optimisation procedure. The subscript SW refers to all incoming short-wave photons – the dependence of the spectral mismatch between the PV and SW irradiance bands on atmospheric water vapour and other factors is discussed in Sect. 2.3.

The parameters $b_i\left(b_i^{\prime }\right)$ and (i=1…5) in Eqs. (1) and (2) depend on nominal power, efficiency, the temperature coefficient for power, and the temperature model parameters and are discussed in more detail in Appendix A, which includes a definition of all parameters listed in Table 2. In practice the module temperature can be either measured or modelled, depending on the availability of measurements and/or meteorological data. Within the PV power models described above, the PV module temperature is a static quantity; i.e. the heat capacity (C) of the PV system is not taken into account. However, when dealing with high-frequency measurements of PV power, it is necessary to employ a dynamic temperature model, as discussed in Barry et al. (2020). The characteristic time constant ($=C/J$, with J being the net thermal energy flux of the PV modules) of typically 10 min means that the large fluctuations in irradiance do not translate directly to module temperature variations; i.e. the temperature response is smoothed out. For simplicity the dynamic temperature model is not included in the present study, since most of the systems had power data collected in 15 min resolution.

2.2 Model calibration under clear-sky conditions

In order to infer the irradiance in the plane of array (GTI) from measured PV power or current, the PV model parameters need to be determined, either from datasheets or with a forward model calibration. This is accomplished using data under clear-sky conditions, together with an accurate simulation of the irradiance, followed by a multi-parameter optimisation to find the parameter values. This section describes the technical details of the clear-sky simulation and the relevant atmospheric input parameters used. Figure 1 displays this procedure graphically, and further explanations are given in the following sections.

Figure 1 Schematic diagram showing the different steps of the calibration procedure, with input data sources in green, model steps and algorithms in blue, simulated parameters in orange, and system parameters (see Table 2) in dark red. Note that in this case only clear-sky days or time periods are considered.

2.2.1 Radiative transfer simulation with libRadtran

The clear-sky simulation of tilted irradiance ${\mathit{G}}_{\mathrm{tot},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}$ is performed with the freely available libRadtran software package (Mayer and Kylling, 2005; Emde et al., 2016), with the input parameters shown in Table 1 and the wavelength range from 300 to 1200 nm for silicon PV applications. The corresponding broadband simulation (${\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}$) is also performed, as an input to the temperature model. Spectral integration is carried out using the Kato parameterisation in order to simplify the effects of water vapour absorption by using the so-called correlated-k approximation (Kato et al., 1999). The DISORT solver allows for an explicit calculation of the diffuse radiance distribution on a predefined lattice of elevation and azimuth angles, and the pseudospherical approximation is employed so that only radiative transfer calculations at solar zenith angles (SZAs) of up to 80^∘ can be reliably performed. The Python package PyEphem (Rhodes, 2022) is used to accurately determine the sun position for the corresponding latitude, longitude, and time coordinates. Consortium for Small-scale Modelling (COSMO) model data (see Sect. 3.2) are interpolated by the package cosmomystic (see Barry et al., 2023a) in both time and space in order to create atmosphere profile files suitable as input for libRadtran, in 15 min resolution. In this way variations in water vapour and other atmospheric trace gases are taken into account, and the atmospheric layers are cut off at the appropriate altitude of each site. Concurrent measurements by an AErosol RObotic NETwork (AERONET) sun photometer (Holben et al., 1998) are used to extract the Ångström exponent α and turbidity coefficient τ_a,1 with the aeronetmystic (see software supplement) software package. Other aerosol optical properties such as the single-scattering albedo and asymmetry factor are taken from the Optical Properties of Aerosols and Clouds (OPAC) species library with the option “continental average” (see Table 3 in Hess et al., 1998).

(Rhodes, 2022)(Baldauf et al., 2011)(Holben et al., 1998)(Hess et al., 1998)

Table 1 Model parameters for the libRadtran simulation of clear-sky days, including information on their source.

In order to speed up the simulation the code is parallelised to run on multiple processors: the simulation times are divided into batches, and libRadtran is then called multiple times as a subprocess from Python. In this way the clear-sky simulation takes approximately 1 s per time step (8 s on an 8-core machine), using a diffuse radiance field of 5^∘ resolution in both elevation and azimuth angle, atmosphere files modified from COSMO data, and modified aerosol inputs from AERONET and OPAC.

2.2.2 Non-linear optimisation for PV system parameters

The simplified parametric model described above can be written as

$$\begin{array}{}\text{(3)}& {\mathit{P}}_{\mathrm{AC},\mathrm{mod}}\equiv \mathit{y}=\mathbf{F}(\mathit{x},\mathit{h})\end{array}$$

for the forward model F described by Eq. (A1) and state space defined by (see Table 2)

$$\begin{array}{}\text{(4)}& \mathit{x}\equiv \left(\mathit{\theta },\mathit{\varphi },n,s,\mathit{\zeta },u_i\right)\end{array}$$

so that the calibration procedure is effectively a non-linear, multi-parameter optimisation problem with eight (for the non-linear temperature model)³ or nine (for the linear temperature model) unknowns. As shown in Table 3, the parameter space h in Eq. (3) contains the irradiance proxy from the libRadtran simulation as well as temperature and wind speed data from either the COSMO model or the measurements, which are interpolated to 15 min resolution. In addition the measured sky temperature (see Sect. 3.1) is used. This inversion problem can be solved with the methods detailed in Rodgers (2000). In this case the Levenberg–Marquardt algorithm is used, with the Jacobian matrix calculated explicitly at each iteration.

(Barry et al., 2020)

Table 2 List of PV model parameters in x. In the calibration procedure, the parameters in x known to a certain degree (from datasheets or other sources) of accuracy are fixed, whereas all others are varied.

Table 3 List of additional inputs in h used for calibration on clear-sky days. The subscripts PV and SW refer to the different wavelength bands used for integration; see the discussion in Sect. 2.2.1.

If one varies all parameters in x it quickly becomes apparent that there are not enough degrees of freedom in the signal to uniquely determine a solution with the Bayesian formalism, since several parameters are highly correlated with each other, for instance the inclination angle θ with the scaling factor s or the orientation angle ϕ with the temperature model parameters or the coefficient ζ (see the discussion in Sect. 4.1). It is thus expedient to extract the temperature model parameters separately using the measured module temperature for different PV system configurations (if available) and then fix those parameters in the optimisation procedure. In Barry et al. (2020) a dynamic model was developed by fitting the measured and modelled module temperatures using three different PV systems with different mountings. These results (for the static model case) are used in the overall optimisation, where appropriate.

The calibration algorithm is designed to allow certain parameters to be fixed if they are known, whereas unknown parameters are varied with a given a priori error, which in turn affects the parameter retrieval error and thus propagates into the uncertainty in the inferred irradiance.

2.3 Model inversion under all sky conditions

The calibrated PV systems can now be used as sensors to extract information about the state of the atmosphere. This section describes the different methods used to infer both irradiance and atmospheric optical properties from PV power data, as summarised in the schematic diagram in Fig. 2. The method employed depends on the prevailing weather conditions, specifically on the amount (and type) of cloud cover.

Figure 2 Schematic diagram showing the different steps of the inversion procedure, with input data sources in green, model steps and algorithms in blue, simulated and retrieved parameters in orange, and system parameters in dark red. Note that in this case all atmospheric conditions (all sky) are considered.

In a nutshell, using a 1D DISORT-based method, one can use GTI to extract AOD or COD during clear or completely overcast periods, respectively, which thus allows the determination of the direct and diffuse irradiance components. Under broken-cloud conditions, a 3D MYSTIC-based method allows one to determine the GHI from GTI directly. The DISORT and MYSTIC lookup tables (LUTs) are provided in an open data repository (Barry et al., 2023b).

2.3.1 Global tilted irradiance from PV model inversion

Once the PV system has been calibrated under clear-sky conditions, the system parameters can be fixed, and the measured PV power can be used to infer the global tilted irradiance (GTI; also denoted as ${\mathit{G}}_{\mathrm{tot},\mathrm{SW}}^{\mathrm{\angle }}$) under all sky conditions. The temperature model makes use of broadband irradiance (see Eqs. 1 and 2), whereas the PV power model uses only the relevant spectral range of silicon PV modules (300–1200 nm) so that the spectral mismatch between the light converted to electricity (${\mathit{G}}_{\mathrm{tot},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}$) and the entire short-wave spectrum needs to be taken into account when inverting the model chain.

The ratio of PV-relevant (${\mathit{G}}_{\mathrm{tot},\mathrm{PV}}^{\mathrm{\angle }}$) to broadband (${\mathit{G}}_{\mathrm{tot},\mathrm{SW}}^{\mathrm{\angle }}$) tilted irradiance is a function of the system geometry, time of day, and local atmospheric conditions, with the largest contributing factor being the precipitable water in the atmosphere. In order to take this into account, the libRadtran clear-sky irradiance simulations (see Sect. 2.2.1) are used to characterise the ratio

$$\begin{array}{}\text{(5)}& {\mathit{\xi }}_{\mathrm{spec},\mathrm{GTI}}\equiv {\displaystyle \frac{{\mathit{G}}_{\mathrm{tot},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}}{{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}}}=f\left(\mathbf{\Theta },\left[{\mathit{H}}_{\mathbf{2}}\mathit{O}\right],{\mathit{\tau }}_{\mathrm{a}}\right)\end{array}$$

as a function of incident angle Θ, precipitable water H₂O, and aerosol optical depth τ_a, for each station and measurement campaign. In this way the available information on the water vapour column and aerosol extinction from the COSMO model and AERONET can be taken into account in the PV model inversion. Details are given in Sect. A3 in Appendix A. The fitting function could in principle be extended to include ozone column abundance, which is however not included here, since this information is not available from the COSMO model data. Note that although this method does not take into account the effect of clouds on the spectral mismatch, it is a good first approximation, which will be improved upon in future work (see also Rivera Aguilar and Reise, 2020, for an alternative method). In addition one could modify this algorithm to include operational satellite retrievals of atmospheric parameters such as ozone concentration, if required.

Once the spectral mismatch factor ξ_spec,GTI has been calculated, the next step is to extract the plane-of-array irradiance from the PV power, which in the case of the models given in Eqs. (1) and (2) is simply the solution to the quadratic equations in ${\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}$; i.e.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\mathit{\xi }}_{\mathrm{spec},\mathrm{GTI}}{b}_{\mathrm{2}}{\left({\mathit{G}}_{tot,SW,\mathit{\tau }}^{\mathrm{\angle }}\right)}^{\mathrm{2}}\\ & +\left(b_{\mathrm{1}}+b_{\mathrm{3}}{\mathit{T}}_{\mathrm{ambient}}+b_{\mathrm{4}}{\mathit{v}}_{\mathrm{wind}}+b_{\mathrm{5}}{\mathit{T}}_{\mathrm{sky}}\right)\\ & {\mathit{\xi }}_{\mathrm{spec},\mathrm{GTI}}{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}-{\mathit{P}}_{\mathrm{AC},\mathrm{meas}}=\mathrm{0}\end{array}\end{array}$$

for the linear temperature model, and

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\displaystyle \frac{{\mathit{\xi }}_{\mathrm{spec},\mathrm{GTI}}}{\left(b_{\mathrm{2}}+b_{\mathrm{4}}{\mathit{v}}_{\mathrm{wind}}\right)}}{\left({\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}\right)}^{\mathrm{2}}\\ & +\left(b_{\mathrm{1}}+b_{\mathrm{3}}{\mathit{T}}_{\mathrm{ambient}}+b_{\mathrm{5}}{\mathit{T}}_{\mathrm{sky}}\right)\\ & {\mathit{\xi }}_{\mathrm{spec},\mathrm{GTI}}{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}-{\mathit{P}}_{\mathrm{AC},\mathrm{meas}}=\mathrm{0}\end{array}\end{array}$$

for the non-linear temperature model. These equations can be solved with the quadratic formula, using the calibrated parameters $b_{\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}}$ ($b_{\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}}^{\prime }$) as defined in Eq. (A5) (Eq. A6) for the linear (non-linear) temperature model; the available data for T_ambient, v_wind, and T_sky; and the spectral mismatch factor for GTI defined in Eq. (5). Note that the inverted ${\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}$ is the irradiance impinging upon the PV module under the glass covering so that the optical model has not been inverted yet. In order to compare this quantity with pyranometers, the transmission of light through the glass τ_PV,rel(Θ) (also known as the “incidence angle modifier”; see Eq. A9) must be taken into account so that the final GTI is given by (see also Eq. A14)

$$\begin{array}{}\text{(8)}& {\mathit{G}}_{\mathrm{tot},\mathrm{SW}}^{\mathrm{\angle }}={\displaystyle \frac{{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}}{{\mathit{\tau }}_{\mathrm{PV},\mathrm{rel},\mathrm{eff}}}}.\end{array}$$

For the direct extraction of GTI an empirical formulation is used to find the effective incidence angle for the diffuse component, whereas for the inversion onto optical properties the refractive index is explicitly taken into account within the radiative transfer simulation. More details are given in Sect. A2 in Appendix A.

2.3.2 Clearness index and irradiance variability classification

Using the global tilted irradiance extracted from the measured PV power data, different methods are used in order to extract atmospheric optical properties and global horizontal irradiance, depending on the prevailing weather conditions. By combining the inverted tilted irradiance with the corresponding clear-sky curve one can calculate a clearness index

$$\begin{array}{}\text{(9)}& {\mathit{k}}_{\mathit{i}}={\displaystyle \frac{{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau },\mathrm{inv}}^{\mathrm{\angle }}}{{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau },\mathrm{clear}}^{\mathrm{\angle }}}}\end{array}$$

for each time step, allowing the data to be separated into clear, overcast, and broken-cloud time periods. The clearness index is then used to estimate the cloud fraction, which is discussed in more detail in Sect. 2.3.3.

On clear days (or during clear time periods) the aerosol optical depth (AOD) can be inferred, whereas under cloudy conditions the cloud optical depth (COD) can be found, depending on the degree of cloud cover. In this work the extraction of COD using a DISORT-based LUT under completely overcast skies is examined in more detail in Sect. 2.3.3. An in-depth analysis of aerosol optical properties will be carried out in future work. For broken-cloud conditions, a MYSTIC-based LUT is used to infer the global horizontal irradiance from tilted irradiance measurements, as discussed in Sect. 2.3.5 (see Chap. 9 of Meilinger et al., 2021b).

2.3.3 Cloud optical depth with DISORT lookup table

Cloud optical properties are functions of microphysical properties such as the droplet size distribution, droplet number concentration, and thermodynamic phase. For water clouds, the absorption and scattering of solar irradiance can be efficiently characterised (Hu and Stamnes, 1993) by the effective radius r_eff and cloud liquid water content (LWC), which can be related to cloud optical depth (COD, τ_c) through

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{\mathrm{c}}={\displaystyle \frac{\mathrm{3}\mathrm{LWP}}{\mathrm{2}{r}_{\mathrm{eff}}{\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}},\end{array}$$

where the liquid water path (LWP) is the integral of the LWC across the height of the cloud. The derivation of this equation (see for instance Petty, 2006) assumes large Mie extinction, which is justified since clouds appear to be (mostly) white in the solar spectrum. Although both τ_c and r_eff can be accurately retrieved from spectral measurements of reflected radiation, the transmission of light through clouds is mostly sensitive to the cloud optical depth. This is due to the fact that changes in transmission due to variations in single-scattering albedo and the asymmetry factor (which depend on r_eff) are small compared to those due to changes in optical depth (Leontyeva and Stamnes, 1994). For illustration, in the two-stream approximation for conservative scattering (no absorption), the transmittance T can be shown to be (see for instance Petty, 2006)

$$\begin{array}{}\text{(11)}& T={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+(\mathrm{1}-g){\mathit{\tau }}_{\mathrm{c}}}},\end{array}$$

where g is the asymmetry factor. For liquid water clouds, scattering is mostly in the forward direction, with g≃0.85, whereas for ice clouds g≃0.7. In both cases the variations in g are small so that τ_c is the primary factor influencing T. The hyperbolic dependency of T on τ_c means that the transmission curve is rather steep at small optical depths but flattens out for COD ≳15. This has implications for the accuracy of ground-based retrievals, as discussed in detail in Sect. 4.4. It must also be noted that in the algorithm described in Sect. 2.3.1, spectral variations in cloud optical properties are not taken into account. In practice this means that variations in the single-scattering albedo at higher wavelengths around 1 µm (silicon PV modules are still sensitive to wavelengths up to 1.2 µm) may be unaccounted for.

In this work a lookup table for the optical depth of a typical stratus cloud is constructed using DISORT in 15 min time intervals, under the assumption of a pseudospherical or plane-parallel atmosphere with horizontally homogeneous liquid water, clouds, and a completely cloudy sky. This means that 3D effects are not taken into account, and the results need to be interpreted with care, especially in situations with broken clouds. In addition, different cloud types such as thicker cirrus clouds, mixed-phase clouds, or multi-layer clouds are not properly represented by the LUT. Due to the pseudospherical approximation only SZAs up to 75^∘ are considered (for SZAs above 75^∘ with cloud cover the pseudospherical DISORT solver is unstable). The cloud parameters in Table 4 are input into libRadtran, and the COD at 550 nm (τ_c,550 nm) is varied on a 16-step logarithmic scale between COD =0.5 and COD =150, using the default Hu parameterisation (Hu and Stamnes, 1993) and 16 streams. Note that the COD LUT also implicitly contains aerosol information as an input, since here the aerosol properties are fixed using the OPAC database (Hess et al., 1998), and the Ångström parameters from AERONET are used.

Table 4 Cloud parameters for the DISORT simulation of a continental stratus cloud (Hess et al., 1998).

As described for the clear-sky simulation in Sect. 2.2, the direct irradiance and diffuse radiance field are calculated with libRadtran, the latter in this case with a coarser resolution of 10^∘ in both azimuth and elevation angles.⁴ The LUT is then used to find the COD by first calculating the plane-of-array irradiance for the corresponding PV system or pyranometer orientation (see Sect. 2.3.1) and then interpolating the COD in time to match the resolution of the measured data. For this purpose the original 1 Hz pyranometer and PV data are smoothed to 1 min resolution, whereas the low-frequency PV data are kept at 15 min resolution (see Sect. 3.1).

In order to determine the exact time points at which a cloud is above the sensor, the cloud fraction is determined by creating a mask based on the clearness index in Eq. (9) using a threshold of 0.8 and overshoot limit of 1.1; i.e.,

$$\begin{array}{}\text{(12)}& \mathit{c}\mathit{f}=\left\{\begin{array}{lcl}\mathrm{1}& \text{if}& {\mathit{k}}_{\mathit{i}}\le \mathrm{0.8}\\ \mathrm{0}& \text{if}& \mathrm{0.8}<{\mathit{k}}_{\mathit{i}}\le \mathrm{1.1}\\ \text{nan}& \text{if}& {\mathit{k}}_{\mathit{i}}>\mathrm{1.1}.\end{array}\right.\end{array}$$

This binary cloud mask (clear =0, cloudy =1) is then also smoothed with a moving average function over 60 min in order to create an estimate of the cloud fraction (〈cf〉₆₀). Varying the thresholds in Eq. (12) shows that the cloud fraction computed in this way depends less on the exact threshold used but more on the window size chosen for the moving average. Indeed, comparison with concurrent cloud camera retrievals shows that 60 min is a reasonable averaging time to use, when averaging a cloud mask created with data at 1 min resolution. However, the algorithm is limited by the viewing angle of the respective PV system or pyranometer, so it can be inaccurate when there are many clouds on the horizon, for instance.

The COD is then only extracted for data points for which cf=1, i.e. by finding the values of τ_c,550 nm for which

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathrm{meas}/\mathrm{inv}}^{\mathrm{\angle }}& ={\mathit{G}}_{\mathrm{dir},\mathrm{SW},\mathrm{cloudy}}^{\mathrm{\angle }}\left({\mathit{\tau }}_{\mathrm{c},\mathrm{550}\mathrm{nm}}\right)\\ & +{\mathit{G}}_{\mathrm{diff},\mathrm{SW},\mathrm{cloudy}}^{\mathrm{\angle }}\left({\mathit{\tau }}_{\mathrm{c},\mathrm{550}\mathrm{nm}}\right)\end{array}\end{array}$$

for all points under a cloud, where meas or inv refers to measured or inverted GTI from pyranometers or PV systems, respectively. The corresponding direct and diffuse components can then also be extracted from the LUT, although in this case the direct irradiance is basically zero (beneath a cloud).

As mentioned above, this approach is limited by the fact that a 1D radiative transfer solver such as DISORT cannot take into account horizontal transport of photons so that 3D effects such as radiative enhancement under broken-cloud conditions (see for instance Schade et al., 2007) are not taken into account. For this reason only situations with overcast conditions will be considered when applying this method. In situations with low overall cloud cover, the COD is not the main determinant of the total irradiance received by the sensor or PV system but rather the cloud fraction and the AOD. To this end a complementary approach using a MYSTIC-based LUT (see Sect. 2.3.5) is used in order to translate measured tilted irradiance to horizontal irradiance under broken-cloud conditions.

2.3.4 Aerosol optical depth with DISORT lookup table

As mentioned above, in this work the extraction of the AOD is not discussed in detail. However the procedure will briefly be described here, since this is used as an alternative method for determining the GHI from tilted irradiance measurements. An AOD-GTI lookup table can be created in a similar way to the COD LUT described in Sect. 2.3.3, where in this case the AOD at 500 nm is varied on a logarithmic scale in 16 steps between AOD =0.01 and AOD =1. In addition, the aerosol properties are fixed to the so-called continental average scheme from the OPAC database (Hess et al., 1998), and for the inversion procedure the AERONET-based Ångström parameters are not used as input. In this context it must be noted that the typical dust event reaching Europe does not have such a high AOD but is rather characterised by small values of the Ångström exponent of less than 1, indicating the presence of coarser dust particles. For example, one study of the climatology of dust events found a mean AOD of 0.155, 0.32, and 0.122 for dust plumes in southern, central, and northern Europe, respectively (Mandija et al., 2018).

Using the AOD-GTI lookup table, the AOD can be extracted on clear-sky days, and from this also the direct and diffuse irradiance as well as the global horizontal irradiance. In this way the AOD is used as an intermediate step for the reverse transposition of GTI to GHI. In Germany and especially in the Allgäu region the AOD is usually very small (during the measurement campaigns it did not exceed 0.5 at 550 nm) so that any errors in the calibration procedure lead to large relative biases in the AOD. This then leads to biases in the direct and diffuse components, but since there are very few absorbing aerosols, these errors have opposite signs and largely cancel out in the determination of the GHI. In Sect. 4, the GHI extracted via both COD and AOD under different conditions is compared to that measured by pyranometers and satellites, as well as the GHI predicted by the COSMO weather model. However the inferred AOD itself is not examined in detail.

2.3.5 From tilted to horizontal irradiance with MYSTIC lookup table

In order to extract the global horizontal irradiance from the tilted irradiance (from pyranometers or PV systems), a MYSTIC-based LUT for the GHI was developed using large-eddy simulation (LES) cloud fields (Črnivec and Mayer, 2019), taking into account various factors such as albedo, water vapour, sensor geometry, and cloud fraction. Detailed 3D radiative transfer simulations were carried out, and the most important factors simply turned out to be the sensor and sun geometry as well as the cloud fraction. A detailed description of the MYSTIC LUT is given in Chap. 9, Sect. 9.1.5, of Meilinger et al. (2021b).

Table 5 shows the limits of applicability of the MYSTIC LUT, for which there are three major reasons. Firstly, despite several optimisations like the reduction in the number of photons used for the Monte Carlo simulations, the computational demand for calculating the LUT is high. For this reason, 20^∘ is chosen as the lower limit for the SZA, since in the latitudes under investigation no smaller values occur. Similar considerations apply to the tilt angle of PV panels – in Allgäu, Germany, one rarely encounters title angles larger than 50^∘. The second limiting factor relates to the derivation of a cloud mask and cloud fraction from the radiation measurements (see Eq. 12). Firstly, this is only possible when there is a direct line of sight between the sun and the module or sensor, which limits the relative azimuth angle between the sun and the PV panel. Secondly, the derivation of cloud fraction from temporally resolved radiation measurements becomes imprecise at large SZAs for geometrical reasons so that the upper limit of the SZA is set to 60^∘. Finally, the cloud fraction limits are determined by two factors: firstly, the LUT model was developed and tested for partly cloudy situations. The special cases of 0 % and 100 % cloud fraction are considered separately with the DISORT-based LUTs, as other parameters like AOD and COD become relevant here. Secondly, the exact cloud fraction limits (0.13 and 0.82) given in Table 5 are constrained by the available cloud scenes from LES simulations.

Table 5 Limits on the input parameters for the MYSTIC LUT.

The measured or inverted tilted irradiance, together with the average cloud fraction over the last hour (as described above; see Eq. 12), is fed into the MYSTIC LUT, along with the sensor and sun geometry. In this way the GHI can be extracted from the GTI under broken-cloud conditions. This method can however not be used to determine the optical depth, nor can the direct and diffuse irradiance components be separated from each other, since the fit was created using the GTI and GHI.

3 Measurement and validation data

3.1 Ground-based measurements

Model calibration and inversion are performed with PV power data recorded over two measurement campaigns in autumn 2018 and summer 2019 as part of the MetPVNet measurement campaign (see Chap. 3 of Meilinger et al., 2021b). There were a total of 24 stations spread out in the region around Kempten (47.715924^∘ N, 10.314006^∘ E), as shown in Fig. 3, with 22 of them equipped with silicon-based pyranometers measuring both GHI and GTI in the plane-of-array of the PV system. Two master stations (MS01 and MS02) were also equipped with secondary standard pyranometers and pyrheliometers in order to measure both components of the incoming short-wave radiation, cloud cameras, and spectrometers to record spectral information. The MS01 station also contained a sun photometer to determine aerosol properties, as part of AERONET, and a pyrgeometer to measure long-wave downwelling irradiance.

Figure 3 Map showing the locations of PV systems used in the MetPVNet measurement campaign (taken from © Google Earth). The top-left corner is at 47.85^∘ N, 10.09^∘ E and the bottom-right corner at 47.38^∘ N, 10.52^∘ E. The yellow line marks the border between Germany and Austria; the grey line is the border between the states of Bavaria and Baden-Württemberg.

The PV power data were for the most part provided by the local distribution network operator Allgäuer Überlandwerk GmbH (AÜW), recorded in 15 min intervals. These data represent the amount of energy generated in the last 15 min so that care needs to be taken to translate them into a measured power corresponding to a specific time stamp. For this purpose the data are simply shifted by half a period and resampled, since by integration of power over 15 min one effectively smooths the power curve. In addition there were five stations equipped by egrid GmbH with high-frequency power measurement devices: for these stations the power was recorded in 1 Hz resolution.

Analysis of the measured data revealed a total of 12 clear-sky days that occurred between 12 September and 14 October 2018, as well as 9 clear-sky days between 25 June and 13 August 2019, as shown in Table 6. COSMO model data for the corresponding days were procured from Germany's national meteorological service, the Deutscher Wetterdienst (DWD), in order to accurately recreate atmospheric conditions using cosmomystic. These days are used for calibration of each PV system.

Table 6 Dates of clear-sky days during the measurement campaigns in autumn 2018 and summer 2019.

The network of PV systems was equipped with low-cost silicon-based pyranometers, with two devices per station: one mounted in the plane of the PV array and one horizontal, with 1 Hz resolution and an overall accuracy of 5 %. An absolute calibration of these sensors had been carried out at the Leibniz Institute for Tropospheric Research (TROPOS) prior to the campaign by comparing their output to that of a secondary standard pyranometer (2 % accuracy). In order to compensate for errors in mounting the plane-of-array pyranometers, the calibration algorithm described in Sect. 2.2 is also applied to the pyranometer data, in this case without an optical model and only using data up to an SZA of 60^∘. Due to the substantial cosine bias, a correction factor is empirically determined:

$$\begin{array}{}\text{(14)}& \mathit{C}\left(\mathit{\mu }\right)=-\mathrm{3.01}{\mathit{\mu }}^{\mathrm{3}}+\mathrm{5.59}{\mathit{\mu }}^{\mathrm{2}}-\mathrm{3.34}\mathit{\mu }+\mathrm{1.45},\end{array}$$

where μ=cos θ₀ for horizontal sensors and μ=cos Θ for tilted sensors (θ₀ is the SZA and Θ the angle of incidence; see Eq. A8). The pyranometer data are used for comparison with the inverted irradiance (both tilted and horizontal) and for finding atmospheric optical properties using the lookup table method.

In order to validate the PV- and pyranometer-based COD retrievals, it would be appropriate to use another ground-based source of cloud optical properties; however unfortunately there are no appropriate meteorological stations in the immediate area that could have been used for this purpose. Although there are several DWD stations in the Allgäu region (in Kempten, Oberstdorf, and Hohenpeissenberg), these provide information on irradiance (direct and diffuse) but not on cloud optical properties (see Becker and Behrens, 2012). Thus, a true validation would have to be done for another dataset with PV systems closer to a measurement station that has ground-based retrievals of COD. For this reason, the COD retrievals are simply compared to the corresponding cloud properties from weather models and satellite data.

3.2 Weather model data

The Consortium for Small-scale Modelling (COSMO) numerical weather model is a nonhydrostatic regional model developed by the DWD (Baldauf et al., 2011). Note that this model was recently replaced by the so-called ICOsahedral Nonhydrostatic (ICON) model, which has been fully operational since January 2021. Since the measurement campaigns took place in 2018 and 2019, in this work the COSMO-EU model with a spatial resolution of 2.2 km is used, as input to the clear-sky irradiance simulation (see Sect. 2.2.1), for PV model calibration (see Sect. 2.2.2), and for validation and comparison of the inverted irradiance with weather model predictions. For the clear-sky simulation, temperature and pressure profiles as well as the water vapour column are extracted from COSMO data, whereas for both the calibration and the inversion procedure the surface temperature and wind speed are used. For comparison and validation both direct and diffuse downward irradiance data are used.

In order to compare COSMO COD data with the cloud optical depths extracted from the PV systems, a 2D COD field must be computed from the 3D cloud variables generated by the COSMO model. For each grid cell, a cloud fraction variable in COSMO indicates which fraction of the cell is covered by clouds. To derive a vertically integrated COD, an assumption needs to be made as to how these clouds overlap in a model column. Following Scheck et al. (2018), the commonly used random-maximum cloud overlap assumption is adopted, along with the method of Matricardi (2005), in order to compute the vertically integrated COD for a number of subcolumns within each model column. From these subcolumn values a mean COD for the cloudy part of the column is derived. A total COD is then computed as the average of the column mean COD over 5×5 columns centred around the column containing the relevant ground station.

3.3 Satellite data

The Copernicus Atmospheric Monitoring Service (CAMS) radiation service (Qu et al., 2017; Schroedter-Homscheidt et al., 2022) is an online satellite-based and numerical-model-based service with radiation as well as cloud and aerosol data available to download for free, covering the period from February 2004 to the present. The spatial coverage is Europe, Africa, the Middle East, the eastern part of South America, and the Atlantic Ocean, interpolated to the point of interest, and with a time resolution of up to 1 min. In this work the global, direct, and diffuse components of irradiance are imported from CAMS (version 4.0), for each station and for all days in the two measurement campaigns. These data are used as a comparison for the irradiance inverted from PV systems.

In addition to irradiance, CAMS provides data on cloud and aerosol properties. In this work, the cloud parameters from the Advanced Very High Resolution Radiometer (AVHRR) Processing scheme Over cLoud, Land, and Ocean Next Generation (APOLLO_NG) analysis (Kriebel et al., 2003; Klüser et al., 2015) are used, using data from the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) instrument on board the Meteosat Second Generation (MSG) satellite. In this case the so-called Stephens method (Stephens et al., 1984) is used to determine the COD, using a two-stream solution of the radiative transfer equation, along with an updated algorithm using a probabilistic approach for cloud detection (Klüser et al., 2015). For comparison with the COD inferred from the PV systems, APOLLO_NG data are extracted for the closest pixel to each station.

4 Results

The calibration and inversion procedure described in Sect. 2 is applied to the data from the measurement campaign described in Sect. 3 in order to extract irradiance and optical properties from the PV systems in the Allgäu region. After a brief summary of the calibration results in Sect. 4.1, the retrievals of tilted and horizontal irradiance are presented in Sect. 4.2 and 4.3, and the inferred COD results are shown in Sect. 4.4.

4.1 Model calibration and uncertainty

The PV models in Eqs. (1) and (2) are used together with the clear-sky simulation described in Sect. 2.2.1 and the clear-sky days (in Table 6) in order to calibrate each system. In each individual case the days that turned out not to be clear are discarded from the calibration dataset, and the data are restricted to the time periods in which the required inputs such as ambient temperature, atmospheric long-wave irradiance, and wind speed are available. In order to validate the calibration results, the retrieved elevation and azimuth angles are compared to ground truth data from the Bavarian Agency for Digitisation, High-Speed Internet and Surveying (LDBV). The so-called Level of Detail 2 (LoD2) database (https://ldbv.bayern.de/produkte/3dprodukte/3d.html, last access: 24 October 2023) contains a 3D building model constructed using airborne laser scanning so that the roof pitch of individual buildings can be extracted. Figure 4 shows a comparison of the retrieved orientation angles with the ground truth values for each system and using the linear temperature model.

Figure 4 Comparison of the retrieved elevation angles (θ_opt, left plot) and azimuth angles (ϕ_opt, right plot) with the corresponding ground truth angles (θ_actual and ϕ_actual) from airborne laser scanning data, for PV systems with 15 min AÜW (red points) and 1 Hz egrid (blue points) data (cf. Table 7 and the description in Sect. 3.1).

In most cases the algorithm finds reasonable values for the angles: the larger deviations can usually be explained for individual cases; for instance for PV04 the inverter MPP tracking algorithm distorts the clear-sky days, whereas for the systems at PV11 the different PV arrays at the site are not well characterised (see the point labelled 11,3 in Fig. 4). In other cases shading effects played a role: in most cases the calibration performed better when using both summer and autumn data, since in summer the sun is much higher, and shading effects play a smaller role. In general the model calibration works best when using as much data as possible, since there is for instance more variation in temperature in order to find more reliable temperature model parameters.

As discussed in Sect. 2.2.2, several parameters are correlated with each other: the size of the PV system (captured by the factor s) correlates with the tilt angle θ, whereas the azimuth angle ϕ shows a large correlation with the parameters of the temperature model, since the warming up and cooling down of the PV system are delayed with respect to the diurnal variation in solar irradiance. In general the use of measured module temperature leads to better calibration results. It turns out that the calibration algorithm presented here does not perform well when using the non-linear Faiman temperature model and 15 min power data, even though this model couples irradiance and wind speed in a more physically correct way (see for instance Faiman, 2008; Barry et al., 2020). The benefit of this model is lost for coarsely resolved 15 min data so that in the end the algorithm proposed here does not always find an optimal solution, specifically if the temperature model parameters are unknown. The bias that then occurs in the final tilted irradiance inversion results can be seen in the plots in Sect. 4.2 as well as in the results in Tables 8 and 9. However, this bias in the tilted irradiance does not always translate into a bias in GHI, as is seen below. Table 7 lists the PV systems used in this work, along with the corresponding time resolution of their data.

Table 7 List of the PV systems used for this work (see also the map in Fig. 3). The data resolution column indicates whether a particular system is used for this analysis or not, with an explanation given in the cases where the system is omitted. Note that stations MS01 and PV02 had no PV systems, only pyranometers and other measurement equipment. PV11 has four separate PV systems.

Table 8 Mean bias error (in W m⁻²) and relative mean bias error (in brackets in %) of GTI retrievals compared to tilted pyranometers. The values marked with ^* are high due to calibration errors using the non-linear temperature model with 15 min data.

Table 9 Root mean squared error (in W m⁻²) and relative RMSE (in brackets in %) of GTI retrievals compared to tilted pyranometers. The values marked with ^* are high due to calibration errors using the non-linear temperature model with 15 min data.

4.2 Global tilted irradiance from PV power data

In this section the plane-of-array irradiance from PV power retrievals is compared to the tilted pyranometer measurements at selected stations during the two measurement campaigns. The results are obtained using two different approaches for module temperature: (i) the linear temperature model (Eq. A3) and (ii) the non-linear Faiman temperature model (Eq. A4). The results are compared in Tables 8 and 9, and scatterplots for both models are shown. Both stations with 15 min power data as well as those with high-frequency data (1 Hz data, smoothed to 1 min resolution) are included in the analysis, and in each case a comparison is made with the measurements from TROPOS silicon-based pyranometers, except for the MS02 master station, where the tilted Kipp & Zonen CMP11 pyranometer is used for validation. In all cases a limit of 80^∘ is imposed on both the solar zenith angle and the incident angle in order to avoid possible errors from both the radiative transfer simulation and the optical model.

Figures 5 and 6 show a comparison between the retrieved GTI and that measured by pyranometers for the 1 and 15 min data, respectively, for each measurement campaign and using both the linear and the non-linear temperature models. The corresponding statistical measures of the mean bias error (MBE), defined by

$$\begin{array}{}\text{(15)}& \text{MBE}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left(X_{\mathrm{inv},i}-X_{\mathrm{ref},i}\right),\end{array}$$

and the root mean squared error (RMSE), defined by

$$\begin{array}{}\text{(16)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\left(X_{\mathrm{inv},i}-X_{\mathrm{ref},i}\right)}^{\mathrm{2}}},\end{array}$$

for the inverted quantities X_inv,i and the reference quantities X_ref,i are shown in Tables 8 and 9, along with the relative error metrics rMBE and rRMSE calculated by normalising the MBE and RMSE with the mean of the reference quantity, 〈X_ref〉. The scatterplots throughout this work are coloured according to a probability density function calculated using the multi-variate Gaussian kernel density function gaussian_kde in the Python toolbox SciPy, with yellow (light grey) for high- and blue (dark grey) for low-frequency points in the colour (black and white) version. In Figs. 5 and 6 one can see that most points lie close to the 1:1 line, for both campaigns, albeit with a positive bias in all cases. The 1 min data show a slightly larger spread of points than the 15 min data, since in the former case there are more outliers caused by (i) temperature effects, (ii) 3D radiative transfer effects, and (iii) spatial effects due to the differences in cloud cover and the sensor position between the PV and pyranometer. In addition the slightly different geometry of flat PV arrays compared to glass-dome-shaped pyranometers could play a role, especially when it comes to their sensitivity to different viewing angles. Another possible reason for the positive bias could be a systematic bias in the tilted pyranometer measurements, even after the bias correction described in Sect. 3.1.

Figure 5 Combined comparison of GTI retrieved from PV power measurements with that measured by tilted pyranometers, using data in 1 min resolution from MS02, PV11, PV12, PV15, and PV19 (see Table 7), for 2018 (a, c) and 2019 (b, d), together with the linear (a, b) and non-linear (c, d) temperature models. The relative mean bias error (rMBE) and relative root mean squared error (rRMSE) are shown in the inset, along with the mean of the reference GTI and the number of data points used, denoted as 〈G_ref〉 and n, respectively.

Figure 6 The same as Fig. 5 using data in 15 min resolution from MS02, PV01, PV06, PV08, PV10, PV11, PV14, PV16, PV17, PV18, and PV21 (see Table 7) for 2018 (a, c) and 2019 (b, d). Note that the large errors in the non-linear model in the bottom row result from errors in the calibration procedure (see the values marked with ^* in Tables 8 and 9).

The two different temperature models achieve similar results for the 1 min data, with the non-linear model showing an MBE of 5.29 W m⁻² (12.25 W m⁻²) in autumn 2018 (summer 2019) compared to an rMBE of 16.23 W m⁻² (21.12 W m⁻²) for the linear model. In general the algorithm performs worse with 15 min data, which has to do with errors from the calibration procedure and uncertainties in the PV power measurements – the systems with high-frequency measurements are thus in general better characterised and deliver more accurate irradiance retrievals, as shown in Sect. 4.3 below. This effect is quite extreme for the non-linear Faiman temperature model (see the values marked with ^* in Tables 8 and 9 as well as the plots in the lower panels of Fig. 6), since in some cases the calibration algorithm cannot find an optimal solution, and the a priori values have to be relied upon, leading to an average MBE of 91.98 W m⁻².

4.3 Global horizontal irradiance from PV power data

In the following the global tilted irradiance (GTI) retrievals are converted to global horizontal irradiance (GHI) and compared to the measurements from pyranometers as well as to the satellite and weather model data. This conversion is performed in three different ways, depending on the prevailing weather conditions, as described in Sect. 2.3. In the case of clear skies, the DISORT-based LUT is used to find the value of aerosol optical depth, and the inferred AOD is then used to calculate the direct and diffuse horizontal irradiance components, which result in the GHI. The same method is used under cloudy skies using the DISORT-based COD LUT but only for times at which the sensor is under a cloud, using the inferred cloud fraction method as described in Sect. 2.3.3. For this reason, this method tends to underestimate the GHI under broken-cloud conditions, which is discussed in detail below.

The third approach for finding the GHI is using the MYSTIC-based LUT described in Sect. 2.3.5, where the LUT input parameters are simply the array geometry, sun position, and cloud fraction. In this case there are certain restrictions on these parameters, as shown in detail in Table 5. This means that not all of the retrieved GTI data points can be transformed into GHI using this method – in particular the SZA is limited to between 20 and 60^∘ and the cloud fraction to between 0.13 and 0.82 so that neither completely overcast nor clear-sky conditions are taken into account. This method thus deals with the case of mixed-/broken-cloud conditions, in which it is more likely that there will be errors due to 3D effects and the sensor position.

4.3.1 GHI retrieval validation with pyranometer measurements

As discussed in Sect. 4.2 above, the PV systems with 1 min data show the best calibration results and the most accurate tilted irradiance retrievals. The scatterplots in Fig. 7 compare the GHI retrieved from these systems to that measured by horizontal pyranometers, using all three methods and for both temperature models. The statistical measures of the different retrievals are shown in Table 10, where it can be seen that the mean bias error reaches the goal of 5 W m⁻² only in certain cases.

Figure 7 Combined comparison of GHI retrieved from PV power measurements with that measured by horizontal pyranometers, using data in 1 min resolution from MS02, PV11, PV12, PV15, and PV19, for 2018 (left two columns) and 2019 (right two columns), and both the linear and the non-linear temperature models. The top row shows GHI retrieved via AOD under clear skies using the DISORT AOD LUT, the middle row is for cloudy periods via the COD using the DISORT COD LUT, and the bottom row is for broken-cloud periods using the MYSTIC 3D LUT.

Table 10 The mean bias error and root mean squared error (in W m⁻²) along with the rMBE and rRMSE (in brackets in %) of GHI retrievals using the 1D DISORT (AOD and COD) and the 3D MYSTIC LUTs compared to horizontal pyranometers, for 1 min and 15 min data.

Under clear conditions (top row of Fig. 7), the linear model applied to 1 min PV data achieves an rMBE for the GHI of 18.15 W m⁻² (3.9 %) in autumn 2018 and 9.44 W m⁻² (1.4 %) in summer 2019. Note that the mean irradiance is higher in summer, but there are less points that can be classified as clear (n≃9400 compared to n≃13 400 in autumn). The non-linear model performs significantly better in autumn, with an rMBE for the GHI of 1.83 W m⁻² (0.4 %), but in summer the bias is similar to the linear model [9.75 W m⁻² (1.5 %)]. Interestingly the linear temperature model performs better in summer than in autumn, whereas the non-linear model performs better in autumn. This could be attributed to differences or uncertainties in the calibration of the temperature models. In general these results show that the DISORT LUT method performs comparably well for extracting GHI from PV power measurements under clear conditions. In all cases the rRMSE is of the order of 5 % to 7 %, on average 33.19 W m⁻².

It is evident from the middle row of Fig. 7 that the bias is greater under cloudy skies, with an average over both campaigns of 11.29 W m⁻² for the non-linear model and 17.50 W m⁻² for the linear model. In autumn the lower average irradiance of 225 W m⁻² leads to a higher relative MBE than in summer, where the average irradiance is 298 W m⁻². At this point it is worth mentioning that the algorithm only finds the COD and thus the GHI when the sensor is under a cloud, hence the lower average irradiance in comparison with that under clear skies, and the RMSE is higher (66.23 W m⁻² on average) than in the case of clear skies, where it is 33.19 W m⁻² on average, as expected. This also means that averaging these results over 60 min can lead to erroneous values for the irradiance, especially under broken clouds, since the periods of cloud enhancement within each 1 h window will not be taken into account.

The GHI retrieved from the 3D MYSTIC LUT shows significantly lower bias in the case of the non-linear temperature model (2.00 W m⁻² compared to 12.71 W m⁻² using the DISORT COD method, averaged over both campaigns), but in the case of the linear temperature model the bias is higher (25.30 W m⁻² compared to 16.09 W m⁻²). The good performance of the non-linear model could indeed be a result of the improved treatment of the PV module temperature during broken-cloud conditions: although neither model contains a dynamic term, the non-linear model couples irradiance and wind speed in a more physically correct way. However, it must be noted that due to the restrictions on the MYSTIC LUT (see Table 5), the number of inferred irradiances included in the statistics is far lower than for the COD method. Another confounding factor could be that the PV panels show better efficiency at higher irradiance, and since the MYSTIC method also takes overshoots into account, the irradiance is on average higher: 603.78 W m⁻² compared to 261.45 W m⁻² for the COD method. The larger fluctuations in irradiance under broken-cloud conditions also lead to a larger RMSE than for the case of cloudy skies, and in summer the RMSE is the highest, as expected.

Figure 8 shows the GHI retrievals from the AÜW systems with 15 min PV power measurements, under clear (top row) and cloudy (bottom row) skies. In this case the MYSTIC 3D LUT is not used, since the determination of the cloud fraction with coarsely resolved data leads to erroneous results, and the rapid fluctuations in irradiance under broken clouds are not properly captured at 15 min resolution. The DISORT 1D LUT performs well under clear skies, as to be expected, with the linear model showing an average MBE of 5.88 W m⁻² and the non-linear model 4.10 W m⁻². Once again the non-linear model outperforms the linear one in autumn 2018, but this trend is reversed in summer 2019. This systematic effect is most probably due to uncertainties in the temperature model calibration. Under cloudy skies the GHI retrievals show a significant positive bias of on average 41.87 W m⁻² (42.60 W m⁻²) for the linear (non-linear) model, which means that the retrieved COD is too small. This is discussed further in Sect. 4.4. Interestingly the large bias errors in tilted irradiance for the non-linear model are not evident in the horizontal irradiance results, which is probably due to the fact that far less points are taken into account in the statistics for GHI (compare the values of n in Figs. 8 and 6). In other words, the outliers have been removed by selecting either clear-sky days or periods for which the cloud fraction is 100 %.

Figure 8 Combined comparison of GHI retrieved from PV power measurements with that measured by horizontal pyranometers, using data in 15 min resolution from MS02, PV01, PV06, PV08, PV10, PV11, PV14, PV16, PV17, PV18, and PV21, for 2018 (left two columns) and 2019 (right two columns), and both the linear and the non-linear temperature models. The top row shows GHI retrieved via AOD under clear skies using the DISORT AOD LUT, and the bottom row is for cloudy periods via the COD using the DISORT COD LUT.

4.3.2 Comparison to satellite and weather model irradiance data

One of the main aims of this work is to show that PV systems can provide a reliable source of information on global horizontal irradiance that is complementary to that from satellite and weather models. Figures 9 and 10 (with corresponding statistical measures in Table 11) show the comparison between GHI retrieved from PV power using the aerosol and cloud optical depths and that from the CAMS retrieval, for 1 and 15 min power data, respectively. Under clear-sky conditions the retrieved GHI shows an average MBE of 15.95 W m⁻² for the linear model and 7.54 W m⁻² for the non-linear model. These values are similar to those found by comparing with ground-based pyranometers, confirming the accuracy of the CAMS data. On the other hand, the GHI retrieved using the DISORT COD LUT under cloudy skies shows a significant negative bias compared to the CAMS retrieval (−37.77 W m⁻² in autumn and −86.87 W m⁻² in summer). There are two possible reasons for this: firstly the simplification to one cloud type means that thinner clouds or multi-layer cloud situations are not properly represented in the model (see the discussion on COD retrievals in Sect. 4.4 below). However, the main reason is related to the retrieval algorithm: by only considering measurements where the PV system is under a cloud, only the periods with lower irradiance values are retrieved, and overshoots are ignored so that at 1 min resolution a large negative bias in irradiance is found. Since the CAMS irradiance retrieval is based on the Heliosat-4 method, in which cloud properties from the APOLLO_NG method are updated every 15 min (Qu et al., 2017), one should expect this bias to reduce at coarser resolution. Indeed, the 15 min data in the bottom row of Fig. 10 confirm this: the average bias is reduced to −3.73 W m⁻² in autumn and −15.07 W m⁻² in summer. Note that an averaging to 60 min is not performed due to the limitations of the DISORT COD algorithm, as discussed in Sect. 4.3.1.

Figure 9 Combined comparison of GHI retrieved from PV power measurements with that from CAMS, using data in 1 min resolution from MS02, PV11, PV12, PV15, and PV19, for 2018 (left two columns) and 2019 (right two columns), and both the linear and the non-linear temperature models. The top row shows GHI retrieved via AOD under clear skies using the DISORT AOD LUT, and the bottom row is for cloudy periods via the COD using the DISORT COD LUT.

Figure 10 Combined comparison of GHI retrieved from PV power measurements with CAMS, using data in 15 min resolution from MS02, PV01, PV06, PV08, PV10, PV11, PV14, PV16, PV17, PV18, and PV21, for 2018 (left two columns) and 2019 (right two columns), and both the linear and the non-linear temperature models. The top row shows GHI retrieved via AOD under clear skies using the DISORT AOD LUT, and the bottom row is for cloudy periods via the COD using the DISORT COD LUT.

Table 11 The mean bias error and root mean squared error (in W m⁻²) along with the rMBE and rRMSE (in brackets in %) of GHI retrievals using 1D DISORT (AOD and COD) and the 3D MYSTIC LUT compared to CAMS, for 1 and 15 min data.

The comparison with COSMO model data is shown in Fig. 11 and Table 12, where in this case the data are averaged over a 60 min period. It is evident that the COSMO model shows a bias under clear-sky conditions: here the assumed AOD is too high so that the irradiance turns out to be too small, with an average bias of 60.92 W m⁻². On the other hand, under cloudy conditions and especially under low-light conditions in summer, the irradiance from COSMO is larger than that retrieved from PV plants, which means that the COD in COSMO is too small. Here the average bias is −38.36 W m⁻². These results confirm the findings of Frank et al. (2018) and Zubler et al. (2011) and are discussed further in connection with the cloud optical depth in Sect. 4.4 below.

Figure 11 Combined comparison of GHI retrieved from 60 min averaged PV power measurements under clear (top row) or completely cloudy (bottom row) conditions using the optical depth via DISORT LUT with that from the COSMO model, for all stations and for 2018 (left two columns) and 2019 (right two columns), using both the linear and the non-linear temperature models.

Table 12 The mean bias error and root mean squared error (in W m⁻²) along with the rMBE and rRMSE (in brackets in %) of 60 min averaged GHI retrievals using 1D DISORT (AOD and COD) and the 3D MYSTIC LUT compared to COSMO model data.

4.4 Cloud optical depth retrievals

As discussed in Sect. 2.3.3, the cloud optical depth is retrieved from both PV systems and pyranometers, using a DISORT-based LUT. In order to avoid errors due to 3D effects, in this work only data with a cloud fraction of 1 are considered, in other words only completely overcast conditions. The results for the linear temperature model are shown in Figs. 12 and 13, compared to the APOLLO_NG and COSMO data, respectively. As can also be seen in Tables 13 and 14, in most cases a smaller COD is extracted from PV systems, except for the comparison between the summer campaign and COSMO data, in which case a positive mean bias of approximately COD =8 is found. Overall, the COD is mostly in the range between 1 and 10, for both campaigns. Taken at face value, the negative bias with respect to APOLLO_NG would imply a positive bias in GHI with respect to CAMS, which is not seen in the 1 and 15 min retrievals. However these results cannot be directly compared due to the effect of both spatial and temporal averaging as well as the limitation of the DISORT COD LUT algorithm, which ignores 3D effects.

Figure 12 Combined comparison of COD retrieved from PV power measurements under completely cloudy conditions using the DISORT LUT with that from APOLLO_NG, for all stations and for 2018 (a) and 2019 (b), using the linear temperature model and averaged over 60 min.

Figure 13 The same as Fig. 12 but compared to the total COD from the COSMO model.

Table 13 The mean bias error and root mean squared error (in W m⁻²) along with the rMBE and rRMSE (in brackets in  %) of COD retrievals from PV systems compared to the APOLLO_NG data.

Table 14 The mean bias error and root mean squared error (in W m⁻²) along with the rMBE and rRMSE (in brackets in %) of COD retrievals from PV systems compared to the COSMO model predictions.

Figures 14 and 15 show the same results using measured pyranometer data to infer the COD. These retrievals show a trend similar to the PV-based ones: once again it is evident that the COD is mostly below 10, and in this range the retrieved data have a positive relative bias. There are several possible reasons for this: firstly it is evident from Eq. (11) that the retrieval is more sensitive to errors in inverted irradiance (transmission) for smaller COD. On the other hand, it must also be noted that the efficiency of both PV modules and silicon-based pyranometers shows a logarithmic dependence on irradiance so that any inaccuracies in the PV model parameters or the pyranometer calibration will have a larger effect on the inverted irradiance under low-light conditions (higher COD). In addition, since the LUT is constructed with water clouds, the effect of optically thin ice clouds is not properly taken into account. Since ice particles scatter slightly less in the forward direction, $\mathrm{1}-g\simeq \mathrm{0.3}$ is larger than for water clouds ($\mathrm{1}-g\simeq \mathrm{0.15}$), and thus a smaller optical depth could lead to similar irradiance at the ground (see Eq. 11). Thirdly, since clouds become more absorbing in the near infrared and considering that silicon PV is sensitive to wavelengths up to 1200 nm, spectral effects could also lead to a bias in the results. In general it must be said that even with measurements at two different wavelengths there exists an ambiguity in the determination of effective radius and COD (Nakajima and King, 1990) so that in the case of spectrally integrated PV-inverted irradiance one cannot expect to have enough information to accurately determine cloud optical properties in all situations.

Figure 14 Combined comparison of the COD retrieved from tilted (plane-of-array – poa) pyranometer measurements under completely cloudy conditions using the DISORT LUT with that from APOLLO_NG, for all stations and for 2018 (a) and 2019 (b), using the linear temperature model and averaged over 60 min.

Figure 15 The same as Fig. 14 but compared to the total COD from the COSMO model.

Notwithstanding the bias in COD retrievals, the ground-based method presented here can still complement satellite retrievals, in particular due to the potentially higher spatiotemporal resolution achievable with large amounts of high-frequency data spread over a large area. Once again, for the summer months the COSMO data show a large mean bias error of COD =7.69, even for large values of COD, confirming the findings of Frank et al. (2018).

5 Conclusions and outlook

In summary, in this work a framework for extracting both global tilted and horizontal irradiance from PV power data has been presented, and a first test for retrieving cloud optical depth is carried out. The algorithm makes use of state-of-the-art radiative transfer solvers in libRadtran, in conjunction with different sources of data for the state of the atmosphere, in particular the aerosol and water vapour content. The calibration procedure uses an explicit calculation of diffuse and direct irradiance, taking into account the spectral response of the relevant PV technology and the optical properties of the glass surface. Where necessary the module temperature is modelled using weather model data for ambient temperature and wind speed input. The PV systems are calibrated using a libRadtran clear-sky simulation with the DISORT solver, with inputs from the COSMO model and AERONET, and the algorithm can be adapted to each system and situation, depending on which parameters are known and which need to be determined by non-linear optimisation.

Once calibrated, the measured PV data are used to extract global tilted irradiance under all sky conditions. In order to take into account the spectral mismatch between the spectral response of PV modules and the entire broadband spectrum, a situation-dependent fit for the dependence of this mismatch on atmospheric conditions is performed, using simulated data for clear-sky conditions and the water vapour and aerosol load of each site. The GTI is then compared to that measured by tilted pyranometers: the retrieved GTI at 1 min (15 min) resolution has a mean bias error of 18.68 W m⁻² (34.47 W m⁻²) averaged over the two measurement campaigns, using the linear temperature model. The non-linear Faiman temperature model achieves a mean bias error of 8.77 W m⁻² for the systems with 1 min data, but for those with 15 min data the calibration algorithm fails to find an optimal solution in several cases so that in the end the mean bias error is 91.98 W m⁻². This shows that an accurate calibration is essential in order to accurately extract irradiance.

The inverted GTI is used to find the global horizontal irradiance using three different methods: (i) under persisting clear or (ii) cloudy conditions a lookup table based on a 1D DISORT simulation is used in order to find either the AOD or the COD and thus the global horizontal irradiance; (iii) under broken-cloud conditions a LUT based on 3D MYSTIC simulations is used to translate the tilted irradiance to horizontal irradiance, using the geometry of the system, sun position, and cloud fraction as inputs. The retrieved GHI is then compared to pyranometer measurements: in the case of the 1D LUT method, with 1 min data under clear (cloudy) skies, the mean bias error is 13.79 W m⁻² (16.09 W m⁻²) for the linear model and 5.79 W m⁻² (12.71 W m⁻²) for the non-linear temperature model. Comparison of the 15 min GHI retrievals with CAMS data reveals a positive bias under clear skies of 3.53 W m⁻² (2.27 W m⁻²) for the linear (non-linear) model, whereas under cloudy skies there is a negative bias of −9.90 and −8.91 W m⁻², respectively. Considering the difference between point measurements and satellite pixels, as well as the inherent bias from considering only periods with 100 % cloud fraction for the DISORT method, these results must be interpreted with care. In the case of clear skies, small errors in the temperature model can lead to large errors in extracted irradiance, and in the case of cloudy skies, simplifying assumptions on the cloud type can lead to errors in COD and extracted irradiance.

The retrieved GHI shows a large bias when comparing it with COSMO model data, thus confirming the results of Frank et al. (2018) and Zubler et al. (2011): under clear skies the 60 min averaged GHI has a mean bias error of 60.92 W m⁻², since COSMO has in general a too large aerosol load, whereas under cloudy skies in summer the MBE is −38.36 W m⁻², since under cloudy conditions the COSMO model generally tends to overestimate the irradiance. The latter result is also confirmed by the COD retrievals: in summer there is a positive bias of 7.69 for the retrieved COD relative to COSMO – the COD in the weather model is thus on average too small.

Overall, the largest source of error in the model chain comes from the PV model itself – an accurate calibration is vital in order to be able to extract irradiance reliably. In general the non-linear temperature model performs the best with high-frequency PV data. In this regard it is helpful to use measured module temperature rather than relying on temperature models. More accurate results could also be achieved by using PV current measurements, since in this case the temperature dependence is almost negligible. This will be explored in future work.

The DISORT LUT is only employed during periods of persistent cloudy or clear-sky conditions in order to infer COD or AOD, respectively, whereas the effect of 3D transport of photons is only taken into account with the MYSTIC LUT for the GHI. This means that the algorithm in its present form can only extract direct and diffuse components reliably under stable conditions for which 1D radiative transfer is still a reasonable approximation. A possible extension to this will be studied in future work, in which explicit 3D simulated cloud fields will be used in conjunction with the 1D DISORT simulation in order to quantify the error that results from neglecting 3D radiative transfer. It is only once the three dimensionality of atmospheric radiative transfer and additional information on cloud type are taken into account that one can accurately extract the direct and diffuse irradiance components under broken-cloud conditions.

Another aspect not taken into account in this study is the possible gain in using several different PV systems that are close enough to each other so as to be able to see the same or similar portions of the sky. In this case it is conceivable that one could extract more information about cloud properties and irradiance compared to that obtained from just one system. Indeed, the rapid proliferation of PV installations could make such multi-sensor inversions an interesting future prospect.

The ultimate goal of this work is to show that PV power data can be used to infer global horizontal irradiance and optical properties of the atmosphere, and the algorithm presented here is the first step in this direction. Although it is clear from the above that this is feasible, moving towards operational use would require several further steps. The largest source of bias comes from the calibration step, in particular the effect of temperature. Access to direct current (DC) data on the inverter level would allow for a much more accurate extraction of the irradiance, without the confounding effect of temperature, since the dependence of MPP current on temperature is an order of magnitude smaller than that of MPP voltage and power. Since inverter data are commonly available in industry, this should be possible provided one has the legal right to access the data. Additionally, the data pre-processing needs to be automated. For instance, one needs to exclude PV systems and/or data subsets that do not meet certain criteria such as shading of PV modules or inverter clipping, and one needs to be able to detect clear-sky periods automatically, even if the system orientation is unknown. This could be done by developing a hybrid approach using both physical modelling and artificial intelligence (AI) algorithms for pattern recognition. In addition, the calibration procedure itself could be augmented with AI to find the unknown parameters more effectively.

Once these aspects are overcome and appropriate agreements with industry partners are made for access to the data, there should be nothing standing in the way of operational use. If this could be achieved at a large scale it would allow for a better characterisation of solar irradiance at the ground and open up several possibilities for improving PV power forecasts at different time horizons. At shorter timescales (sub-hourly) one could use the additional information on irradiance variability as further input to empirical forecasts based on statistical methods, whereas at longer timescales (more than 3 h) these data could be assimilated into weather models. In order for this to make a difference one needs a much larger dataset of PV systems for the analysis, which then requires further automation. The first steps in this direction are currently being explored.

Appendix A: Modelling details

A1 Parametric power model

In Buchmann (2018) a simple model is proposed, based on the combination of several different modelling steps from the literature (see for instance Skoplaki and Palyvos, 2009, or Dows and Gough, 1995). Here the PV power output is written as

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& {\mathit{P}}_{\mathrm{AC},\mathrm{mod}}\\ & =A{\mathit{\eta }}_{\mathrm{DCAC}}{\mathit{\eta }}_{\mathrm{module}}\left(\mathit{T},{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }},{\mathit{v}}_{\mathrm{wind}}\right){\mathit{G}}_{\mathrm{tot},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}\\ & =A{\mathit{\eta }}_{\mathrm{DCAC}}{\mathit{\eta }}_{\mathrm{module}}\left(\mathit{T},{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }},{\mathit{v}}_{\mathrm{wind}}\right)\\ & [{\mathit{G}}_{\mathrm{dir},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}+{\mathit{G}}_{\mathrm{diff},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}],\end{array}\end{array}$$

where A is the surface area of the PV system; η_DCAC is the converter efficiency; and the direct and diffuse components of the irradiance in the plane of array and underneath the glass covering (see Sect. A2) are given by ${\mathit{G}}_{\mathrm{dir},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}$ and ${\mathit{G}}_{\mathrm{diff},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}$, respectively. The temperature-dependent module efficiency is defined by (Evans and Florschuetz, 1977)

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& {\mathit{\eta }}_{\mathrm{module}}\left(\mathit{T},{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }},{\mathit{v}}_{\mathrm{wind}}\right)\\ & ={\mathit{\eta }}_{\mathrm{module},\mathrm{n}}[\mathrm{1}-\mathit{\zeta }\left({\mathit{T}}_{\mathrm{module}}-T_{\mathrm{n}}\right)],\end{array}\end{array}$$

where η_module,n and ζ are the module efficiency and temperature coefficient under standard test conditions (STC); i.e. at T_n=25 ^∘C, $G_{\mathrm{tot},\mathrm{SW},\mathrm{n}}^{\mathrm{\angle }}=\mathrm{1000}$ W m⁻² with an air mass of 1.5. Note that in principle one could include a logarithmic dependence of the module efficiency on irradiance (Sauer, 1994), which is however not considered here. The module temperature is modelled using both the linear model (TamizhMani et al., 2003) defined by

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{\mathit{T}}_{\mathrm{module}}& =u_{\mathrm{0}}{\mathit{T}}_{\mathrm{ambient}}+u_{\mathrm{1}}{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}+u_{\mathrm{2}}{\mathit{v}}_{\mathrm{wind}}\\ & +u_{\mathrm{3}}{\mathit{T}}_{\mathrm{sky}}\end{array}\end{array}$$

and the non-linear model (Barry et al., 2020; Faiman, 2008) defined by

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\mathit{T}}_{\mathrm{module}}& ={\mathit{T}}_{\mathrm{ambient}}+{\displaystyle \frac{{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}}{u_{\mathrm{1}}+u_{\mathrm{2}}{\mathit{v}}_{\mathrm{wind}}}}\\ & +u_{\mathrm{3}}\left({\mathit{T}}_{\mathrm{sky}}-{\mathit{T}}_{\mathrm{ambient}}\right),\end{array}\end{array}$$

where v_wind is the wind speed at 10 m above the ground, ${\mathit{T}}_{\mathrm{sky}}^{\mathrm{4}}={\mathit{G}}_{\mathrm{LW}}^{\downarrow }/\left(\mathit{ϵ}\mathit{\sigma }\right)$ defines the sky temperature, and $u_{\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}}$ are model parameters. Here an emissivity of ϵ=1 is assumed, and σ is the Stefan–Boltzmann constant. Note that Eqs. (A1), (A2), (A3), and (A4) can be combined into the general non-linear expressions given in Eqs. (1) and (2) (see for instance Skoplaki and Palyvos, 2009; Dows and Gough, 1995), where in this special case the coefficients are given by

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& b_{\mathrm{1}}=s(\mathrm{1}+\mathit{\zeta }\mathrm{25}),b_{\mathrm{2}}=-u_{\mathrm{1}}s\mathit{\zeta },b_{\mathrm{3}}=-u_{\mathrm{0}}s\mathit{\zeta },\\ & b_{\mathrm{4}}=-u_{\mathrm{2}}s\mathit{\zeta },b_{\mathrm{5}}=-s\mathit{\zeta }{u}_{\mathrm{3}}\end{array}\end{array}$$

for the linear and

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& b_{\mathrm{1}}^{\prime }=s(\mathrm{1}+\mathit{\zeta }\mathrm{25}),b_{\mathrm{2}}^{\prime }=-{\displaystyle \frac{u_{\mathrm{1}}}{s\mathit{\zeta }}},b_{\mathrm{3}}^{\prime }=s\mathit{\zeta }(u_{\mathrm{3}}-\mathrm{1}),\\ & b_{\mathrm{4}}^{\prime }=-{\displaystyle \frac{u_{\mathrm{2}}}{s\mathit{\zeta }}},b_{\mathrm{5}}^{\prime }=-s\mathit{\zeta }{u}_{\mathrm{3}}\end{array}\end{array}$$

for the non-linear temperature model, where $s\equiv A{\mathit{\eta }}_{\mathrm{DCAC}}{\mathit{\eta }}_{\mathrm{module},\mathrm{STC}}$ is a constant scaling factor. The model equations in Eqs. (1) and (2) are used in this work, and the coefficients $u_{\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}}$, ζ, and s are allowed to vary freely, unless their a priori values are known from datasheets and/or from temperature modelling. In the cases where measured temperature is available, the parameters $u_{\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}}$ fall away.

A2 Optical model of glass covering

A2.1 Model formulation

In order to model the optics of the glass surface of the PV modules the following equation for the transmission of photons as a function of the incident angle Θ is used (De Soto et al., 2006; Sjerps-Koomen et al., 1996):

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{\mathit{\tau }}_{\mathrm{PV}}\left(\mathbf{\Theta }\right)& =\exp \left(-{\displaystyle \frac{\mathit{\kappa }L}{\cos {\mathbf{\Theta }}_{\mathrm{r}}}}\right)\\ & \left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{{\sin }^{\mathrm{2}}\left({\mathbf{\Theta }}_{\mathrm{r}}-\mathbf{\Theta }\right)}{{\sin }^{\mathrm{2}}\left({\mathbf{\Theta }}_{\mathrm{r}}+\mathbf{\Theta }\right)}}+{\displaystyle \frac{{\tan }^{\mathrm{2}}\left({\mathbf{\Theta }}_{\mathrm{r}}-\mathbf{\Theta }\right)}{{\tan }^{\mathrm{2}}\left({\mathbf{\Theta }}_{\mathrm{r}}+\mathbf{\Theta }\right)}}\right]\right),\end{array}\end{array}$$

where Θ_r is the angle of refraction from Snell's law (nsin Θ_r=sin Θ), n is the index of refraction of glass, κ is the glazing extinction coefficient, and L is the glazing thickness. The incident angle Θ is the angle between the solar position vector and normal vector of the PV array, defined by

$$\begin{array}{}\text{(A8)}& \cos \mathbf{\Theta }=\cos {\mathit{\theta }}_{\mathbf{0}}\cos \mathit{\theta }+\sin {\mathit{\theta }}_{\mathbf{0}}\sin \mathit{\theta }\cos \left({\mathit{\varphi }}_{\mathbf{0}}-\mathit{\varphi }\right),\end{array}$$

where θ is the angle of inclination of the PV array, ϕ is its orientation, θ₀ is the solar zenith angle, and ϕ₀ is the solar azimuth.

In principle one should take into account the wavelength dependence of n and κ; however for most materials they are relatively constant, with n increasing slightly at lower wavelengths, and since in practice the exact properties of the glass covering for each system are unknown, it suffices to use the effective values for all wavelengths.

The so-called incidence angle modifier (see also Duffie and Beckman, 2013) is defined by the ratio of the transmission τ_PV(Θ) and the transmission at normal incidence; i.e.,

$$\begin{array}{}\text{(A9)}& {\mathit{\tau }}_{\mathrm{PV},\mathrm{rel}}\left(\mathbf{\Theta }\right)\equiv {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{PV}}\left(\mathbf{\Theta }\right)}{{\mathit{\tau }}_{\mathrm{PV}}\left(\mathrm{0}\right)}},\end{array}$$

where

$$\begin{array}{}\text{(A10)}& {\mathit{\tau }}_{\mathrm{PV}}\left(\mathrm{0}\right)\begin{array}{l}=e^{-\mathit{\kappa }L}\left[\mathrm{1}-{\left(\frac{n-\mathrm{1}}{n+\mathrm{1}}\right)}^{\mathrm{2}}\right]\\ =e^{-\mathit{\kappa }L}\frac{\mathrm{4}n}{(n+\mathrm{1}{)}^{\mathrm{2}}}.\end{array}\end{array}$$

The use of a relative transmission coefficient is justified by the fact that the absolute transmittance is already taken into account when characterising the solar cell under standard conditions (Sjerps-Koomen et al., 1996). The normalisation with τ_PV(0) means that the result is less sensitive to the product κL and more to the angle of incidence. In the forward model the wavelength-integrated direct irradiance and the diffuse radiance beams are each multiplied with the factor τ_PV,rel(Θ) in order to take into account the attenuation due to the glass surface. The values of the extinction coefficient and thickness of the glass are fixed to κ=4 m⁻¹ and L=0.002 m, respectively (De Soto et al., 2006), whereas the effective refractive index n is allowed to vary (n=1.526 with an a priori error of 1 %). In principle one could also vary the product κL that controls the absorption, but as mentioned above the incident angle modifier approach is applicable to a wide range of glass covers (Duffie and Beckman, 2013; Klein, 1979).

A2.2 Optical model in forward model calibration

In the forward model calibration, the transmission function can then be used to calculate the direct and diffuse irradiance ${\mathit{G}}_{\mathrm{dir},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}$ and ${\mathit{G}}_{\mathrm{diff},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}$ as

$$\begin{array}{}\text{(A11)}& {\mathit{G}}_{\mathrm{dir},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}={\displaystyle \frac{\cos \mathbf{\Theta }}{\cos {\mathit{\theta }}_{\mathbf{0}}}}{\mathit{G}}_{\mathrm{dir}}^{\downarrow }{\mathit{\tau }}_{\mathrm{PV},\mathrm{rel}}\left(\mathbf{\Theta }\right)\end{array}$$

and

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{\mathit{G}}_{\mathrm{diff},\mathrm{PV},\mathit{\tau }}^{\mathrm{\angle }}& =\underset{\mathit{\varphi }}{\overset{\mathrm{2}\mathit{\pi }+\mathit{\varphi }}{\int }}\underset{-\mathit{\theta }}{\overset{\mathit{\pi }/\mathrm{2}-\mathit{\theta }}{\int }}{\mathbf{L}}_{\mathrm{diff}}^{\mathrm{\Omega }}\cos {\mathbf{\Theta }}^{\mathbf{\prime }}{\mathit{\tau }}_{\mathrm{PV},\mathrm{rel}}\left({\mathbf{\Theta }}^{\mathbf{\prime }}\right)d\mathrm{\Omega },\\ & \text{for}\cos {\mathbf{\Theta }}^{\mathbf{\prime }}\ge \mathrm{0},\end{array}\end{array}$$

where ${\mathbf{L}}_{\mathrm{diff}}^{\mathrm{\Omega }}$ is the diffuse radiance distribution calculated by DISORT. In this way there is no need for an empirical incidence angle modifier, since the direction of each diffuse photon is explicitly described. The same formulation is used to calculate ${\mathit{G}}_{\mathrm{dir},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}$ and ${\mathit{G}}_{\mathrm{diff},\mathrm{SW},\mathit{\tau }}^{\mathrm{\angle }}$, i.e. over the whole wavelength range.

A2.3 Inversion of the optical model

For the inversion of the PV model, two different methods are used: for the extraction of cloud optical depth with DISORT, the optical model can be explicitly taken into account in the radiative transfer simulation, whereas for the direct extraction of GTI and its translation to GHI with the MYSTIC lookup table, the empirical formulation (Duffie and Beckman, 2013) for the effective angle for diffuse photons as a function of tilt angle θ,

$$\begin{array}{}\text{(A13)}& {\mathrm{\Theta }}_{\mathrm{diff}}=\mathrm{59.7}-\mathrm{0.1388}\mathit{\theta }+\mathrm{0.001497}{\mathit{\theta }}^{\mathrm{2}},\end{array}$$

is used for all time points with a clearness index below 0.3 so that the final inverted GTI is given by

$$\begin{array}{}\text{(A14)}& {\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathrm{inv}}^{\mathrm{\angle }}={\displaystyle \frac{{\mathit{G}}_{\mathrm{tot},\mathrm{SW},\mathrm{inv},\mathit{\tau }}^{\mathrm{\angle }}}{{\mathit{\tau }}_{\mathrm{PV},\mathrm{rel},\mathrm{eff}}}},\end{array}$$

where

$$\begin{array}{}\text{(A15)}& {\mathit{\tau }}_{\mathrm{PV},\mathrm{rel},\mathrm{eff}}=\left\{\begin{array}{lcc}{\mathit{\tau }}_{\mathrm{PV},\mathrm{rel}}\left({\mathbf{\Theta }}_{\mathrm{dir}}\right)& \text{if}& {\mathit{k}}_{\mathit{i}}\ge \mathrm{0.3}\\ {\mathit{\tau }}_{\mathrm{PV},\mathrm{rel}}\left({\mathrm{\Theta }}_{\mathrm{diff}}\right)& \text{if}& {\mathit{k}}_{\mathit{i}}<\mathrm{0.3}\end{array}\right.\end{array}$$

is the effective incidence angle modifier for clearness index k_i.

A3 Spectral mismatch fitting procedure

The ratio of silicon PV irradiance to broadband irradiance as defined in Eq. (5) is a function of the atmospheric composition (primarily water vapour content) and angle of incidence of the incoming solar beam. Indeed, the shape of the diurnal variation in ξ_spec,GTI is dependent on the azimuth angle of the PV plant, as shown in Figs. A1 and A2, where the ratio is plotted for both GHI and GTI, for two different PV systems with different orientations, using libRadtran clear-sky simulations from the different measurement campaigns. The points in the upper (lower) plots in each figure are coloured according to precipitable water (AOD), and it is evident that for a given incident angle the water vapour plays the largest role in determining the ratio.

The mismatch ratio remains at roughly 0.83 throughout the day, depending on water vapour column. In the case of PV12 in Fig. A1, one can see that in the mornings the ratio is smaller, i.e. more red light, whereas in the evenings the ratio is larger, since the panel is facing more to the south-east. This means that the panel detects more diffuse light in the evenings. The behaviour at PV15 (Fig. A1) is the opposite – in the morning the ratio is larger, since the panel looks more to the south-west (mornings have a higher diffuse component, i.e. more blue light), whereas in the evenings it is smaller. However in the summer time when the sun goes down far to the north, the diffuse component again plays a role so that the ratio becomes larger (right-hand plot). In general the variation is the greatest for larger SZAs, as seen in the greyed-out areas in the plots.

Figure A1 The spectral mismatch ratio at PV12 (azimuth angle ϕ≃133^∘) as a function of time (mean incident angle) as well as the water vapour column (a, b) and AOD (c, d), for both 2018 (a, c) and 2019 (b, d) measurement campaigns. Θ_min is the minimum incident angle. The solid and dashed black lines represent the mean mismatch ratio for GTI and GHI, respectively. The fit for GTI corresponding to Eq. A16 is shown as a dashed red line and is performed up to an SZA of 80^∘.

Figure A2 The same as Fig. A1 but for PV15 (azimuth angle ϕ≃228^∘).

As shown in Figs. A1 and A2, the data are first grouped by the time of the day in order to calculate the mean value, then split into two halves on either side of the minimum incident angle (Θ_min). Each branch can then be fitted with the function

$$\begin{array}{}\text{(A16)}& \langle {\mathit{\xi }}_{\mathrm{spec},\mathrm{GTI}}\rangle =x_{\mathrm{0}}\exp \left(-{\displaystyle \frac{x_{\mathrm{1}}}{\cos \mathbf{\Theta }}}-{\displaystyle \frac{x_{\mathrm{2}}}{\cos {\mathbf{\Theta }}^{\mathrm{2}}}}\right),\end{array}$$

shown as the dashed red line in the plots. This shows the general form of the diurnal variation in spectral mismatch and that for silicon PV the ratio is about 83 % for most of the day.

In order to capture the effect of precipitable water and AOD, the function

$$\begin{array}{}\text{(A17)}& {\mathit{\xi }}_{\mathrm{spec},\mathrm{GTI}}=p_{\mathrm{0}}\exp \left(-p_{\mathrm{1}}\mathbf{\left[}{\mathit{H}}_{\mathbf{2}}\mathit{O}\mathbf{\right]}-p_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{a}}\right)\end{array}$$

is fitted to the irradiance data for each time step (averaged over all days of the respective measurement campaign), where [H₂O] is the precipitable water from COSMO and τ_a the AOD at 500 nm from AERONET. The fit coefficients are then interpolated over the entire dataset in order to calculate ξ_spec,GTI at any time of the day. Future work will examine the effects of clouds on the spectral mismatch factor – here the clear-sky fit is applied to all data, which due to whitening of the skylight by clouds could lead to a bias in the final result.