High-quality direct and diffuse radiation data are essential for the development and utilization of solar energy resources in a region and serve as important parameters to characterize its climatic features (WMO, 2014). However, in reality, direct radiation measurements require sun-tracking devices, and diffuse radiation measurements necessitate shading devices, making the cost expensive and maintenance challenging. As a result, many observation sites only conduct total solar radiation observations, making it difficult to obtain long-term series of direct and diffuse radiation data. The diffuse fraction represents the proportion of diffuse radiation relative to the total solar radiation in the sky and is closely related to both direct and diffuse radiation. The Ångström-Prescott formula indicates a linear relationship between the clear-sky index and corresponding relative sunshine hours, with a negative correlation between cloud cover and sunshine duration (Ångström, 1924; Prescott, 1940). Studying the diffuse radiation model based on known or easily accessible parameters such as the clear-sky index allows for the analysis of the diffuse fraction pattern. This research holds significant importance for the planning, design, and operation of solar power plants, as well as for studying climate change. Additionally, it plays a crucial role in enhancing the yields of agricultural and forestry crops like rice, wheat, and apples, improving building energy efficiency, and promoting the development of regional tourism and leisure industries, thereby presenting essential economic and social value.

Currently, research on the regularity of diffuse fraction mainly focuses on constructing diffuse fraction models using piecewise linear and polynomial fitting methods. Shang et al. investigated the direct and diffuse fraction models of horizontal solar radiation in low-latitude areas and established a piecewise linear regression model between atmospheric diffuse transmittance ratio and clear-sky index (Shang et al., 2017). Feng et al. compared five representative direct and diffuse fraction models in Beijing and identified the main influencing factors and the most suitable model for Beijing (Feng et al., 2015). Ma et al. improved the Erbs direct and diffuse solar radiation models by introducing the sunshine percentage and establishing a polynomial model, which enhanced the model accuracy to some extent (Ma et al., 2016). Fu et al. found that the solar radiation duration was highly correlated with the diffuse radiation ratio in regions with severe atmospheric pollution, with correlation coefficients exceeding −0.8, an increase in the concentration of aerosol particles in the air has a significant impact on sunshine duration and SDRF over China. (SDRF—defined as the diffuse solar radiation divided by the total solar radiation in conditions of sunny weather) (Fu and Dan, 2018). Zhang and Zhao conducted a study on the transmission and extinction characteristics of mid-infrared radiation in the atmosphere, analyzing the influence of different atmospheric profiles, molecular attenuation, and aerosols on radiation transmission (Zhang et al., 2017; Zhao et al., 2018). Yvonne et al. investigated the use of total radiation to invert sunshine duration through a piecewise linear model, which could be used when there was no sunshine observation data available (Yvonne and Wouter, 2007; Chong et al., 2017). Ahwide et al. studied the linear relationship between the daily average diffuse fraction and the clear-sky index, predicting the daily direct radiation and sunshine duration in Sabhā, Ghadāmis, and Tripoli in Libya (Ahwide et al., 2013). Basharat and Abid proposed a generic model for estimating diffuse fraction in India using clear-sky index and sunshine duration, based on linear and polynomial fitting methods (Basharat and Abid, 2017). Li et al. established a linear relationship between the inclined total solar radiation and sunshine duration in Hong Kong (Li et al., 2008). Manuel et al. proposed a piecewise linear effect of atmospheric diffuse, absorption, and optical quality on the relationship between total solar radiation and diffuse radiation (Manuel et al., 1991). Muneer and Munawwar developed a second-order polynomial empirical formula for the components of daily diffuse radiation using statistical methods and data from India, Japan, Spain, and the UK, including sunshine duration, cloud distribution coefficient, and clear-sky index (Muneer and Munawwar, 2007). Fu et al. analyzed the trends in sunshine duration and diffuse fraction under clear sky conditions in seven city clusters in China between 1960 and 2005 and pointed out that anthropogenic activities and air pollution, such as aerosol optical thickness and tropospheric NO2 concentration, had a highly negative correlation with sunshine duration and diffuse fraction (Fu et al., 2015). Khorasanizadeh et al. estimated the diffuse radiation in Kerman, Iran, using polynomial and exponential fitting functions, finding that a first-order linear formula performed best for daily diffuse fraction and monthly average diffuse coefficient simulations (Khorasanizadeh et al., 2016). Fan et al. estimated horizontal diffuse radiation in China’s five climatic regions using polynomial and exponential fitting functions and concluded that different function models were suitable for different climatic regions (Fan et al., 2019). Alexandre et al. proposed a fourth-order polynomial model for calculating hourly, daily, and monthly average diffuse radiation using the clear-sky index (Alexandre et al., 2014). Finally, Kambezidis et al. (2021) have determined the upper and lower limits of the diffuse fraction for classitying skies into overcast, intermediate and clear with a wordwide validation.

