Numerical Assessment of Downward Incoming Solar Irradiance in Smoke Influenced Regions—A Case Study in Brazilian Amazon and Cerrado

Smoke aerosol plumes generated during the biomass burning season in Brazil suffer long-range transport, resulting in large aerosol optical depths over an extensive domain. As a consequence, downward surface solar irradiance, and in particular the direct component, can be significantly reduced. Accurate solar energy assessments considering the radiative contribution of biomass burning aerosols are required to support Brazil’s solar power sector. This work presents the 2nd generation of the radiative transfer model BRASIL-SR, developed to improve the aerosol representation and reduce the uncertainties in surface solar irradiance estimates in cloudless hazy conditions and clean conditions. Two numerical experiments allowed to assess the model’s skill using observational or regional MERRA-2 reanalysis AOD data in a region frequently affected by smoke. Four ground measurement sites provided data for the model output validation. Results for DNI obtained using δ-Eddington scaling and without scaling are compared, with the latter presenting the best skill in all sites and for both experiments. An increase in the relative error of DNI results obtained with δ-Eddington optical depth scaling as AOD increases is evidenced. For DNI, MBD deviations ranged from −2.3 to −0.5%, RMSD between 2.3 and 4.7% and OVER between 0 and 5.3% when using in-situ AOD data. Overall, our results indicate a good skill of BRASIL-SR for the estimation of both GHI and DNI.


Introduction
Brazil has a vast solar energy resource [1][2][3] and has experienced a boost in photovoltaic deployment in recent years due to government incentives and technological advances [4,5]. Several studies have shown that solar energy, as part of the diversification of the renewable energy mix, could be decisive to increase energy security, counter-balancing the vulnerability imposed by the high dependency on the hydro-power [6][7][8][9]. In particular, concentrating solar power (CSP) technologies have shown a noteworthy potential for Brazil in scenarios of climate change mitigation ( [10][11][12][13]), especially as a complementary heat supply for industrial processes or hybrid power generation [14,15]. It should be noted, however, that some potential areas for CSP development, like the Central-West and the Southeast regions, are often affected by biomass burning haze during the dry season as a result of long-range transport [16][17][18][19].
Atmospheric aerosols are the most important factor for solar radiation extinction in cloudless conditions, followed by water vapor. In particular, the direct normal irradiance (DNI) is 2-4 times more sensitive to the presence of aerosol than the global horizontal irradiance (GHI) [20]. The impact of dust aerosols on the DNI has been assessed for several arid and semi-arid sites [21][22][23]. Meanwhile, similar assessments taking into account the biomass burning aerosols are lacking. Although lower in magnitude than the dust aerosols impact due to the comparatively moderate AOD values, the impact of biomass burning aerosols in the DNI is unlikely to be negligible in regions where burning activity is seasonally intense (in number and burned area). The large loads of aerosols typically injected into the atmosphere during the dry season in Brazil can result in high aerosol optical depths (AOD). During this period of the year, in Central Amazon, values of AOD at 500 nm between 0.75 and 1.0 are frequently found, while in the Southern Amazon, values of AOD at 500 nm well above 1.0 are common, and values up to 5.0 have been reported in years with more intense biomass burning activity [24].
An accurate evaluation of the downward solar irradiance at the surface, chiefly the DNI, is thus essential for analysing of the bankability of large-scale CSP projects in regions affected by biomass burning plumes and haze. Long-term quality-assured ground measurements would be the most reliable way to identify and assess CSP and photovoltaic (PV) solar power installations prospect sites. However, surface solar irradiances ground data is only measured at a few locations sparsely and heterogeneously distributed, and measurements of the direct component are even more scarce globally and in Brazil. In this context, numerical models came into action as necessary and powerful tools to improve the spatial monitoring of solar irradiance. In the energy sector, numerical models can support bankability analysis of solar power projects and short-term forecasting for daily plant operation [25]. Other knowledge areas, like climate and weather prediction and agrometeorology, can benefit from high spatial resolution and accurate solar irradiance estimates provided by radiative transfer models.
The Brazilian Atlas for Solar Energy [26,27] has been used as a reference document to support the Brazilian solar power sector and energy planning. The BRASIL-SR [27][28][29][30] is a physically based radiative transfer code that follows a two-stream approximation with δ-Eddington scaling, and was used to provide the surface downward solar irradiance solar resource mapping for both editions. Despite the high confidence (low uncertainties) of previous solar resource mapping conducted with this model, earlier versions of BRASIL-SR used an aerosol parameterization based on meteorological horizontal visibility range [30] that was not able to adequately represent the actual atmospheric load of aerosol from biomass burning, often concentrated in higher atmospheric levels [31,32]. In such conditions, the model systematically overestimated the clear-sky downward solar irradiance at the surface, especially the DNI, increasing uncertainties during the dry season.
In this work, a new version of the BRASIL-SR clear-sky model was developed to improve the representation of aerosol radiative attenuation, and reduce the uncertainties of the surface solar irradiance estimates in cloudless hazy conditions and clean conditions. The study area comprises the Amazon and Cerrado biomes, regions which experience large biomass fire activity during the dry season but hold distinct amounts of annual precipitation and water vapour content in the column.
The paper is organized as follows: the model BRASIL-SR is described in Section 2, and Seciton 3 describes the available database and methodology used to improve and evaluate the solar irradiance estimates provided by the model. Seciton 4 discuss the model outputs for sites located in the Amazon and Cerrado regions of Brazil. Finally, Seciton 5 includes the discussion of results, conclusions and future research.

