A New Model for Estimating the Diffuse Fraction of Solar Irradiance for Photovoltaic System Simulations

Abstract

We present a new model for the calculation of the diffuse fraction of the global solar irradiance for solar system simulations. The importance of an accurate estimation of the horizontal diffuse irradiance is highlighted by findings that an inaccurately calculated diffuse irradiance can lead to significant over- or underestimations in the annual energy yield of a photovoltaic (PV) system by as much as 8%. Our model utilizes a time series of global irradiance in one-minute resolution and geographical information as input. The model is validated by measurement data of 28 geographically and climatologically diverse locations worldwide with one year of one-minute data each, taken from the Baseline Surface Radiation Network (BSRN). We show that on average the mean absolute deviation of the modelled and the measured diffuse irradiance is reduced from about 12% to about 6% compared to three reference models. The maximum deviation is less than 20%. In more than 80% of the test cases, the deviation is smaller 10%. The root mean squared error (RMSE) of the calculated diffuse fractions is reduced by about 18%.

1. Introduction

Adapting the common terminology in energy meteorology to differentiate between the power and energy of the solar radiation, the word ‘irradiance’ is used in this work to denote the instantaneous solar power per square meter in W/m², whereas the word ‘irradiation’ refers to the integral of the irradiance over time, thus denoting the energy of the solar radiation in Ws/m² or kWh/m² [1].

In photovoltaic (PV) system simulations, the global horizontal irradiance and the ambient temperature are the two most important inputs in order to determine the PV system’s energy output. The global horizontal irradiance is split up in its direct and diffuse components. These components are then separately translated to a tilted plane if the PV system in question has a module orientation that differs from the horizontal plane. In simple terms the global irradiance incident on a tilted module is then calculated as the sum of the direct and the diffuse irradiance on the tilted plane.

The model for estimating the diffuse fraction of the global horizontal irradiance is hence the first element in a chain of models that is necessary to simulate the electrical output of a PV system. As such, it has strong influence on the final output of the simulation, which is demonstrated by the following comparative simulations for two locations: Lindenberg, Germany and Gobabeb, Namibia. The analysed PV system is a standard 8 kWp (kilo Watt peak) grid connected system, the simulation is conducted in one-minute resolution with measurement data from the Baseline Surface Radiation Network (BSRN) [2].

During the four exemplary days in June in Lindenberg, Germany, chosen for Figure 1, the model used for this comparison (reduced version of Reindl et al. [3]) underestimates the diffuse irradiance by 18%. In the next step of the simulation, the global irradiance on the tilted plane is calculated. In this case the modules are facing south and are elevated by 30° from the horizontal. The model applied for this step is from Hay and Davies [4]. When using the modelled horizontal diffuse irradiance, the resulting global irradiance on the modules is still 9% lower than using the measured horizontal diffuse irradiance.

The rest of the PV system model chain is then simulated with the help of the simulation core of PV software provider Valentin Software (Berlin, Germany) [5].

Table 1. Measured and modelled diffuse irradiation for Lindenberg, Germany (LIN), and Gobabeb, Namibia (GOB). The top section refers to the figures above, time ranges from 10–13 May 2003. The modelled underestimation of the diffuse irradiation by −17.9% leads to an underestimation of the global irradiation on the tilted PV module by −9.3%, hence leading to an underestimation of the simulated PV energy yield by 8.3%. When considering the whole year (2nd section), the modelled diffuse irradiation differs from the measured value by −7.2%, leading to an underestimation of the irradiation on module surface by −3.1%. The difference in the annual PV energy yield is −2.7%. In the 3rd and 4th section, the results of the same analysis are presented for Gobabeb, Namibia. Here, the PV module faces North and is tilted at 23°. The modelled sum of diffuse irradiation for selected days (23–26 July 2014) is 20.4% higher than the measured sum, leading to a deviation in the PV energy of 4.4%. Over the whole year of 2014, the deviation of the diffuse irradiation is even 42%, which causes a difference in the annual PV yield of 7.6%.

LIN, 10–13 May 2003	Unit	with Measured Data	Modelled	Deviation
Global horizontal irradiation	kWh/m²	19.7	-	-
Diffuse irradiation	kWh/m²	12.0	9.9	–17.9%
Global irradiation on tilted surface	kWh/m²	22.6	20.5	–9.3%
PV energy yield	kWh	65.6	60.2	–8.3%
LIN, whole year 2003	Unit	with Measured Data	Modelled	Deviation
Global horizontal irradiation	kWh/m²	1185.1	-	-
Diffuse irradiation	kWh/m²	555.9	515.7	–7.2%
Global irradiation on tilted surface	kWh/m²	1467.0	1422.0	–3.1%
PV energy yield	kWh	4339.2	4221.4	–2.7%
GOB, 23–26 July 2014	Unit	with Measured Data	Modelled	Deviation
Global horizontal irradiation	kWh/m²	18.2	-	-
Diffuse irradiation	kWh/m²	4.6	5.5	20.4%
Global irradiation on tilted surface	kWh/m²	22.5	23.7	5.0%
PV energy yield	kWh	66.2	69.1	4.4%
GOB, whole year 2014	Unit	with Measured Data	Modelled	Deviation
Global horizontal irradiation	kWh/m²	2433.1	-	-
Diffuse irradiation	kWh/m²	454.9	645.8	42.0%
Global irradiation on tilted surface	kWh/m²	2401.9	2600.8	8.3%
PV energy yield	kWh	6808.7	7325.4	7.6%