The Boland Ridley Laurent (BRL) model, proposed by Boland and Ridley, is a commonly used S-shaped function, namely the Sigmoid function, with values ranging between 0 and 1. Its basic form is f(x) = 1/(1 + e^ (-x)). As the independent variable approaches positive or negative infinity, the function output tends to a smooth state of 0 and 1 respectively, which is similar to the distribution of diffuse fractions, making it suitable for constructing diffuse fraction models. In 2008, Boland and Ridley utilized hourly and 15-min data from the Geelong region in Australia to propose the BRL model for diffuse radiation values (Boland et al., 2008). In 2019, they further employed this model to estimate diffuse fractions in Chile, Costa Rica, and Australia, evaluating its correlation with climate factors such as temperature, precipitation, and latitude (Rojas et al., 2019). Based on previous research (Ridley et al., 2010), compared to the traditional piecewise linear fitting and fourth-order polynomial fitting, the BRL model performs better in calculating low diffuse fractions under high clear-sky index conditions, as the former two methods exhibit underfitting issues when the diffuse fractions approach zero. Considering that the residuals between the calculated and the observed values from different models do not follow a normal distribution, the evaluation metrics selected are the median absolute percentage error for measuring the error magnitude and the median deviation for comparing the performance of different models. Research in the southern hemisphere shows that, in terms of median absolute percentage error, the BRL model is 3.07% better than the piecewise linear fitting and 15.48% better than the fourth-order polynomial fitting. Regarding the median deviation, the BRL model is 0.01 smaller than the piecewise linear fitting and 0.08 smaller than the fourth-order polynomial fitting. It is evident that the BRL diffuse fraction model has advantages in fitting performance, providing a new research direction for the study of diffuse solar radiation.

In this paper, based on the BRL model and atmospheric radiation transmission theory, the Taylor series expansion is utilized to represent the diffuse fraction relative to the clear-sky index and atmospheric optical quality. Independent model training datasets and validation datasets are established at different time scales, including minute, 10-min, and hourly data. Nonlinear fitting is performed on the training datasets to construct the binary BRL model for diffuse fractions relative to the clear-sky index and atmospheric optical quality. The validation datasets are used to calculate and analyze the correlation coefficient, average deviation, root mean square error, and t-statistic between the observed and calculated values of the diffuse fractions, thereby evaluating the effectiveness and applicability of the binary BRL model.

This paper is based on the BRL model, and the ratio of the ground-level direct horizontal radiation component to the total horizontal solar radiation at the ground surface is defined as the direct fraction, as shown in Equation 5, according to the principle that the total horizontal solar radiation at the ground surface is equal to the sum of the normal direct radiation component and the scattered radiation, as given in Equations 3, 4. I g = I d + I b (3) I b = I n ⁡ sin ⁡ α (4) K b = I b I g (5)

In the equation, I b represents the ground-level direct horizontal radiation component, I n represents the normal direct radiation reaching the ground, and K b represents the direct fraction. Therefore, the relationship that the sum of the diffuse fraction and the direct fraction equals 1 can be easily derived, as shown in Equation 6. K d + K b = 1 (6)

According to the theory of atmospheric radiation transfer, solar radiation reaching the ground undergoes a complex process of attenuation as it passes through the atmosphere, particularly due to Rayleigh scattering and absorption by various components such as gases, water vapor, and aerosols present in the atmosphere. Assuming that the extinction processes, including diffuse and absorption, are mutually independent, the Lambert-Beer law is extended to a broad spectral range, and the normal direct radiation reaching the ground is given as (Ruiz–Arias et al., 2010). I n = ∏ i I 0 e − m i δ i (7)

In the equation, m i represents the optical mass corresponding to the i -th extinction process, and δ i represents the optical thickness corresponding to the i -th extinction process. By substituting Equations 1, 4, 7 into Equation 5, the direct fraction, denoted as K b , can be transformed into: K b = 1 K t ⋅ ∏ i e − m i δ i = 1 K t ⋅ e ∑ i − m i δ i (8)

According to Equations 6, 8 and assuming that the extinction processes are categorized into two parts, namely diffuse and absorption, the diffuse fraction can be transformed into: K d = 1 − 1 K t ⋅ e − m r δ r + m a δ a (9)