BRASIL-SR Model
BRASIL-SR [27][28][29][30] is a satellite-based model that estimates the downward surface solar irradiance. The core of the BRASIL-SR is a physically-based radiative transfer model that is executed for two atmospheric conditions: cloudless (although a high aerosol optical depth is allowed) and overcast with a very high cloud optical depth. Then, the solar irradiance components at the surface for any cloud cover conditions are obtained by the interpolation between both solutions using the effective cloud cover index obtained from visible satellite imagery. The model runs using time steps that are typically coincident with satellite image times, but this is not a constraint for studies in a clear-sky atmosphere.
For clear sky assessments, the model BRASIL-SR requires the following regional input data for each grid cell: longitude, latitude, altitude, surface temperature, relative humidity, total precipitable water vapor (PWV), total ozone in the column (O 3 ), AOD in 550 nm and Angström's exponent (AE), biome classification, and the Moderate Resolution Imaging Spectroradiometer (MODIS) bi-directional reflectance distribution functions (BRDF) kernel parameters. Additionally, local data from observations (e.g., from the Aerosol Robotic Network-AERONET [33]) of PWV, O 3 , AOD, and AE can be entered as input for a particular grid cell, overriding regional data.
The radiative transfer calculations follow a two-stream approximation with δ-Eddington scaling to estimate the downward surface GHI. The DNI is obtained with and without δ-scaling. Likewise, the diffuse horizontal irradiation (DHI) can be obtained as GHI, considering δ-scaling, or from the difference between GHI and the direct horizontal irradiation obtained without δ-scaling. There is no estimation of the circumsolar irradiation. The model uses 37 spectral intervals distributed from 200 to 3500 nm, with narrower intervals in the UV and visible spectral ranges and wider in the near-infrared region. Absorption of solar radiation by water vapor, O 3 , CO 2 , and oxygen is calculated using exponential-sum fits to transmittances (ESFT). The current model version uses updated absorption coefficients for water vapor calculated based on the method proposed by Wiscombe and Evans [34] for convergence of ESFT and transmittance data obtained with Py4CAtS-PYthon [35] and HITRAN2016 spectroscopic database [36]. The solar spectral data at the top of the atmosphere (TOA) follows Gueymard [37]. The vertical atmospheric profiles of air temperature, atmospheric pressure, and gases follow Anderson et al. [38], and are available for the five standard atmospheres, to be selected according to the surface temperature of the grid cell.
The surface boundary condition uses the ground spectral albedo for direct and diffuse radiation obtained using the BRDF kernel parameters and polynomial formulae with coefficients that are wavelength-independent, as described by Schaaf et al. [39]. The BRDF kernel parameters depend only on the wavelength and the soil and vegetation characteristics in the grid cell. BRDF surface kernel parameters derived for MODIS are linearly interpolated to the 37 spectral wavelengths used in the model BRASIL-SR and later used to calculate the spectral albedos for each wavelength.
Aerosol optical depth at 550 nm is vertically distributed for each grid cell and interpolated for the remaining 36 spectral interval wavelengths using the Angström exponent entered as input, which is assumed to be applicable for all spectral intervals. BRASIL-SR uses a fixed aerosol profile for heights between 5 km and 50 km with the aerosol optical depth within these heights given by τ 5−50km = 0.0216. Below 5 km-height, two heightexponentially decaying profile options are available to the user: (a) the maximum aerosol extinction coefficient occurs in the first km above the surface (the first atmospheric layer), or (b) the maximum aerosol extinction coefficient occurs in the second layer (2 km above the surface) with the first (1 km) and third (3 km) layers presenting equal extinction coefficients. The second option simulates the typical aerosol profile observed during the biomass burning season [32], and it is selected for this study.
The selection of aerosol optical properties is based on the prevailing biome at each grid cell. The optical properties for the Brazilian Amazonia and Cerrado regions were obtained as described in Appendix A. The model BRASIL-SR assumes a default set of optical properties corresponding to the polluted continental aerosol for other biomes and locations, although additional sets of optical properties can be easily introduced in the future. Local data of precipitable water vapor, when available, is corrected to the model base height using the scale height for the profile. Afterward, the ratio between local (height-corrected) and calculated precipitable water vapor is used to correct the water vapor profile.
The broadband outputs of GHI, DNI, and DHI are available for the whole area (all domain grid points) established in the input dataset. The BRASIL-SR outputs include the spectral downward surface solar radiation for specific locations defined by the user.
BRASIL-SR is less complex and demands less computational resources than models like LibRadtran [40] and DISORT [41]. However, it demands more computational resources than broadband site-specific models. Memory requirements can be particularly demanding. Although BRASIL-SR radiative transfer calculations are typically run for the time of satellite images, we expect to tackle this limitation in future versions of BRASIL-SR.

Experimental Data and Methods
As mentioned, the current study focuses on improving the representation of the impact of biomass burning aerosols on the downward surface solar irradiance provided by the model BRASIL-SR. Therefore, the cloudless condition is mandatory to isolate the aerosol radiative effect. Ground measurement sites operating in the Amazon and central region of Brazil were carefully selected to manage the input database required to feed the model. These regions are particularly interesting due to the high number of fire spots within their domain associated with biomass burning during the dry season (from May to October), and to a greater extent within the biomass burning season (from August to October).
Although several AERONET stations in the Amazon and Cerrado regions make available long time-series of high-quality data, only a few have co-located instruments for solarimetric observations, particularly for DNI data acquisition. Ground data of incoming solar irradiance is fundamental to validate the model outputs and quantify their deviations. The Green Ocean Amazon (GOAmazon) [42] was a successful experiment and provided the scientific community with good quality solarimetric data and AOD acquired by colocated AERONET sites. The GOAmazon extended through the wet and dry seasons from January 2014 through December 2015. Unfortunately, during the biomass burning season of 2014, there was no AOD data from level 2.0 in the AERONET database for either ARM_Manacapuru or Brasilia_SONDA. On the other hand, 2015 was a dryer and warmer year in the Amazon region, associated with a strong El Niño-Southern Oscillation (ENSO) activity [43], and biomass burning activity observed in the Amazon and Central regions of Brazil was more intense than in previous years [44,45]. Considering both the data availability and aerosol load, the period from July to December 2015 was selected for the study. It includes the months with significant aerosol load in the region of study as well as higher density of observations, as shown in Figure 1 (top panel and bottom panel, respectively).