lists the results of the comparison.

During these four days, the total PV energy yield would be 65.6 kWh when using the measured horizontal diffuse irradiance values. With the diffuse irradiance modelled by Reindl et al. [3], the total PV energy yield is only 60.2 kWh—an underestimation of 8.3%. The comparison was also conducted for the whole year 2003 in Lindenberg, where the annual deviation of the modelled diffuse irradiance is −7.2%, leading to a deviation of the annual PV energy output of −2.7%.

The second half of Table 1 lists the results of the same comparison that was conducted for the location of Gobabeb in Namibia, for the year 2014. Here, the deviation of the annual diffuse irradiation is as high as 42% which leads to an overestimation of the global irradiation of the module surface of 8.3% and to an overestimation of the annual PV energy yield of 7.6%.

These examples highlight the importance of a more accurate estimation of the horizontal diffuse irradiance. An inaccurately calculated diffuse irradiance can lead to significant over- or underestimations in the annual energy yield of a PV system. This is especially relevant in the price sensitive market of PVs, where only few percent more or less of PV energy output can render a project possible or uneconomical [6].

2. Measurement Data and Methodology

In the following section the measurement data and the methodology used in this contribution are presented.

2.1. Data basis (Baseline Surface Radiation Network)

As a source of high quality measurement data the data base of the BSRN is used [2]. The BSRN comprises 59 stations worldwide, 44 of which provide one-minute measurements of global horizontal and diffuse horizontal irradiance. The time range of the measurements starts in 1992 for the first stations and is still running until now. For this study the following criteria were applied for selecting the datasets:

• High annual completeness of one-minute measurements of global and diffuse irradiance;
• Between 60° North and −60° South;
• No leap years.

Table 2. Overview over the 28 datasets that were used in this work. The locations are spread over the globe between −45° South and 52° North. Height above sea level, surface, topography and climate zones (according to Köppen [7]) show a high level of variation. The years of measurement were chosen to provide a high annual completeness of measurement (ACM), i.e., as few missing data points as possible. The total resulting amount of data points that is used in the following analysis exceeds 14.7 million (or approximately 7 million when omitting night time).

ID	Name	Country	Latitude in °	Longitude in °	Height in m	Time Zone	Surface	Topography	Climate Zone	ACM in %
ASP 2005	Alice Springs	Australia	–23.798	133.888	547	9.5	grass	flat, rural	BWh	99.4
BIL 2003	Billings	USA	36.605	–97.516	317	–6	grass	flat, rural	Cfa	99.6
BOU 2009	Boulder	USA	40.05	–105.007	1577	–7	grass	flat, rural	BSk	99.0
BRB 2010	Brasilia	Brasil	–15.601	–47.713	1023	–3	concrete	flat, rural	Aw	96.6
CAB 2009	Cabauw	Netherlands	51.971	4.927	0	1	grass	flat, rural	Cfb	99.1
CAM 2003	Camborne	UK	50.217	–5.317	88	0	grass	flat, rural	Cfb	90.4
CLH 2013	Chesapeake Light	USA	36.905	–75.713	37	–5	water, ocean	flat, rural	Cfa	99.8
CNR 2011	Cener	Spain	42.816	–1.601	471	1	asphalt	mountain valley, urban	Cfb	99.8
COC 2011	Cocos Islands	Cocos Islands	–12.193	96.835	–1	6.5	n.a.	n.a.	Af	95.6
DAA 2003	De Aar	South Africa	–30.667	23.993	1287	2	sand	flat, rural	BSk	88.1
DAR 2011	Darwin	Australia	–12.425	130.891	30	9.5	grass	flat, rural	Aw	100
FUA 2011	Fukuoka	Japan	33.582	130.375	3	9	asphalt	flat, urban	Cfa	99.9
GOB 2014	Gobabeb	Namibia	–23.561	15.042	407	1	n.a.	flat rural	BWh	100
IZA 2011	Izaña	Spain	28.309	–16.499	2372.9	0	rock	mountain top	Csb	96.1
LAU 2005	Lauder	New Zealand	–45.045	169.689	350	12	grass	flat, rural	Cfb	98.1
LER 2003	Lerwick	UK	60.133	–1.183	84	0	grass	hilly, rural	Cfb	100
LIN 2003	Lindenberg	Germany	52.21	14.122	125	1	cultivated	hilly, rural	Cfb	100
PAL 2011	Palaiseau	France	48.713	2.208	156	1	concrete	flat, urban	Cfb	99.7
PAY 2009	Payerne	Switzerland	46.815	6.944	491	1	cultivated	hilly, rural	Cfb	99.9
REG 2009	Regina	Canada	50.205	–104.713	578	–6	cultivated	flat, rural	BSk	100
SAP 2011	Sapporo	Japan	43.06	141.328	17.2	9	asphalt	flat, urban	Dfb	99.9
SBO 2009	Sede Boqer	Israel	30.905	34.782	500	2	desert rock	hilly, rural	Cwb	98.2
SMS 2007	São Martinho da Serra	Brasil	–29.443	–53.823	489	–3	concrete	flat, rural	Cfa	91.5
SOV 2001	Solar Village	Saudi Arabia	24.91	46.41	650	3	desert, sand	flat, rural	BWh	100
TAM 2006	Tamanrasset	Algeria	22.78	5.51	1385	1	desert, rock	flat, rural	BWh	99.9
TAT 2006	Tateno	Japan	36.05	140.133	25	9	grass	flat, urban	Cfa	99.9
TOR 2010	Toravere	Estonia	58.254	26.462	70	2	grass	flat, rural	Dfb	100
XIA 2009	Xianghe	China	39.754	116.962	32	8	desert, rock	flat, rural	Dwa	100

