Artefact in Ineichen clearsky model for sun close to horizon

[kdebrab]

Consider following example:

import pandas as pd
import pvlib
times = pd.date_range('2018-02-16', '2018-02-17', freq='min')
location = pvlib.location.Location(latitude=60., longitude=0.)
for model in ['ineichen', 'haurwitz', 'simplified_solis']:
    location.get_clearsky(times, model=model)['ghi'].plot(label=model, legend=True)

which results in following figure:

The (default) 'ineichen' clearsky model seems erratic at sunrise and sunset. The effect is especially visible for high latitudes during Winter when the sun slowly approaches the horizon.

Upon closer look in the code, I found the following comment:

pvlib-python/pvlib/clearsky.py

Lines 98 to 103 in 04bd753

	# This equation is found in Solar Energy 73, pg 311. Full ref: Perez
	# et. al., Vol. 73, pp. 307-317 (2002). It is slightly different
	# than the equation given in Solar Energy 73, pg 156. We used the
	# equation from pg 311 because of the existence of known typos in
	# the pg 156 publication (notably the fh2-(TL-1) should be fh2 *
	# (TL-1)).

which refers to the formula:

pvlib-python/pvlib/clearsky.py

Lines 122 to 126 in 04bd753

	ghi = (np.exp(-cg2*airmass_absolute*(fh1 + fh2*(tl - 1))) *
	np.exp(0.01*airmass_absolute**1.8))
	# use fmax to map airmass nans to 0s. multiply and divide by tl to
	# reinsert tl nans
	ghi = cg1 * dni_extra * cos_zenith * tl / tl * np.fmax(ghi, 0)

However, apart from the typo mentioned, there is another important difference between both formula of Solar Energy 73, and that's the np.exp(0.01*airmass_absolute**1.8) factor, which is only present in the formula of pg 311, but not in the formula of pg 156!

It seems like the artefact is exactly a result of that np.exp(0.01*airmass_absolute**1.8) factor. This factor slowly increases from 1.01 at solar zenith to 1.95 at a zenith angle of 85°, but then shoots to 1000 when the sun reaches the horizon.

So, maybe there was an error also in the formula of Solar Energy 73, pg 311, and that factor wasn't supposed to be there?

Nevertheless, in reference [3] (Reno et al., 2012), this model was found to excel for US sites with the same formula as used by pvlib. The US sites were all below the 49th parallel. I wonder whether this model is still recommendable for higher latitudes?

[kdebrab]

There is actually another version of the paper of Ineichen, which doesn't contain the typo nor the factor causing the artefact: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.535.4248&rep=rep1&type=pdf.

[adriesse]

Hi Karel! I saw this effect too but simply didn't have the time to investigate and report--so thank you for this! One hypothesis I had was that this might have worked ok for hourly averages.

[kdebrab]

Hi Anton, The effect indeed mostly disappears in hourly averages.

Unless you go even further north. Then the whole curve gets really biased. See e.g. the result for 70° latitude:

[wholmgren]

Thanks for the detailed report. This artifact and the discrepancies among the papers and code is part of what motivated me to add the simplified solis model to the library.

[cwhanse]

And is an example of the original motivation for PVLib. We intended the open source libraries to expose and hopefully resolve this kind of ambiguity in models. @kdebrab do you see a resolution among the various papers? When I can find time I'll pull the Solar Energy papers and see what we can determine on our end. One of us should write Pierre and see if he can enlighten us.

[mikofski]

Hi @kdebrab thanks for catching this. How does Ineichen compare to haurwitz and solis when the mystery factor np.exp(0.01*airmass_absolute**1.8) is removed? And how do they compare to Bird?

[adriesse]

Pierre will retire in May, so don't wait too long to ask!

[mikofski]

Gasp!

[adriesse]

I checked with Pierre Ineichen, and he explained that Richard Perez added the exponential term in the later paper based on further analysis of measurements at 10 continental US sites. I guess no one checked the effect on locations outside this geographic area.

I would suggest that function be modified to accept a flag to indicate whether the Perez enhancement is to be used or not, with a default value of False.

[wholmgren]

Thanks for checking. I support your proposed solution. I'm not going to get to it anytime soon, though, and would welcome a pull request.

[cwhanse]

I concur with @adriesse proposal. I would add a comment to the function signature warning that the Perez enhancement may produce biased results at high latitude, and link the False/True condition for the flag to the appropriate reference.

[martinherrerias]

One way to remove the artifacts while introducing less bias is to clip the airmass_absolute values to the point where the exponential correction becomes problematic. Something like:

airmass_critical = (cg2*(fh1+fh2*(tl-1))/0.018).^1.25
airmass_clipped = np.fmin(airmass_absolute, airmass_critical)
ghi = (np.exp(-cg2*airmass_clipped*(fh1 + fh2*(tl - 1)) + 0.01*airmass_clipped**1.8))
...
ghi = cg1 * dni_extra * cos_zenith * tl / tl * np.fmax(ghi, 0)

airmass_critical is a minimum for the combined exponential terms, which shoot up afterwards. Loosely:
airmass_critical = argmin(f(x)), for f = exp(-cg2*x*(fh1 + fh2*(tl - 1)) + 0.01*x^1.8))

This threshold starts at AM = 2.6 (for turbidity = 1 and altitude = 0), equivalent to ~22.5° solar elevation, and goes up to O(10) for higher turbidities and altitudes. For most cases, it corresponds to solar elevation angles below the 4° cutting threshold of the CIE 1994 quality control guidelines that Ineichen & Perez seem to have used, i.e. they probably didn't see beyond it.

Clipping at the minimum makes the adjustment smooth (C1), and leaves the results of the model untouched for most operating conditions.

[cwhanse]

@martinherrerias your proposal makes a lot of sense to me, now that I better understand the the role that correction plays. Rather than continuing in this old issue, would you like to open a new issue? If you are also inclined a pull request is welcome.

[adriesse]

I would suggest controlling this new feature by the same flag or a new one. It should be possible to run the model as published.

This reminds me of some discussions about using models outside of the range of inputs for which they were developed and/or validated.

[mikofski]

Are there any references for this new feature? @martinherrerias any plans to present or publish this idea? I think that would be useful. Thanks!