Observational Data
Four sites provided the observational data used for this work: ARM_Manacapuru, Manaus_EMBRAPA, Brasilia_SONDA and Palmas_SONDA. ARM_Manacapuru, Man-aus_EMBRAPA and Brasilia_SONDA have co-located AERONET stations, with the same name, that provide level 2.0 AOD , PWV, and O 3 column content data. Additionally, spectral irradiance data were available from multifilter rotating shadowband radiometers (MFRSR), operating at Manaus_EMBRAPA and ARM_Manacapuru. These two sites were part of GoAmazon experiment under the classification of time point zero (T0e) and three (T3), respectively. While the former one was installed and operated since february of 2011 [46] at an unity of the Brazilian Agricultural Research Corporation (EMBRAPA), the later was specifically deployed for the GoAmazon experiment. Brasilia_SONDA, and Palmas_SONDA sites belong to the Brazilian Environmental Data Organization System (SONDA) network of measurement sites distributed in the Brazilian territory aiming to provide reliable meteorological and solarimetric data to support the Brazilian renewable energy sector [47,48]. There is no aerosol data acquisition system in Palmas_SONDA.
The Data Availability section at the end of the paper summarizes how to access all data sources used.
The top panel of Figure 1 displays the time-series of AOD 550nm for the year 2015, indicating the selected period of study (from July to December) with a thick black dashed line rectangle. The inset Figure shows the probability distribution for AOD 550nm for each AERONET station for the selected period of study. It can be seen that the highest values of AOD 550nm happened in October for all three sites, although high values can be found in November as well for ARM_Manacapuru and Manaus_EMBRAPA stations. The highest probabilities are found for the 0.  Measurements of GHI, DHI, and DNI at ARM_Manacapuru were collected by SKYRAD instrument [49] at 1 min frequency. MFRSR spectral irradiance data sampled with a frequency of 20 s at ARM_Manacapuru is also available at the ARM site [50]. The broadband and solar spectral data acquired at ARM_Manacapuru are open-access at the Atmospheric Radiation Measurement (ARM) website [51].
At Brasilia_SONDA site, DNI observations were stored at 1 min frequency using a NIP Eppley pyrheliometer until 18 September 2015, when it was replaced by a Kipp&Zonen CHP 1. The GHI and DHI data were obtained using a pyranometer CM22 Kipp&Zonen and automatic sun tracker 2AP BD. At Palmas_SONDA, the GHI and DHI data were obtained using a CM11 Kipp&Zonen pyranometer combined with shadowing ring CM121B + adapter CV2. GHI and DHI at Brasilia_SONDA and Palmas_SONDA are also stored at 1 min frequency. The DNI was calculated from (GHI − DH I) and using the solar zenith co-sine correction. The solar radiation data acquired at Brasilia_SONDA and Palmas_SONDA is available in open access mode at the SONDA network website ( [48]). The data quality control procedure follows the criteria used by WMO for the Baseline Solar Radiation Network (BSRN) [47,52].
The GHI, DNI, and DHI data at Manaus_EMBRAPA were derived from MFRSR data ( [53]) since there was no pyrheliometer or pyranometer on site. Spectral irradiance at Manaus_EMBRAPA, sampled every 60 s, was also obtained from the MFRSR [53].
Additionally to the quality check adopted by the solar data owners, the following constrains were imposed on the records for downward surface solar irradiance [54,55]: GH I/GH I 1 = 1.0 ± 0.08 for (DHI + DN I * cos(SZA)) > 50 Wm −2 , SZA < 75 • (1) where GH I 1 is DH I + DN I * cos(SZA) and SZA stands for solar zenith angle. GHI, DNI, hemispheric spectral irradiance and direct normal spectral irradiance data were compared to BRASIL-SR output data, mainly for consistency check of the spectral result of the model. An example of BRASIL-SR spectral output and how it compares with observations can be found in Appendix B. Table 1 summarize the location of each site, as well as the quantities measured and instruments used in each site. On field, for a regularly maintained sensor, the uncertainty of DNI measurements is about 2.0%, while for GHI and DHI the uncertainty is about 3.0% for large irradiance signals, but could be as large as 10% for low solar elevation angles [25]. In the case of measurements with a rotating shadowband radiometer (of which MFRSR are a special case), the uncertainty of field GHI and DHI measurements are in the order of 4%, respectively, while for the DNI is around 5% [25]. In addition, when only two components are measured (like in MFRSR measurements and in Palmas_SONDA site), some closure tests that could further reduce the uncertainty can not be applied [25].

BRASIL-SR Input Data and Configuration
Two configurations of the model BRASIL-SR were used: (a) the in-situ experiment used to assess GHI and DNI in specific locations shown in Figure 2 using local measurements for aerosol data (AOD and Angström exponent), precipitable water vapor, and ozone; and (b) the regional experiment, used to assess GHI and DNI in the regional domain (see Figure 2) using aerosol, precipitable water vapor, and ozone data provided by the reanalysis database. In both configurations, the time step was 5 min. In this study, there was no use of site-adaptation techniques. The former uses the local observation dataset available for the locations with operational AERONET data acquisition system. Details on the availability of AERONET data and characteristic values for the study period are described in Section 3.1. In all AERONET sites, AOD 500nm is interpolated to 550 nm using the Angström exponent for 500-675 nm range and the Angström law. AOD 550nm , Angström exponent for the range 440-870 nm, precipitable water vapor, and ozone are linearly interpolated to each time step. It is expected that the in-situ experiment is representative of the actual skill of the model, although some uncertainty is introduced in the time-interpolation of the AERONET dataset, particularly if the number of data points for a given day is low. To minimize such uncertainty impact, the in-situ experiments were performed only for days presenting at least five valid AERONET measurements.
The configuration used in the regional experiment assumes that the local observation dataset is not available. The BRASIL-SR runs for a spatial domain presenting 533 × 299 grid cells over the Brazilian Amazon and central regions with latitudes ranging from −2.0 • to −17.0 • and longitudes ranging from −62.0 • to −46.0 • , as depicted in Figure 2. Amazonia and Cerrado are the predominant biomes in the domain, except for the southwest corner, where the Pantanal biome (located in the western border of Brazil, between Bolivia, Paraguay and the Cerrado region) and a small region in Bolivia are present ( Figure 2). Monthly mean surface BRDF kernel parameters derived for MODIS [56] were used to derive the spectral surface albedo. Biome data is obtained using the Brazilian Institute of Geography and Statistics (IBGE) biome map 1:250,000 [57].
Regional AOD and Angström exponent data were obtained from MERRA-2 Reanalysis hourly M2T1NXAER collection (variables TOTEXTTAU and TOTANGSTR) [58,59], linearly interpolated to each BRASIL-SR time step. The ozone profile from Anderson et al. [38] was corrected using monthly data from the Copernicus Climate Change Service [60]. Global Forecast System (GFS) reanalysis [61] provides surface temperature and relative humidity, and column precipitable water vapor.