In the equation, m r represents the optical mass corresponding to the extinction process caused by Rayleigh scattering, δ r represents the optical thickness corresponding to the extinction process caused by Rayleigh scattering, m a represents the optical mass corresponding to the extinction process caused by atmospheric absorption, and δ a represents the optical thickness corresponding to the extinction process caused by atmospheric absorption.

Due to the clearness index K t ∈ 0 , 1 , the decomposition according to the Taylor series is as follows: 1 K t = 1 1 − 1 − K t = ∑ j = 0 ∞ 1 − K t j (10)

Equation 10 only converges when 1 − K t < 1 is considered, and by taking j = 0 and normalizing the absorption extinction effect with the diffuse extinction effect, represented by τ , Equation 9 can be expressed as follows: K d = 1 − e − m r δ r 1 + τ (11)

Based on the logarithmic identity and Taylor series expansion, Equation 11 can be transformed into: K d = 1 − e − e ln m r δ r 1 + τ = 1 − e − e ln m r δ r + ln 1 + τ = 1 − e − e C + τ − 1 2 τ 2 + 1 3 τ 3 − 1 4 τ 4 + ⋅ ⋅ ⋅ (12)

Equation 12 represents the derived diffuse fraction BRL general model, where C is a constant, and the polynomial in the exponential part approximates the extinction processes related to atmospheric diffuse and absorption. The specific model should incorporate variables that contribute significantly to atmospheric extinction in the research.

The present study proposes two main contributions. Firstly, the model incorporates the clearness index K t , which represents the ratio of ground total radiation to extraterrestrial solar radiation. Secondly, considering the significant impact of the relative path length of direct radiation in the atmospheric extinction process, the model introduces the atmospheric optical mass as the second variable alongside the clearness index K t . Therefore, the proposed model is a binary BRL model for the diffuse fraction, as shown in Equation 13. K d K t , m = a 1 − a 2 ⋅ e − e a 3 + a 4 ⋅ K t + a 5 ⋅ m (13)

In the equation, a1∼a5 are the models parameters, and m is the atmospheric optical mass.

The atmospheric optical mass here refers to the relative atmospheric optical mass, which represents the ratio of the distance that solar radiation travels through the atmosphere to reach the sea level to the distance it would travel through the atmosphere to reach the sea level when the sun is at the zenith. The atmospheric optical mass is closely related to the solar zenith angle and the elevation of the site, and its precise calculation requires corrections for air pressure and atmospheric refraction, which are determined through Equations 14–16 (Rigollier et al., 2000; Chen et al., 2017; Kasten et al., 1989). m = p / p 0 sin ⁡ α s + 0.50572 α s + 6.07995 − 1.6364 (14)

In the equation, α s represents the solar zenith angle after atmospheric refraction correction, which can be calculated using Equation 16. p is the atmospheric pressure at the observation site, and p 0 is the average atmospheric pressure at sea level. The latter is involved in the air pressure correction, which is particularly crucial for sites located at higher elevations such as mountainous areas, and it can be calculated using Equation 17. α s = α + 0.061359 180 π × 0.1594 + 1.123 180 π α + 0.065656 180 π 2 α 2 1 + 28.9344 180 π α + 277.3971 180 π 2 α 2 (15) p / p 0 = e − h / h z (16)

In the equation, α represents the solar zenith angle, h represents the elevation of the observation site, and h z represents the altitude of the nearby Rayleigh atmosphere, taken as 8,434.5 m. A more accurate estimation of the solar elevation due to atmosphere refraction is given by Kambezidis (2021).

Next, the model parameters are solved. Using Matlab programming, the training dataset is imported, and through nonlinear fitting, the binary BRL model function is derived. Table 1 lists the solved parameters of the binary BRL model, with one set for each of the hourly, 10-min, and minute-scale radiation data, making a total of three sets for each model. The binary BRL model surfaces are shown in Figure 1, where (a)-(c) represent the surfaces fitted to the hourly, 10-min, and minute-scale data, respectively.