gives an overview of the locations and years that were used for validation. In total, 28 locations with one year of measurement each were selected. The datasets feature a high geographic and climatological diversity. The last column of the table lists the annual completeness of the measurements (ACM) in %. The validation datasets comprise more than seven million data points (nights omitted) on which the following analysis is based.

2.2. Description of Quantities and Models

For the calculation of the position of the sun, the solar position algorithm provided by the National Renewable Energy Laboratory (NREL; Golden, CO, USA) is used [8]. The clear sky irradiance $E_{\mathrm{clear}}$ is calculated on the basis of an adaption of the approach of Bourges [9]:

$$E_{\mathrm{clear}}=0.78E_{\mathrm{ext}}\sin {\left({\mathsf{\gamma }}_{\mathrm{S}}\right)}^{1.15}$$

(1)

where ${\mathsf{\gamma }}_{\mathrm{S}}$ is the elevation of the Sun and $E_{\mathrm{ext}}$ is the extra-terrestrial irradiance. The extra-terrestrial irradiance was calculated using Maxwell’s approach [10]. In 2014 this formula was identified as the most accurate for another subset of BSRN data by Hofmann et al. [11]. The clearness index $kt$ is the fraction of the measured global irradiance to the clear sky irradiance:

$$kt=\frac{E_{\mathrm{global},\mathrm{measured}}}{E_{\mathrm{clear}}}$$

(2)

In the models for calculating the diffuse fraction of the global irradiance that are presented in Section 2.3 a simpler approach of calculating the clear sky index is used:

$$E_{\mathrm{clear}}=E_{\mathrm{ext}}\sin \left({\mathsf{\gamma }}_{\mathrm{S}}\right)$$

(3)

This causes the $kt$ value at clear sky to be around 0.8 instead of 1 in the existing models. For the comparison of the results in Section 4, the calculation of the clearness index occurs according to the respective model description.

All PV system simulations are conducted using the simulation core of PV*SOL, a commercial PV system planning and simulation software by Valentin Software. More information about the models that are relied on in the simulation core can be found at PV*SOL [5].

2.3. Presentation of Existing Models

For the estimation of the diffuse fraction of the global horizontal irradiance, several algorithms were developed in the past. Most of them can be categorized as models with one or two parameters as input. The one-parameter models feature a simple dependency of the diffuse fraction ($df$) on the clearness index ($kt$), cf. Figure 2.

These models typically define three functions for different ranges of $kt$. These $kt$ ranges can be more or less referred to as different cloud situations. A $kt$ value of less than 0.4 means only 40% of the possible global irradiance is measured, which is a good indicator for overcast skies. The maximum $kt$ value that is detected at clear sky conditions is around 0.78–0.8. In between those areas, i.e., for $kt$ values between 0.4 and 0.78, broken cloud situations are most likely [1,12]. Values of $kt>1$ are possible due to broken cloud enhancement, firstly stated for the ultraviolet by Nack and Green [13] and later confirmed by Seckmeyer et al. [14]. Values of $kt>1$ are possible at all wavelengths of the solar spectrum.

The first model, a one-parameter approach, was presented by Liu and Jordan in 1960 [15], but it soon became apparent that it was not able to produce good results in other locations than it was designed for (Blue Hill, MA, USA) [16,17].

In consequence, other models were developed that can also be categorized as one-parameter models: Orgill and Hollands [18], Erbs, Klein and Duffie [19], Reindl, Beckman and Duffie [3] and Boland and Ridley [20]. A schematic overview of those models is provided in Figure 2, along with sample measurement data of Lindenberg, Germany, 2003, for reference. Other one-parameter approaches include the model by Oliveira et al. [21] that provides varying clearness index polynomials for three periods per year (December–January, April–August and September–March).