Model Validation
The following statistical indices, in the form proposed by Gueymard [62], were considered to evaluate the model performance in assessing the clear-sky downward surface irradiance at the four ground measurement sites described earlier: where p i and o i are the model outputs and observed values, respectively, O m is the mean observed value and N cs is the number of clear sky observations. • Root Mean Square Difference (RMSD) • Kolmogorov-Smirnov integral (KSI)-defined as the integrated differences between the cumulative distribution functions (CDFs) of the model (p) and observational (o) data sets.
where D n is the absolute difference between the two normalized distributions within irradiance interval, h min and h max are the minimum and maximum solar irradiance values, and Ac is a characteristic quantity of the distribution • OVER-statistical parameter similar to KSI, but the integration is calculated only for those CDFs' differences that exceed the critical limit ( 1.63

√ N cs
) of the Kolmogorov-Smirnov method. From Equation (7), it can be notice that OVER is 0 (zero) if the CDFs' differences always remains below the critical value [62].
In addition to the comparison with observational data, BRASIL-SR results were compared with two broadband clear sky models, McClear [55,63] and REST2 [64,65].
McClear, version 3.1, 1-minute data were downloaded from the Copernicus Atmosphere Monitoring Service (CAMS) [66]). In this model, a combination of a look-up table of pre-calculated radiative transfer results and an interpolation scheme is used. The aerosol input comes from MACC reanalysis AOD data. The aerosol type is derived from the partial optical depths of the MACC aerosol species, while Angström exponent is calculated from MACC AOD at 550 and 1240 nm. All MACC variables are corrected for altitude before entering the McClear model. REST2, version 5, 1-min clear-sky time series were obtained through the use of irradpy Python library, as described by Bright et al. [65]. REST2 v5 input MERRA-2 total AOD and Angström exponent are also downloaded and interpolated to the 1-minute time step by the use of the same library irradpy. REST2v5 is based on the spectral model SMARTS [67], from which two sets of transmittance parameterizations for broad spectral bands 0.29-0.70 µm and 0.70-4.0 µm were derived.

Selection of Clear-Sky Periods
In this work, we assumed clear-sky conditions as cloudless conditions, but high load of aerosol pollutions scenarios were allowed. The methods developed by Bright et al. [68] and Inman et al. [69] for clear-sky detection (hereafter, Bright-Sun and Inman15, respectively) were considered to classify sky condition as cloudless or cloudy based on local solar data observations. For the clear-Sun disk situations, i.e. those cases when there is a clear line of sight to the Sun, Ineichen et al. [70] (hereafter, Ineichen06), and Ineichen et al. [71] (hereafter, Ineichen09) methodologies were also tested. Inman15, Ineichen06, and Ineichen09 received a high score compared to other robust methods investigated by Gueymard et al. [72]. On the other hand, the more recent Bright-Sun algorithm makes use of the previous experience in clear-sky detection to propose an efficient and widely applicable algorithm.
The implementations of the four methods were made on Python, based on Bright [73]. For those methods that required clear-sky irradiances, time series from McClear clear-sky model [55,63] were used. McClear 1 min-resolution data can be obtained from the Copernicus Atmosphere Monitoring Service (CAMS) [66]). The solar zenith angles SZA > 85 • were discarded to avoid uncertainties regarding solar data acquisition at low solar elevation angles.
The Inman15 method is based on five criteria that are applied on ten minutes sliding windows Reno et al. [74]. However, Reno et al. [74] relaxed the threshold values for GHI and added new thresholds for DNI. If the imposed criteria are all met within a given window for GHI and DNI, then all data records belonging to that window are qualified as clear-sky [69].
Ineichen06 method assumes two criteria based on the direct solar irradiance to screening the clear-sky condition: (a) a DNI threshold-DN I > 0.9 * ETH/cos(SZA) * exp[−2 * m/(9.4 + 0.9 * m)], where the right side represents the DNI in cloudless condition DN I cs , ETH is the the extraterrestrial solar irradiance, and m is the air mass; and (b) a DNI variability threshold-the DNI variability remains ±10% of DN I cs variability.
The method Ineichen09 assumes a minimum and maximum threshold for the modified clearness index-above 0.65 and below 1.0 are classified as clear-sky. The method performs poorly for detecting clear-sky situations but presents a remarkably high score for the classification of clear-Sun disk situations (i.e., cases in which there is a clear sight of the sun) even in polluted conditions [72]. However, a high tendency to give false positives was observed.
Bright-Sun algorithm starts with preliminary results from Reno model [74], that are used in an optimization of the GHI, DNI, and DHI clear-sky curves. Afterward, a tricomponent analysis that uses normalized thresholds for the three variables is performed, followed by the verification of two duration criteria: a maximum of 9 periods of cloudiness in 90 min-window and no cloudiness in a 30 min-window. This somewhat complex methodology gave overall better results than earlier methodologies, as showed in Bright et al. [68]. However, in our study clear-sky screening using Bright-Sun occasionally failed if the clear-sky curve presented a different shape compared with the curve of observational data. Such false results could happen, for instance, when MACC reanalysis AOD (used in the obtaining of McClear clear sky time series) varied through the day in a way that did not match AERONET observations. Table 2 summarize the percentage of observational solar data records classified as clear-sky/clear-Sun. It can be seen that the Bright-Sun delivered the lowest percentage of observations flagged as clear-sky. Inman15 and Ineichen06 detected similar quantity, around 60% lower than the Ineichen09 clear-sky identifications.   Regarding the study period, a similar behavior was observed. Remarkably, even when Bright-Sun reported the lowest percentage of clear-sky periods, it also appears to have the lowest number of false positives and false negatives. Ineichen06 typically gave false negatives for low SZA in polluted conditions, most likely due to a violation of the first condition for higher AOD values, as well as false positives for high SZA, when the imposed conditions appeared to not being stringent enough. Ineichen09 showed the highest rate of detection for true positives but also gave a significant number of false positives and false negatives, as illustrated in Figure 3 and already reported by Gueymard et al. [72]. Finally, Inman15 showed a tendency to give false negatives at higher SZA. In highly polluted conditions (here arbitrarily defined as AOD 500nm ≥ 0.5), all methods frequently failed to detect apparent true positives, even in cases when BRASIL-SR clear-sky irradiances were very close to observation. Such fails could be partially attributable to inaccuracies in the McClear clear-sky irradiances propagated from the MACC AOD data.
Although an all-sky camera would be desirable for a more accurate assessment of false positive/negative occurrences, the comparison with solar irradiances provided by clear-sky models using local input data can indicate possible false negatives. For instance, Figure 4 shows results for September 10th, at ARM_Manacapuru, when the BRASIL-SR in-situ experiment displays a remarkable resemblance with observations from 10:30 UTC and up to 15:00 UTC. On the other hand, the Bright-Sun algorithm flagged data as clear-sky only a short period close to 21:00 UTC, when the McClear outputs are close to observations. It can be seen that, given the shape of the clear-sky curve, a single coefficient of optimization is unlikely to improve the clear-sky curve for the whole day.
In general, Bright-Sun performed best among the four tested algorithms, with a better skill for detecting clear-sky and clear-Sun even in moderately polluted situations. Therefore, this was the method used for both clear-sky and clear-Sun screening procedures.