The distribution relationship between the diffuse fraction and the clearness index ( K t ) and atmospheric optical mass (m) was investigated across three different time scales. A comparison of the parameters obtained from the 10-min and minute-scale data in Table 1 reveals a close similarity, whereas they significantly differ from the parameters derived from hourly data. These discrepancies in parameters correspond to different model forms. As observed from Figure 1, in the BRL model, the clearness index ( K t ) ranges between (0, 1), the atmospheric optical mass (m) ranges between (0, 6), and the function values range from (0, 1). The model shows a monotonic decreasing trend concerning the clearness index ( K t ), while the atmospheric optical mass (m) does not significantly affect the monotonicity of the diffuse fraction. The clearness index (x) is divided into three intervals: 0 < K t < 0.3 represents overcast weather, 0.3 ≤ K t < 0.7 represents cloudy weather, and 0.7 ≤ K t < 1 represents sunny weather. Under overcast and sunny weather conditions, the binary BRL model exhibits a relatively gradual change in the diffuse fraction with varying clearness index ( K t ). However, under cloudy weather, the diffuse fraction displays nearly linearly decreasing behavior concerning the clearness index ( K t ), albeit with different slopes at the transition points for different model parameters. It is crucial to note that the binary BRL surface is not a flat surface but exhibits certain distortions associated with the atmospheric optical mass (m). Additionally, the diffuse fraction is primarily concentrated within the range of atmospheric optical mass ( m < 2 ).

To assess the performance of the binary BRL model in calculating the diffuse fraction, validation computations were conducted using a validation dataset. Since the proposed binary BRL model in this study consists of three sets of parameters obtained from hourly, 10-min, and minute-scale radiation data, the validation process required investigating the impact of different time scales on model construction. Additionally, the validation dataset considered the model’s applicability across various seasons and regions. To achieve this, a separate set of 10-min and hourly radiation data measured in May 2022 at the Lijiang station and hourly radiation data measured in July 2022 and 10-min radiation data measured in December 2022 at the Kazuo station were selected for comparison with the training data. As minute-scale data were unavailable, the 10-min data was used for validation due to its similarity to minute-scale data in model construction. The validation process employed a series of statistical analyses, including correlation coefficient r , mean bias, root mean square error, t-statistic, skewness, and kurtosis (Xing et al., 2022; Chen et al., 2017). These analyses were used to statistically compare the measured diffuse fraction with the diffuse fraction calculated by the binary BRL model, thereby evaluating the model’s computational performance.

The correlation coefficient r was utilized to analyze the linear relationship between the measured diffuse fraction and the diffuse fraction calculated by the binary BRL model, as calculated using Equation 17: r = ∑ i = 1 N K dc i − K dc ¯ K dm i − K dm ¯ ∑ i = 1 N K dc i − K dc ¯ 2 ∑ i = 1 N K dm i − K dm ¯ 2 (17)

In equation, N represents the data quantity, K dc i is the i -th diffuse fraction calculated by the S0 or S1 model, K dc ¯ is the average value of N K dc , K dm i is the i -th measured diffuse fraction, and K dm ¯ is the average value of N K dm . As r approaches 1, it indicates that the model’s computational performance is better.

The Mean Bias Deviation (MBD) is used to analyze the systematic error performance of the binary BRL model in calculating the diffuse fraction. The closer the MBD value is to 0, the smaller the systematic error in the model’s calculation of the diffuse fraction. When MBD > 0, it indicates an overestimation phenomenon in the calculated diffuse fraction; when MBD < 0, it indicates an underestimation phenomenon. The MBD is calculated using Equation 18: M B D = ∑ i = 1 N K dc i − K dm i n (18)

The Root Mean Squared Deviation (RMSD) is used to analyze the random error performance of the binary BRL model in calculating the diffuse fraction. The RMSD can be computed using Equation 19, and a smaller value indicates a lower random error in the model’s calculation of the diffuse fraction, indicating higher precision. R M S D = ∑ i = 1 N K dc i − K dm i 2 N (19)

Although MBD and RMSD can evaluate the performance of the binary BRL model in calculating the diffuse fraction from the aspects of accuracy and precision, they have inherent limitations. For example, in cases of overall overfitting or underfitting, it is possible to observe very small RMSD values with significantly large MBD values. Additionally, when fitting a large number of discrete values that can cancel out positively and negatively, small MBD and large RMSD can coexist. Therefore, this study introduces the t-statistic, as shown in Equation 20, to comprehensively evaluate the model’s performance in calculating the diffuse fraction, considering both MBD and RMSD. A smaller t-statistic indicates better performance. t = N − 1 M B D 2 R M S D 2 − M B D 2 (20)

Furthermore, the skewness and kurtosis of the residuals between the measured diffuse fraction and the diffuse fraction calculated by the binary BRL model were calculated in this study. Skewness and kurtosis provide information about the symmetry and peakness of the residual distribution relative to a standard normal distribution. Ideally, the skewness and kurtosis values of a standard normal distribution are both 0. Therefore, smaller values of skewness and kurtosis indicate better performance of the model in calculating the diffuse fraction.

