## 2 Materials and Methods  
This Section gives a full account of the new (mathematical) method that accurately determines the upper and lower limits of \(\rm{k_{d}}\) at a site. To apply the methodology, a typical meteorological year (TMY) must preferably be used instead of a long period of measurements or model simulations. This peculiarity is based on the fact that a TMY refers to all representative situations of the included parameters that occur within a year, while a long-term series of measurements or model simulations contains extreme values of the included parameters; the second, therefore, option may result in a broader dispersion of the calculated \(\rm{k_{d}}\) values. On the other hand, TMYs are nowadays used more often in many climatological studies, e.g., [11; 12; 13]. The selected solar radiation parameter to vary with \(\rm{k_{d}}\) is the direct horizontal solar radiation (\(\rm{H_{b}}\)) because of much less dispersion of the \(\rm{H_{b}}\)-\(\rm{k_{d}}\) paired values in comparison with the cases of the functions \(\rm{H_{g}}\) (global horizontal solar radiation) vs. \(\rm{k_{d}}\) or \(\rm{H_{d}}\) (diffuse horizontal solar radiation) vs. \(\rm{k_{d}}\). Even in the case of the direct horizontal solar radiation, the scatter plot of \(\rm{H_{b}}\) vs. \(\rm{k_{d}}\) is found to be greater than that of the relationship \(\rm{H_{bn}}\) (direct-normal solar radiation) vs. \(\rm{k_{d}}\). Figure 1 shows this difference; here hourly mean \(\rm{H_{b}}\), \(\rm{H_{bn}}\) and \(\rm{k_{d}}\) values have been plotted for the typical month of January of the TMY-BOU (Boulder, CO, USA) derived from the PV-GIS platform [14; 15] (see below for description). Therefore, the \(\rm{H_{bn}}\) solar component is considered in the rest of the analysis. The mathematical determination of the upper and lower limits of \(\rm{k_{d}}\) in a diagram as that of Figure 1b is fully described here. The first step is to find the best-fit curve to the \(\rm{H_{bn}}\)-\(\rm{k_{d}}\) data points. A second-order polynomial of the form \(\rm{H_{bn}}=a\;{k_{d}}^{2}+b\;{k_{d}}+c\) was used, where \(\rm{a}\), \(\rm{b}\), \(\rm{c}\) are the polynomial coefficients. This type of best-fit curve was selected as the results obtained showed that it is sufficient. Figure 2 shows the best-fit black curve to the scatter plot of Figure 1b, which has the expression \(\rm{H_{bn}}=902.30\;{k_{d}}^{2}-2087.44\)\(\rm{k_{d}}+1188.58\), \(\rm{R^{2}}=0.996\). As WMO recommends that the sunshine duration is measured when \(\rm{H_{bn}}>120\;{Wm^{-2}}\), the green line represents this threshold. This straight line crosses the best-fit curve at a certain point (\(\rm{k_{d}}\),\(\rm{H_{bnWMO}}=120\)). Since the y-axis value is known (120 \(\rm{Wm^{-2}}\)), solving the above second-order equation for \(\rm{k_{d}}\), two roots are obtained and  
Figure 1: Variation of (**a**) \(\rm{H_{b}}\) and (**b**) \(\rm{H_{bn}}\) as function of \(\rm{k_{d}}\) for January of the PV-GIS TMY-BOU. All data points are hourly values. Dispersion is greater in (**a**). the lower-than-1 is considered because k\({}_{\rm d}\) cannot go higher than 1 by definition; in this case, the accepted solution is k\({}_{\rm du}\) = 0.76. The next step is to find the equation of the straight line that is tangent to the best-fit curve at the point C (0.76,120); this point is the crossing point between the best-fit curve and the WMO-defined threshold for SSD measurements. This line expresses the first derivative of the best-fit curve at the point C; in other words, it represents the slope of the best-fit line at C. This is done by taking the first derivative of the H\({}_{\rm bn}\) expression (H\({}^{\prime}\)\({}_{\rm bn}\)) in respect to k\({}_{\rm d}\) and replacing k\({}_{\rm d}\) with 0.76. It is found that H\({}^{\prime}\)\({}_{\rm bn}\) = \(-\)715.94 Wm\({}^{-2}\). Then, H\({}_{\rm bn}\) - 120 = \(-\)715.94 (k\({}_{\rm d}\) - 0.76) from which the expression for the tangent line is H\({}_{\rm bn}\) = \(-\)715.94 k\({}_{\rm d}\) + 664.12.