Global Horizontal Irradiation
The validation procedure for BRASIL-SR outputs was accomplished using observational data at the time of model results. Thus, only observations with available BRASIL-SR results and classified as clear sky were considered. Figure 5 shows the scatter plots between observed and model-estimated global horizontal irradiation for the four sites and the two experiments: in situ and regional. In general, a noteworthy overall correspondence between BRASIL-SR outputs and observations, with data pairs aligned close to the 1:1 line. Table 3 shows the benchmarking results for GHI data. There is no local data for AOD and Angström exponent in Palmas_SONDA, so only the regional experiment results are shown. The relative bias (MBD) was low for all sites and for the two model experiments, being lower (−0.3% < MBD < 1.3%) for the in-situ experiment, as expected. The MBD ranged from −1.0% to 2.1% for the regional experiment.  Figure 5. Scatter plot comparing observational GHI data and BRASIL-SR GHI outputs running the regional (grey dots) and in-situ experiments (red dots) for the four ground measurement sites. The root mean square deviation were also lower for the in-situ experiment with values between 2.1% and 2.9%, while 3.2% ≤ RMSD ≤ 5.7% for regional experiment. MAD values were also higher for the regional experiment, with values between 2.6% and 4.4%, while in the in-situ experiment, the lowest value was 1.6% in Brasilia_SONDA, and similar values were obtained in ARM_Manacapuru and Manaus_EMBRAPA (2.1% and 2.0%, respectively). The low values obtained for KSI and OVER parameters indicate the high similarity between ob-served and modeled cumulative distribution functions (CDFs), except for Palmas_SONDA, where the value of OVER was larger than zero in the regional experiment, although still presenting a very low value (1.9%).
As discussed in Section 3.2, the in-situ experiment uses interpolated aerosol data based on at least five AERONET data records, from level 2.0, acquired in the same day. This procedure is less stringent than considering a restricted time difference (typically a few minutes) between AERONET data and irradiance observations, as used in other studies [21,22,75]. However, due to the low availability of co-located AERONET and DNI observational data in clear sky conditions for the region of study, using only quasi-simultaneous data would imply in much lower number of data records to compare and validate the model results, compromising the robustness of the statistic analysis performed.
Statistics for McClear and REST2 running for the same clear-sky dataset are presented for benchmark comparisons. McClear and REST2 models had similar performance in ARM_Manacapuru and Manaus_EMBRAPA. The REST2 performed better in Brasilia, while McClear was better in Palmas_SONDA. The statistical parameters show that the BRASIL-SR regional experiment provided the best performance, except for Palmas_SONDA, where McClear delivered the lowest deviations. Figure 6a shows the shape of the deviation distribution for ARM_Manacapuru using violin plots layered with strip plots. The violin plot is similar to a box plot, with the addition of a rotated kernel density plot on each side. The vertical line shows the median and interquartile range (IQR). The strip plot complements the violin plot, showing all observations along with some representation of the underlying distribution. Since measurement uncertainties increase for high SZA and model deviations are amplified for low irradiance values, a threshold SZA ≤ 75 • was used.
As expected, the best results were obtained when using interpolated local AERONET data (in-situ experiment), with a narrower IQR and median value very close to zero. Nevertheless, a more fair comparison can be made among the three other cases where all models used reanalysis data as input for AOD and Angström exponent. All three presented similar distribution shapes for Manacapuru site. The deviations of the BRASIL-SR outputs in the regional experiment displayed a median value closer to zero than McClear and REST2. BRASIL-SR also achieved the narrowest IQR and distribution shape. Some outliers, however, can be noticed with deviation above 15%. Among the three models, REST2 displayed the highest number of outliers.
For Manaus_EMBRAPA, BRASIL-SR results in both in-situ and regional experiments displayed a positive bias that is, however, smaller than the bias obtained by McClear and REST2 (Figure 6b). The medians of model output deviations were higher than MBD values in Table 3, and the three models overestimated the observed values. However, irradiance data for Manaus_EMBRAPA was acquired with an MFRSR on field, which presents somewhat higher uncertainties (see Section 3.1 for details).
The error distributions at Brasilia_SONDA also displayed the lowest medians and narrowest IQR for BRASIL-SR results as shown in Figure 6c). While McClear results were symmetric, all other three presented a widening in the direction of positive errors, with the highest number of outliers by BRASIL-SR regional experiment and REST2.  Figure 7 shows that the highest deviations between modeled and observational GHI occurred in Palmas_SONDA (the reader should note that the vertical scale used is different from the one used in Figure 6, extending up to 80%). The best GHI estimates were obtained with McClear, displaying median values close to zero. BRASIL-SR regional experiment results presented a negative median and a significant number of positive outliers presenting errors above 20%. On the other hand, REST2 presented a positive bias and even higher GHI errors. The violin plots are consistent with MBD values presented in Table 3. Given the similarities between the three cases using reanalysis in the other three locations, the results for Palmas_SONDA suggest that aerosol data from MACC reanalysis gives more consistent results than MERRA-2 for this location. However, particularities of each model are also expected to influence.