Another category of algorithms is formed by the two-parameter models that—in addition to the clearness index $kt$—also make use of the sun height, ${\mathsf{\gamma }}_{\mathrm{S}}$. Two-parameter models include the approaches by Reindl, Beckman and Duffie [3], Skartveit and Olseth [12] and Maxwell [10]. The reduced version of the two-parameter model of Reindl, Beckman and Duffie [3] is presented in Figure 3.

Not classifiable as one- or two-parameter model is the noteworthy approach by Furlan et al. [22] who developed a multi-parameter regression model for data from Sao Paolo, Brazil. Another important contribution was achieved by the model by Perez and Ineichen [23], which features a dynamic time-series approach to model the direct normal from the global irradiance based on the DISC model by Maxwell [10].

A good overview and an approach of global validation of the above mentioned models for calculating the diffuse fraction, also using BSRN data, is given by Zernikau [17]. In this thesis it was also shown that all analysed one- and two-parameter models showed relative mean absolute errors (rMAE) of (10.4 ± 0.4)% for the 24 BSRN locations that were included in the study. The author also analysed the minimal rMAE that can be achieved with a one- or two-parameter model by generating global medians of measurements of the diffuse fraction and the clearness index. According to this study, the minimal globally achievable rMAE for any two-parameter model is 8.9%.

Another one-location comparison of the models was conducted for Athens, Greece, by Kambezidis [24]. Similarly, a study comparing ten models was presented by Jacovides et al. for validation data of Athalassa, Cyprus [25]. A model-to-model comparison for Hong Kong without validation on measurement data is provided by Wong [26]. A comparison of eight models for the location of Vienna, Austria, was conducted by Dervishi [27], resulting in findings that are in general agreement to the above mentioned studies.

3. Presentation of a New Model for the Diffuse Fraction of Solar Irradiance

In order to reduce the uncertainty of PV system simulations, a new model for calculating the diffuse fraction of global horizontal irradiance is presented in this section. The model consists of three parts that are calculated independently and then combined depending on statistic features of the clearness index. Each part is presented and afterwards the combination of the three parts into a single resulting diffuse fraction $df$ is explained.

3.1. Part One. Diffuse Fraction as Function of Clearness Index

Like existing models with one parameter, this part of the new model makes use of the relation between the clearness index $kt$ and the diffuse fraction. Instead of parameterized functions, a matrix of probabilities is utilized. For the generation of the matrix, the one-minute time series of global and diffuse horizontal irradiance are converted into value pairs of the clearness index $kt$ and the diffuse fraction $df$, following the equations in Section 2.2. Each value pair is then stored into a matrix with $kt$ ranging from 0 to 1.5 and $df$ ranging from 0 to 1, both with a step size of 0.01. The frequency of occurrence of a specific $df$ value for a given $kt$ value is then converted into a probability value, so that for every value of $kt$ a function of cumulated probabilities can be calculated. Figure 4 shows an example of such a probability matrix. In the matrix shown here, measurement values from Alice Springs, Australia (2009), Billings, USA (2005), Boulder, USA (2010), Brasilia, Brazil (2011), Cabauw, The Netherlands (2011), Cener, Spain (2010), De Aar, South Africa (2003), Fukuoka, Japan (2013), Gobabeb, Namibia (2013), Lauder, New Zealand (2007), Lerwick, UK (2002), Lindenberg, Germany (2002), Payerne, Switzerland (2010), Regina, Canada (2011), Tateno, Japan (2003) and Xianghe, China (2006) were incorporated.

In order to determine a value for $df$ for a given $kt$, the procedure is as follows. Since this is the first part of the new model, the diffuse fraction of this part is referred to as $d{f}_1$:

(1)

Select column of probability matrix that corresponds to the $kt$ value;

(2)

Generate a Markov number $r_{\mathrm{M}}$ (randomized number between 0 and 1) [28];

(3)

Select the row where $r_{\mathrm{M}}$ is smaller than the cumulated probabilities for the first time;

(4)

The selected row corresponds to $d{f}_1$ value.

The usage of real measurement values, incorporated into a matrix of probabilities, holds the advantage of preserving the natural relationship of the diffuse fraction and the clearness index and additionally resulting in a more realistic variability of the modelled $df$ value series.

3.2. Part Two. Change of df as Function of Change of kt

By analyzing the extensive BSRN measurement database, a strong correlation has been found between the relative changes of the clearness index (from one minute to the next) to changes of the diffuse fraction. In Figure 5 this correlation is shown for Lindenberg, Germany, for the year 2003.

It was observed that for positive relative changes of $kt$ (when the current $kt$ is higher than the $kt$ one minute before), the diffuse fraction will most likely show a negative relative change. If the relative change of $kt$ is negative, the change of the diffuse fraction will be positive.

There are situations, however, where $kt$ is changing from one minute to the next without an observable change of $df$ (compare the horizontal value accumulation at $d_{df}=0$). These situations are typical for days with rapid irradiance enhancements due to moving broken clouds.