Table 2 presents the results of the statistical analysis between the calculated diffuse fraction using the binary BRL model and the measured diffuse fraction. All correlation coefficients are above 0.8, MBD values are within ±0.2, RMSD values are within 0.25, and the minimum t-statistic is as low as 0.1074.

When validated using the hourly and 10-min data from the Lijiang station in May 2022, the binary BRL model constructed using the 10-min training data exhibits the best overall performance, with correlation coefficients of 0.9411 and 0.9440, and t-statistic values as low as 3.2119 and 3.6541, respectively. The binary BRL model constructed using minute-scale training data follows closely, while the model constructed using hourly data shows the least satisfactory performance among the three.

When validated using the July hourly radiation data from the Kazuo station, the binary BRL model constructed using hourly training data demonstrates the best overall performance, while validated using the December 10-min data from the Kazuo station, the binary BRL model constructed using 10-min training data demonstrates the best overall performance, with correlation coefficients of 0.9257 and 0.9603, t-statistic values of 5.7016 and 22.7655, respectively. For the July hourly validation data from the Guangzhou station, the maximum correlation coefficient is only 0.8807, indicating an overall underestimation of the calculated diffuse fraction compared to the measured values. The correlation coefficient for the December 10-min validation data is above 0.9, however the corresponding t-statistic is relatively large. The validation results for the hourly and 10-min data from the Wuhan station are generally good, with correlation coefficients above 0.9. The hourly validation data from the Urumqi station shows similar performance to the hourly data from the Guangzhou station, with the maximum correlation coefficient at 0.8575. However, the 10-min data from the Urumqi station for validation exhibit correlation coefficients above 0.9, with smaller mean bias and root mean squared error compared to the 10-min data, but the t-statistic is relatively larger.

In conclusion, the binary BRL model performs well when validated with data from the Lijiang station, while the calculated diffuse fraction for data from the Kazuo, Guangzhou, Wuhan, and Urumqi stations is relatively less accurate, with MBD values generally below 0, indicating an overall underfitting state. This is attributed to the training dataset being exclusive to radiation data from the Lijiang station. Therefore, validation using data from the Lijiang station is more effective in evaluating the model’s performance. For future research, it is necessary to include cross-regional station data in the training dataset to enhance the model’s applicability.

Using the validation data from the Lijiang and Kazuo stations as examples, Figure 2 illustrate the comparative validation plots of different-parameter binary BRL models based on the aforementioned station validation datasets. The three sets of different parameters for the binary BRL model correspond to those obtained from hourly, 10-min, and minute-scale radiation data training.

Figures 2A–F present the validation comparison plots of the diffuse fraction calculated by the different-parameter binary BRL models using the Lijiang station’s hourly and 10-min radiation data for validation. Specifically, the binary BRL models trained with 10-min and minute-scale data yield better representations of the diffuse fraction distribution. Figures 2G–L show the validation comparison plots of the diffuse fraction calculated by the different-parameter binary BRL models using the Kazuo station’s hourly radiation data in July and 10-min radiation data in December, respectively. From the figures, it is evident that the diffuse fraction calculated by the binary BRL models generally falls below the measured diffuse fraction, indicating an underfitting state.

The influence of the diffuse fraction on sunshine duration measurement was investigated in this study. Based on atmospheric radiation transfer theory and utilizing Taylor series expansion, a binary BRL model was established to relate the diffuse fraction to the clear-sky index and atmospheric optical thickness. Training datasets were created using radiation measurement data at different time scales, namely hourly, 10-min, and minute-scale, from the Lijiang station. Through non-linear fitting in Matlab, the binary BRL model was developed to estimate the diffuse fraction with respect to the clearness index and atmospheric optical thickness, providing corresponding model parameters and their corresponding two-dimensional surfaces. Validation datasets from two distinct stations, namely the Lijiang and Kazuo stations, were selected to calculate the diffuse fraction using the binary BRL model. Statistical analyses were performed to assess the correlation coefficient, mean bias deviation (MBD), root mean squared deviation (RMSD), and t-statistic between the measured and computed diffuse fraction values.

The conclusions and prospects are as follows:

(1) The diffuse-fraction values computed by the binary BRL model demonstrated correlation coefficients above 0.8, MBD within ±0.2, RMSD within 0.25, and minimum t-statistic as low as 0.1074. This indicates the suitability of the model for constructing the sky diffuse fraction and provides new insights into the study of the distribution pattern of diffuse radiation and its impact on sunshine duration measurement.