Direct Normal Irradiation
The two-stream approximation with δ-Eddington scaling has been widely used to solve the radiative transfer equation in several numerical models because its computational time demand is low and yet reasonably accurate, at least for estimating global horizontal downward surface irradiance. However, it may also lead to an overestimation of the direct normal surface irradiance. Mathematically, this is expected because the scaled optical depth is always smaller than the actual optical depth. However, the first studies on the two-stream method suggested that the neglected circumsolar irradiance could compensate for the DNI overestimation [76].
In this study, results for DNI obtained using δ-Eddington scaling and without scaling are compared. Figure 8 displays the scatter plots between observed and estimated DNI obtained with the traditional use of δ-Eddington scaling the optical depth for the four sites and the two experiments performed. Again, for the Palmas_SONDA site, results are available only for the regional experiment. The corresponding statistics are presented in Table 4. A positive bias is apparent for all cases in Figure 8, as confirmed by MBD values that range from 4.0 to 5.4% in the in-situ experiment and between 2.5 and 8.4% in the regional experiment.  Figure 8. Scatter plot for DNI between observational data and BRASIL-SR (obtained using δ-Eddington scaling) for regional experiment (grey) and in-situ experiment (red) for the four ground measurement sites. On the other hand, Figure 9 gathers the scatter plots between observed and estimated DNI obtained without scaling the optical depth for the four sites and two experiments considered. The corresponding performance parameters are presented in Table 5. While the positive bias decreased, becoming slightly negative in the in-situ experiment (expected result since BRASIL-SR does not consider the circumsolar irradiance), the MBD got similar values or even slightly higher than in the regional experiment.  The MAD is closer to zero, and the RMSD and KSI values are lower than the results achieved in the scaled AOD case. Remarkably, not using scaling in the optical depth reduced not only the MBD and MAD but also the RMSD in all cases for both experiments. KSI and OVER parameters were also significantly improved, with shallow values (OVER ≤ 6.2) for all cases, except for Palmas_SONDA in the regional experiment, where the highest OVER was obtained (OVER = 26.9), but still, this value is lower than the OVER = 91.9 value obtained using δ-Eddington scaling for the same site. It should be considered that uncertainties in the observed DNI are higher in Manaus_EMBRAPA and Palmas_SONDA than in ARM_Manacapuru and Brasilia_SONDA since DNI is not obtained independently in the former sites, but from the difference between global and diffuse horizontal irradiance using a MFRSR and a CM11 Kipp&Zonen pyranometer, respectively [77].
The violin plots of the DNI error obtained with and without δ-Eddington scaling are displayed in Figure 10, layered with strip plots. Again, in this Figure and all following DNI error figures, we used a threshold SZA ≤ 75 • . In both ARM_Manacapuru and Brasilia_SONDA, the improvement of the results for DNI is evident when not using optical depth scaling, reducing the width of the distributions and IQRs, and bringing median values close to zero. In Manaus_EMBRAPA, again, all performance metrics improved without δ-Eddington scaling. In addition, the number and magnitude of outliers were reduced. However, a widening of the output distribution towards negative DNI errors and an apparent increase in the number of negative outliers is observed without δ-Eddington scaling. As previously discussed, at least part of this behavior could be attributed to the higher uncertainties in MFRSR irradiances. The variability of DNI error with AOD is represented in Figure 11, where a threshold of SZA ≤ 75 • was considered. Results presented in Figure 11 should be analyzed considering that higher AOD values are less frequent, as shown in the histograms of Figure 1, inset panel. Situations with high AOD but apparent clear sky conditions were also frequently falsely flagged as cloudy (see Section 3.4 for details). Therefore, the variability expressed by the box-and-whisker plot for higher AOD values can be misleading due to the lower number of samples. The median error in ARM_Manacapuru was 13.8% with δ-Eddington scaling and −3.7% without scaling at AOD 550nm = 0.7 (71 samples were included in the 0.6-0.7 AOD 550nm bin; a total of 42 data points presented AOD 550nm > 0.7, while 2223 data points had AOD 550nm < 0.6). In Manaus_EMBRAPA, the median error was 11.0% with δ-Eddington scaling and −7.3% without scaling at AOD 550nm = 0.7 (37 samples were included in the 0.6-0.7 AOD 550nm bin; a total of 7 data points presented AOD 550nm > 0.7, while 1462 data points had AOD 550nm < 0.6). In Brasilia_SONDA, the median error was 13.9% with δ-Eddington scaling and −1.7% without scaling at AOD 550nm = 0.5 (26 samples were included in the 0.4-0.5 AOD 550nm bin; no data points presented AOD 550nm > 0.5 for the in-situ experiment and clear sky conditions, and 2108 data points presented AOD 550nm < 0.4). In summary, an increase in the relative error of DNI results obtained with δ-Eddington optical depth scaling as AOD increases is evidenced for the three sites. A dependency somewhat weaker and opposite in sign was also obtained for results with no scaling (Figure 11, panels d-f). In Table 5 are also presented, for comparison, the performance parameters for DNI outputs from McClear and REST2 for the same data points used to run BRASIL-SR. Results with BRASIL-SR in-situ experiment displayed the best metrics, thus emphasizing again the impact of using more accurate aerosol input data. BRASIL-SR regional experiment results for ARM_Manacapuru and Manaus_EMBRAPA were similar to those obtained with McClear in terms of MBD, RMSD, and MAD, but presenting lower values for KSI and OVER. As for GHI in Palmas_SONDA, McClear displayed the best performance, while BRASIL-SR regional results were more accurate in Brasilia_SONDA site. REST2 presented the poorest performance in all sites, except for Brasilia_SONDA where it presented a better performance than McClear.
The shape of the DNI error distribution for ARM_Manacapuru is displayed in Figure 12a using violin plots layered with strip plots, considering SZA ≤ 75 • . The reader should note that the vertical axis scale used here is different from the one employed in Figure 10, extending up to 160%. The three model versions used aerosol data from reanalysis databases, and their output distributions are quite similar in shape, but the REST2 presents a higher positive bias. Most data points display errors between −10 and 25%, although outliers can reach errors of up to 150%. Among the three, McClear presented the lower number and magnitude of outliers, and there is a clear simililarity between the populations of outliers obtained with BRASIL-SR regional experiment and REST2, suggesting that MACC aerosol data is more consistent than MERRA-2 for ARM_Manacapuru. For Manaus_EMBRAPA, the results from the three models using reanalysis aerosol data are similar (Figure 12b), with higher positive median values for the DNI error for the McClear and REST2 outputs. However, the mean deviation values are partially compensated by the negative outliers that are higher in magnitude for these two models. The number and magnitude of both positive and negative outliers are much lower than for ARM_Manacapuru for all models.
In the violin plots for Brasilia_SONDA (Figure 12c) can be seen a similarity between the distributions for BRASIL-SR regional and REST2, with a broadening of the distribution towards positive DNI errors compared to the more symmetrical distribution presented by McClear output deviations. The MBD and median were closer to zero (Table 5) for BRASIL-SR regional experiment. A distinct group of seven negative outliers with errors of −25% occurred in BRASIL-SR regional experiment and REST2, suggesting they are linked to aerosol data provided by MERRA-2 used by both models. McClear displayed the highest median deviation values, in agreement with the performance metrics presented in Table 5.
Finally, the performance of McClear was distinctly better than the other two models in Palmas_SONDA (Figure 13), suggesting that MACC aerosol data is more consistent for Palmas_SONDA than MERRA-2, although this assessment is difficult to verify without AERONET data. For both BRASIL-SR regional and REST2, there is a significant number of positive outliers with large positive errors. It should be noted that the vertical axis scale is different for Palmas_SONDA violin plots than the scale used in Figure 12, extending up to 260%. Negative errors are also common in the three models, being more abundant for McClear at this particular site.