In correspondence to part 1, the relationship between the relative change of $df$ and the relative change of kt ($d_{kt}$) is also expressed in a matrix of probabilities, displayed in Figure 6. This matrix is only used in the $d_{kt}$ range of −0.5 to 1, since the amount of measurement values outside of this range is too small, which results in unwanted noise. The procedure to retrieve a value for the diffuse fraction in this part, $d{f}_2$, is as follows:

(1)

Calculate the relative change of $kt$ as:

$$d_{kt}=k{t}_{\mathrm{now}}/k{t}_{\mathrm{before}}-1$$

(4)

(2)

For $d_{kt}$ greater than −0.5 and smaller than 1

• Select the column of the probability matrix that corresponds to $d_{kt}$;
• Generate a Markov number $r_{\mathrm{M}}$ (randomized number between 0 and 1) [28];
• Select the row where $r_{\mathrm{M}}$ is smaller than the cumulated probabilities for the first time;
• The selected row corresponds to change of $df$, that is:

  $$d_{df}=d{f}_{\mathrm{now}}/d{f}_{\mathrm{before}}-1$$

  (5)

(3)

For $d_{kt}$ smaller than −0.5, $d_{df}$ is not taken from the matrix, but extrapolated as:

$$d_{df}=0.5d_{kt}{}^4-1.23d_{kt}{}^3+1.1d_{kt}{}^2-0.87d_{kt}$$

(6)

(4)

For $d_{kt}$ greater than 1, $d_{df}$ is extrapolated as:

$$d_{df}=-0.35-0.15d_{kt}$$

(7)

(5)

The diffuse fraction for part 2, $d{f}_2$, can now be calculated as:

$$d{f}_2=d_{df}d{f}_{\mathrm{before}}$$

(8)

3.3. Part Three. Geometric Calculation for Days with Clear Sky

3.3.1. Calculation of the Daily Course of $df$

In the case of clear sky, $df$ is mainly dependent on the air mass relative to its daily minimum. For that reason, in this part of the model a geometric approach has been chosen capable of reproducing the characteristic daily course of the diffuse fraction for clear sky days. The diffuse fraction of this part is calculated as:

$$d{f}_3=\frac{AM}{A{M}_{\min }}d{f}_{\min }$$

(9)

The air mass can be modelled as a function of the elevation of the sun. The minimal air mass $A{M}_{\min }$ is calculated for each day by using the maximum elevation angle ${\mathsf{\gamma }}_{\mathrm{S},\max }$:

$$A{M}_{\left(\min \right)}=\frac{1}{\sin {\left({\mathsf{\gamma }}_{\mathrm{S},\left(\max \right)}\right)}^{1.15}}$$

(10)

Figure 7 displays the measured (blue) and modelled (green dashed) course of the diffuse fraction over an exemplary day in Tateno, Japan (13 February 2006). While the clearness index $kt$ (black) is relatively stable around 1, the diffuse fraction is around 0.5 shortly after sunrise and before sunset and is falling down to a minimal diffuse fraction $d{f}_{\min }$ at noon, to 0.136 in this example.

However, the main challenge in modelling the diffuse fraction over the course of a clear sky day is to find a good approximation for the minimal diffuse fraction of the day, since this factor is subject to strong variations in every possible respect: from location to location, from season to season and even from day to day.

3.3.2. Daily Variation of $d{f}_{\min }$

In order to illustrate the daily variation of $d{f}_{\min }$, seven consecutive days in Tamanrasset, Algeria, are shown in Figure 8. Even for this non-cloudy site $d{f}_{\min }$ may vary significantly from day to day: The minimal value of the diffuse fraction (grey, bottom plot) of each day is varying between 0.328 on the first day (21 March 2006) and 0.062 on the third day (23 March 2006).

3.3.3. Seasonal Variation of $d{f}_{\min }$

In addition to daily variations, $d{f}_{\min }$ also shows seasonal variation on some locations. Figure 9 displays the minimal diffuse fractions of all clear days in Tamanrasset, Algeria, in 2006 (black plus symbols) over the course of a year. While in wintertime $d{f}_{\min }$ ranges mostly between 0.05 and 0.15, it almost never falls below 0.1 in summertime and features values between 0.15 and 0.5. When looking at the daily mean $kt$ values (grey crosses), no significant correlation can be observed which implies that other factors must have influence on the minimal daily diffuse fraction. The monthly means of aerosol optical depth (red) and the water vapor (blue dashed) taken from the NASA Terra/MODIS satellite [29,30] however feature a seasonal behavior similar to $d{f}_{\min }$.

3.3.4. Summary of Factors Influencing $d{f}_{\min }$

This leads to the conclusion that $d{f}_{\min }$ is dependent on a series of factors. A list of factors that proved influential on $d{f}_{\min }$ is given below:

(1)

The clearness index $kt$. The values of $kt$ are averaged in a range of 120 min around noon:

$$\overline{kt}=\frac{1}{t_{\mathrm{range}}}{{\displaystyle \sum }}_{i=t_{\mathrm{noon}}-\frac{t_{\mathrm{range}}}{2}}^{t_{\mathrm{noon}}+\frac{t_{\mathrm{range}}}{2}}k{t}_i$$

(11)

(2)

The variability of the clearness index, $k{t}_{\mathrm{Var}}$. For the same period of time, the changes of $kt$ are registered:

$$k{t}_{\mathrm{Var}}={{\displaystyle \sum }}_{i=t_{\mathrm{noon}}-\frac{t_{\mathrm{range}}}{2}}^{t_{\mathrm{noon}}+\frac{t_{\mathrm{range}}}{2}}\left|\frac{k{t}_i}{k{t}_{i-1}}-1\right|$$

(12)

(3)

The maximum elevation of the sun during the day, ${\mathsf{\gamma }}_{\mathrm{S},\max }$ and the minimum air mass during the day, $A{M}_{\min }$, compare to top of this section.

(4)

The aerosol optical depth (AOD) and the water vapour (wv) of the respective month. These values are taken from the NASA Terra/MODIS satellite [29,30] and averaged on a month per month basis between 2001 and 2015. Figure 10 gives an impression of the worldwide seasonal characteristics of AOD and wv.

(5)

The up and down time of $kt$ in the morning and in the evening. As a measure of the steepness of the $kt$ curve, the time span is determined between sunrise and when $kt$ first reaches the threshold of 1 in the morning. A second time span between the moment when $kt$ is at last above 1 in the evening and sunset is measured as well. The two values are averaged and are a good indicator for $d{f}_{\min }$in places with high day-to-day variation of $d{f}_{\min }$: The longer the up/down time, the higher $d{f}_{\min }$ will be.

3.3.5. Resulting Equations for $d{f}_{\min }$

The factors that influence $d{f}_{\min }$ mentioned in the above section are combined in a series of posynomials, depending on the location and availability of data. The coefficients and exponents of the following posynomials were fitted with the help of the Levenberg-Marquardt algorithm implementation of the software gnuplot 5.0 [31].

The datasets used for the posynomial fits are taken from the same locations as in Table 2, but for different years of measurement: asp 2009, bil 2005, bou 2010, brb 2011, cab 2011, cam 2002, clh 2014, cnr 2010, coc 2007, daa 2003, dar 2010, fua 2013, gob 2013, iza 2010, lau 2007, ler 2002, lin 2002, pal 2010, pay 2010, reg 2011, sap 2013, sbo 2011, sms 2006, sov 2002, tam 2003, tat 2003, tor 2005, xia 2006. For Case 1 only the subset bou 2010, iza 2010, sbo 2011, sov 2002, tam 2003 and xia 2006 was used.

Case 1: AOD and water vapor data available, location features strong seasonal changes of AOD (indicator for seasonal aerosol concentrations e.g., due to sandstorms in desert regions, see

Table 3. Values for a, b and c factors of the $d{f}_{\min }$ fit, used to model $d{f}_{\min }$ for given $kt$, $k{t}_{\mathrm{var}}$, AM, AOD, water vapour and up/down time. The RMS of residuals is 0.0528.

Factors	0	1	2	3	4	5
a	−4.29127	0.09656	−1.26822	0.05940	−0.30991	0.00043
b	0.19589	0.93797	0.03795	1.48181	0.08588	0.79801
c	6.01645	-	-	-	-	-

the values for factors a, b and c).

$$d{f}_{\min }=a_{0}k{t}^{b_0}+a_{1}k{t}_{\mathrm{var}}{}^{b_1}+a_{2}A{M}^{b_2}+a_{3}AO{D}^{b_4}+a_{4}w{v}^{b_4}+a_5{t}_{\mathrm{up}/\mathrm{down}}{}^{b_5}+c$$

(13)

Case 2: AOD and water vapor data available, no strong seasonal changes of AOD (see

Table 4. Values for a, b and c factors of the $d{f}_{\min }$ fit, used to model $d{f}_{\min }$ for given $kt$, $k{t}_{\mathrm{var}}$, AM, AOD, water vapour and up/down time. The RMS of residuals is 0.0427.

Factors	0	1	2	3	4	5
a	−2.49013	0.08345	0.00673	0.14107	−0.05853	0.00158
b	0.15065	0.72204	2.25298	0.75615	0.37413	0.67690
c	2.58895	-	-	-	-	-

the values for factors a, b and c).

$$d{f}_{\min }=a_{0}k{t}^{b_0}+a_{1}k{t}_{\mathrm{var}}{}^{b_1}+a_{2}A{M}^{b_2}+a_{3}AO{D}^{b_4}+a_{4}w{v}^{b_4}+a_5{t}_{\mathrm{up}/\mathrm{down}}{}^{b_5}+c$$

(14)

Case 3: AOD and water vapor data are not available (see

Table 5. Values for a, b and c factors of the $d{f}_{\min }$ fit, used to model $d{f}_{\min }$ for given $kt$, $k{t}_{\mathrm{var}}$, AM and up/down time. The RMS of residuals is 0.0480.

Factors	0	1	2	3	4	5
a	−0.75568	0.10744	0.02533	-	-	0.01203
b	0.16313	0.58318	1.26937	-	-	0.45174
c	0.71854	-	-	-	-	-

the values for factors a, b and c):

$$d{f}_{\min }=a_{0}k{t}^{b_0}+a_{1}k{t}_{\mathrm{var}}{}^{b_1}+a_{2}A{M}^{b_2}+a_5{t}_{\mathrm{up}/\mathrm{down}}{}^{b_5}+c$$

(15)

Case 4: AOD, water vapor and up/down time data are not available (see

Table 6. Values for a, b and c factors of the $d{f}_{\min }$ fit, used to model $d{f}_{\min }$ for given $kt$, $k{t}_{\mathrm{var}}$, AM and up/down time. The RMS of residuals is 0.0542.

Factors	0	1	2	3	4	5
a	−2.28942	0.23589	0.02445	-	-	-
b	0.27308	0.19371	1.26262	-	-	-
c	2.23274	-	-	-	-	-

the values for factors a, b and c).

$$d{f}_{\min }=a_{0}k{t}^{b_0}+a_{1}k{t}_{\mathrm{var}}{}^{b_1}+a_{2}A{M}^{b_2}+\mathrm{c}$$

(16)

3.4. Combination of the Three Parts

The three parts of the algorithm generate the values $d{f}_1$, $d{f}_2$ and $d{f}_3$. Depending on the current weather situation, expressed by characteristics and statistical features of $kt$, they are combined to one single, resulting $df$:

$$df=w_{1}d{f}_1+w_{2}d{f}_2+w_{3}d{f}_3$$

(17)

The mean absolute deviation of $kt$ at a given time of day $t_{\mathrm{x}}$ that is used as condition above is calculated as follows:

$$ma{d}_{kt}=\frac{1}{t_{\mathrm{range}}}{{\displaystyle \sum }}_{i=t_{\mathrm{x}}-t_{\mathrm{range}}}^{t_{\mathrm{x}}}\left|\frac{k{t}_i}{k{t}_{i-1}}-1\right|$$

(18)

with $t_{\mathrm{range}}=30\min$. For illustration of the weighing conditions mentioned in

Table 7. Weighing factors for the combination of $d{f}_1$, $d{f}_2$ and $d{f}_3$ to one single $df$, depending on $kt$ characteristics.

Name	Condition	w1	w2	w3
Clear Sky	$ma{d}_{kt}<ma{d}_{kt,\mathrm{lower}}$ $k{t}_{\mathrm{clear},\mathrm{lower}}<kt<k{t}_{c\mathrm{lear},\mathrm{upper}}$	0	0.2	0.8
Transition	$ma{d}_{kt,\mathrm{lower}}<ma{d}_{kt}<ma{d}_{kt,\mathrm{upper}}$ $k{t}_{\mathrm{clear},\mathrm{lower}}<kt<k{t}_{c\mathrm{lear},\mathrm{upper}}$	0.2	0.2	0.6
Standard	Else	0.2	0.8	0

With $k{t}_{\mathrm{clear},\mathrm{lower}}=0.95$, $k{t}_{\mathrm{clear},\mathrm{upper}}=1.2$, $ma{d}_{kt,\mathrm{lower}}=0.005$ and $ma{d}_{kt,\mathrm{upper}}=0.05$.

, Figure 11 displays all-sky camera images from the Institute of Meteorology and Climatology of the Leibniz University Hannover [32,33]. Picture A shows a moment where no clouds are visible. It is classified as “Clear Sky” since $kt=1.03$ and $ma{d}_{kt}=0.0025$. Picture B shows a moment where $kt=0.147$ and $ma{d}_{kt}=0.107$, hence being classified as “Standard”. In picture C some light clouds are visible around the sun. This moment is classified as “Transition” as $kt=1.01$ and $ma{d}_{kt}=0.028$. The “Transition” condition can be interpreted as clear sky with only few light clouds.

The generated matrices and other model data can be obtained from the authors upon request.

4. Results

In this section, the results of the validation of the new algorithm are presented. As mentioned in Section 2.1, the validation is conducted for 28 locations with one year of one-minute values each, basing the validation on more than seven million data points worldwide.

The overall results are then compared to the results of three existing models for the diffuse fraction: the model of Orgill and Hollands [18] (OH), a one-parameter model, the reduced version of the two-parameter model of Reindl et al. [3] (RR), and the model by Perez and Ineichen [23] (PZ), all introduced in Section 2.3. The first two models were also identified as two of the three best performing models among the eight investigated approaches by Dervishi [27]. The model by Perez and Ineichen [23] is still popular in the community and widely made use of. The model by Skartveit [12] was not used for the model comparison since no indications were found that show a significant advantage of this model over Orgill and Hollands [18], Reindl et al. [3] or Perez and Ineichen [23].