Discussion, Conclusions and Future Research
Overall, our results indicate a good skill of BRASIL-SR for the estimation of both GHI and DNI (the later for the case without optical depth scaling). For the period of study, the model skill showed a comparable or superior skill than the obtained with broadband models McClear and REST2 at the measurements sites, except for Palmas_SONDA where McClear presented the best skill. The great performance of McClear over Palmas_SONDA could be related to an improved representation of aerosol loading in the model for that portion of the domain. However, such a hypothesis requires more investigation once observational data of aerosol optical properties are scarce in the Palmas region, therefore, an issue for validation of aerosol products from both models and satellite algorithms.
The RMSD deviations of the GHI provided by BRASIL-SR using local AOD data (insitu experiment) was below 2.9% for all ground sites, that is, in the range of the instruments total uncertainty, between 3 to 5% [25]. Also for GHI, the MBD and MAD obtained for GHI were below 1.3% and 2.1% for the in-situ experiment, and below 2.1% and 4.4% for the regional experiment, and the RMSD was below 5.7% for the regional experiment.
Our results confirmed former research outputs and demonstrated that δ-Eddington scaling in the two-stream approximation led to an overestimation of the DNI. Mathematically, this behavior was expected because the scaled optical depth is always smaller than the actual optical depth. It was also observed that the higher AOD, the larger positive deviations of the DNI estimates based on the δ-Eddington scaling. On the other hand, our findings suggest that although DNI results without scaling seem more accurate than when using δ-Eddington scaling, there is still a negative bias in results without scaling, that increases in magnitude for larger AOD values. Such negative bias could be related to the contribution of the circumsolar irradiation not being evaluated by the model BRASIL-SR and the use of climatological aerosol intensive optical properties within the model. Joseph et al. [76] suggested that the underestimated circumsolar irradiance could compensate for the DNI overestimation. Our findings, however, suggest that the DNI overestimation due to scaling is typically larger in magnitude than the circumsolar underestimation, in consonance with Sun et al. [78]. Also supporting the use of an optical depth without scaling, Räisänen and Lindfors [79] showed that this is the best approach for an atmosphere with the presence of aerosols, whereas the δ-scaling could still be appropriate for cloudy atmospheres.
For the in-situ experiment, the MBD for the DNI without scaling varied from −2.3 to −0.5%. MAD varied between 1.6 and 3.6%, RMSD ranged from 2.3 to 4.7% and OVER between 0 and 5.3%. For the regional experiment, bias were slightly positive, with MBD ranging from 0.1 to 2.1%; MAD varied from 4.3 to 7.3%, RMSD ranged from 5.9 to 9.6% and OVER was near 1.7% for ARM_Manacapuru and Brasilia_SONDA, increasing to 6.2% for Manaus_EMBRAPA and reaching 26.9% for Palmas_SONDA. RMSD for DNI estimations with the in-situ experiment were comparable to the instrument uncertainty on field conditions for Brasilia_SONDA and Manaus_EMBRAPA (for pyrheliometers, the uncertainty on field is estimated to vary between 2% to 2.5%, while for RSI is around 5% [25]). For ARM_Manacapuru, the obtained RMSD was about twice the observations' uncertainty, but still about half the RMSD values obtained with the regional experiment, REST2 and McClear.
The deviations of DNI estimates provided by BRASIL-SR in-situ experiment using an optical depth without scaling are comparable in magnitude with those obtained by Ruiz-Arias et al. [22] in their study Control Experiment using SMARTS [80] and AERONET data, where the authors obtained negative MBD with values from −7.1% to 0.2% and RMSD ranged from 1.5% to 7.7%, depending on the location. On the other hand, when using MERRA-2 data, Ruiz-Arias et al. [22] obtained MBD values that ranged from −12.1% to −5.7%, and the RMSD ranged from 8.7% to 17.3%. While the comparison with these metrics seems to favor BRASIL-SR regional experiment results, the sites studied by Ruiz-Arias et al. [22] present, on average, higher AOD levels than the sites covered by our study and are dominated by coarse mode aerosol, both factors possibly contributing to larger errors in MERRA-2 aerosol data. The results of the regional experiment for DNI using an optical depth without scaling and MERRA-2 aerosol data were comparable or superior to those of McClear version 3.1 using MACC and REST2 version 5 using MERRA-2, for the same locations and times, again with the except for Palmas_SONDA where McClear presented the best performance with lower RMSD, MAD, KSI and OVER.
In summary, the aerosol emitted to the atmosphere by biomass burning events can get AOD values as high as 5.0 (although such extreme values were not observed in the sites and period covered by our study), and intensely attenuate the downward surface solar irradiance. The new improved version of the model BRASIL-SR provides GHI and DNI outputs presented low uncertainties for cloudless conditions for all aerosol loads considered. Since DNI and aerosol data are very scarce in the Amazon and Cerrado regions, future studies should cover more extensive timeframes and geographical areas, allowing a more comprehensive and detailed performance benchmark.
There is work in progress to improve cloud representation in BRASIL-SR. The major goal is to develop a reliable spectral model providing GHI, DNI for any atmospheric condition concerning cloud or aerosol optical thickness in a tropical region. The fine spatial resolution of the GOES imagery and its cloud and aerosol products can help to overcome the ground data scarcity. In addition, a parameterization for circumsolar irradiation such as the one proposed by Sun et al. [78] should be included in future versions of BRASIL-SR. Code optimization is also in progress to allow for more efficient memory use, facilitating running the BRASIL-SR for shorter time steps or larger spatial domains.  Data Availability Statement: BRASIL-SR model can be downloaded from http://labren.ccst.inpe. br/brasil-sr.html (last accessed on 10 October 2021), together with a set of example data for a day. The broadband and solar spectral data acquired at ARM_Manacapuru can be downloaded at the Atmospheric Radiation Measurement (ARM) website ( [51]). Brasilia_SONDA and Palmas_SONDA broadband solar data can be obtained from SONDA network [48]. MFRSR data from Manaus_EMBRAPA can be accessed upon request to Nilton E. Rosário (nrosario@unifesp.br). Finally, AERONET version 3.0 level 2.0 data can be downloaded from the AERONET site [81].