Figure 12 displays two weeks in Alice Springs, Australia, 2005. The measured global horizontal irradiance is plotted on top (green); the resulting clearness index $kt$ is plotted for reference underneath (black). In the three following plots, the measured diffuse fraction (black) is displayed, together with the diffuse fraction that was modelled with the new approach (blue), with the model from Orgill and Hollands [18] (grey) , the model from Reindl et al. [3] (orange) and the model from Perez and Ineichen [23].

While the new model is able to reproduce the diffuse fraction in good accordance to the measurement values most of time, the inherent problem of models with static one- or two-parameter relationships between the clearness index and the diffuse fraction becomes apparent. Especially on clear sky days the existing models fail to reproduce the characteristic behavior of the diffuse fraction.

In order to evaluate the performance of the new model in statistical terms, the root mean squared errors (RMSEs) are calculated for the new model as well as for the three reference models. Figure 13 shows the RMSE for the four models over all test data sets. The RMSE produced by the new model is smaller than those produced by the models of Orgill and Hollands [18], Reindl et al. [3] and Perez and Ineichen [23], in parts significantly, except for one case in Izaña, Spain, 2011 (iza 2011). The overall RMSE, averaged over all test data sets, can be reduced from 0.138 (OH), 0.134 (RR) and 0.139 (PZ) to 0.116 for the new model, which equals an amelioration of 16%–20%.

A further validation is conducted by comparing the annual diffuse irradiation values that are estimated by the models with the measured value. Figure 14 lists the relative deviations of the modelled from the measured annual diffuse irradiation. In most of the cases, the deviation resulting from the new model is significantly smaller than the deviation resulting from the models of OH, RR or PZ. There are few cases where the model leads to higher deviations than the existing ones, e.g., for Billings, USA (bil 2003), Solar Village, Saudi Arabia (sov 2001) or Tamanrasset, Algeria (tam 2006). Extreme deviations of more than 20%, however, as apparent in some of the test cases for the two existing models, do not occur when using the new model. The average of the absolute (i.e. unsigned) relative deviations for all test cases can be reduced by nearly 50% from 11.9% for OH, 12.7% for RR and 10.9% for PZ to only 6.4% for the new model.

The histogram of the mean absolute deviations of the annual diffuse irradiation displayed in Figure 15 illustrates the frequency of the deviations each model produces. While the model of Reindl et al. [3] (RR) has its peak in the class of 0 to 5%, it still has several outliers of 35%–55%. The model of Orgill and Hollands [18] (OH) features only one extreme outlier at 35%–40%, but has most of its results lying in the class of 10 to 15%. The model by Perez and Ineichen [23] shows no outliers of more than 25% but has its results evenly distributed between 0 and 15%. The new model does not produce any outliers and has its peak in the class of 0%–5%, covering 50% of the test cases alone.

This histogram can be converted into confidence probabilities, depicted in Figure 16. The new model results in mean absolute deviations of less than 5% in 50% of the cases. In more than 80% of the cases, the deviation is less than 10%. Compared to the other three models this is a significant improvement, where the probability of producing less than 5% deviation is only 25% (OH), 36% (RR) and 25% (PZ) and the probability of producing less than 10% deviation is 36% (OH), 54% (RR) and 50% (PZ).

5. Conclusions

The newly developed model for the diffuse fraction of solar irradiance on PV systems provides significantly better agreement with measurements than the other models published so far. This is achieved by the following features: the first part utilizes the dependency of the diffuse fraction $df$ on the clearness index $kt$, in analogy to existing one-parameter models. In the new model, the correlation is expressed as probability matrices rather than single functions, leading to realistic, more natural diffuse fraction characteristics. Also taking advantage of probability matrices, the second part uses the relation of the relative changes of $df$ over the relative changes of $kt$. The third part takes into account the diffuse fraction characteristics of days with clear sky only using a geometrical approach. The crucial factor for the third part is the minimum daily diffuse fraction for which a posynomial model has been introduced.

The presented new model was analyzed and compared to two other models for 28 locations worldwide with one year of one-minute measurement data each. It was shown that the new model has a high quality of modeling the diffuse irradiance. The mean RMSE over all test cases was reduced by 16%–20%, whereas the mean absolute deviation of the annual diffuse irradiation was found to be nearly 50% smaller compared to the reference models. In more than 80% of the test cases, the deviation of the annual diffuse irradiation is smaller than 10%, with an overall maximum deviation of 20%.

With the new model, the diffuse irradiance can be calculated with much lower uncertainty, hence significantly reducing the uncertainty of PV energy yield simulations. Possible future work for the improvement of the model will include further investigations on the minimal daily diffuse fraction that has a very decisive influence on the model quality for clear sky days. Such days may be identified by cloud cameras that will allow for a much better estimation of the cloud fraction compared to satellite images [33].