Acknowledgments:
We thank the PI investigators and their staff for establishing and maintaining the AERONET and ACONVEX(LFA-IFUSP) sites used in this investigation. We also thank the research team from the Laboratory of Modeling and Studies on Renewable Energy Resources for maintaining the SONDA network and providing quality surface data. We acknowledge the Federal University of São Paulo and the National Institute for Space Research (INPE) for supporting the research team. Finally, we thank Silvia Vitorino Pereira for the preparation of schematic of the modelling domain, presented in Figure 2 and the graphical abstract.

Conflicts of Interest:
The authors declare no conflict of interest. The funding institutions had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Biomass Burning Aerosol Optical Properties
Average optical properties for Cerrado and Amazonia regions were obtained for the biomass burning season using AERONET data version 3.0, level 2.0, ( [33,[82][83][84]) acquired between 1993 and 2019. Background and biomass burning conditions are separated using a threshold for AOD 500nm equal to the mean plus standard deviation value obtained from data acquired outside of the biomass burning season, following Pérez-Ramírez et al. [85]. In this manner, for Amazonia, biomass burning conditions were selected considering AOD 500nm ≥ 0.24, while the selection threshold was AOD 500nm ≥ 0.14 for Cerrado.
Average aerosol refractive real and imaginary index obtained for the AERONET wavelengths considering the aforementioned threshold for biomass burning conditions are interpolated between AERONET wavelengths using the pchip interpolate algorithm as implemented in SciPy Python library [86]. Refractive real and imaginary indexes for wavelengths lower than 440 nm are set to the value for this wavelength. Similarly, values for higher wavelengths are set to the values obtained for 1040 nm. The refractive indexes thus interpolated are combined with the average aerosol size distribution using Mie code [87] to obtain single scattering albedo and asymmetry factor to each one of the 37 spectral intervals of BRASIL-SR. The resulting spectral single scattering albedo and asymmetry parameter for Amazonia and Cerrado can be observed in Figures A1 and A2, as well as the AERONET values (mean and standard deviation). It can be seen that interpolated values agree with the AERONET data and are within the uncertainty range.

Appendix B. Spectral Irradiances with BRASIL-SR
Spectral GHI, DNI and DHI are available as output of BRASIL-SR. In this Appendix, an example of such output is presented. Figure A3 presents the GHI, DNI and DHI for 19 September 2015, at ARM_Manacapuru. In this example, DHI data is obtained from the difference between GHI and direct horizontal irradiance obtained without optical depth scaling since the latter was shown to be more accurate (Section 4.2). A point at 18:30 UTC was selected for the spectral output comparison.  Figure A4 displays BRASIL-SR spectral output between 200 and 2400 nm and MFRSR irradiation data. In the comparison, 5 min averages of MFRSR spectral irradiances, centered at the model output time, were used. A rather conservative uncertainty of 1.5% was considered for the spectral irradiance in all channels ( [53]). Modeled spectral GHI is in almost perfect agreement with MFRSR observations, as can be seen in Figure A4, left panel. Spectral DNI displays a slight underestimation in the visible range of the spectrum, although considering the uncertainties of the model (resulting from a fixed set of aerosol optical properties, interpolation of input data and others, in addition to the uncertainties of the model methods) and the differences in resolution between model spectral grid and MFRSR, the correspondence can be considered adequate. The underestimation of DNI results in a similar overestimation of DHI since the spectral DHI is calculated from the difference between GHI and the unscaled DNI.