Improving Clear-Sky Solar Power Prediction over China by Assimilating Himawari-8 Aerosol Optical Depth with WRF-Chem-Solar

: Although the Weather Research and Forecasting model with solar extensions (WRF-Solar) is tailed for solar energy applications, its ofﬁcial version lacks the consideration of the online aerosol-radiation process. To overcome this limitation, we have coupled the aerosol module online with the radiation module, then assimilated the high-resolution aerosol optical depth (AOD) from the Himawari-8 next-generation geostationary satellite using a three-dimensional variational (3DVAR) AOD data assimilation system to optimize the irradiance predictions with the better aerosol–radiation interaction. The results show that data assimilation can signiﬁcantly eliminate the AOD underestimations and reasonably reproduce the AOD temporal distributions, improving 51.63% for biases and 61.29% for correlation coefﬁcients. Compared with the original WRF-Solar version, coupled online with an advanced aerosol module miniﬁes the bias value of global horizontal irradiance (GHI) up to 44.52%, and AOD data assimilation contributes to a further reduction of 17.43%.


Introduction
Solar power is being used increasingly and exponentially as a kind of new and clean energy source and is expected to have a bright application prospect under the carbonneutral scenario to mitigate climate change in China affordably. To achieve the carbonneutral goal, photovoltaic power generation in China is expected to increase by 40% from the current 3.1% [1]. As a result, accurate forecasting of clear-sky solar power is becoming increasingly crucial to estimate the amount of energy available and the share of renewable sources in advance [2,3], particularly for clear-sky solar power prediction, during which sharing the most abundant solar energy.
Global horizontal irradiance (GHI) is the first and one of the most essential variables in most solar power prediction systems [4]. There are two typical approaches for forecasting solar irradiance [5,6]. One is statistical ways by using historical ground-based observations and satellite observations data with traditional mathematical models or updated machine learning techniques [7][8][9][10]. It saves calculation time but can only provide relatively good results for a temporal range between 30 min and 6 h [11][12][13]. The other is the numerical weather prediction (NWP) method, which usually shows better performances from 4 to 6 h onward by considering the dynamic phenomenon and cloud microphysics [14][15][16]. Some recent studies have shown that NWP forecasts are more competitive than satellitebased approaches to solar power forecasting, even on short time scales [4,17]. Typically, forecasts based on global, multiscale, and mesoscale NWP models are widely used and

Three-Dimensional Variational Data Assimilation
For computational affordability in the operational application, the 3DVAR technique in Grid-point Statistical Interpolation (GSI v3.7) is used to improve the aerosol forecasting with the optimal aerosol initial conditions (ICs) by assimilating satellite observed AODs in this study [61,62].
Fundamentally, data assimilation finds the best estimate of the analysis as the ICs for the subsequent forecast by minimizing the following scalar cost function:

.2. Three-Dimensional Variational Data Assimilation
For computational affordability in the operational application, the 3DVAR technique in Grid-point Statistical Interpolation (GSI v3.7) is used to improve the aerosol forecasting with the optimal aerosol initial conditions (ICs) by assimilating satellite observed AODs in this study [61,62].
Fundamentally, data assimilation finds the best estimate of the analysis as the ICs for the subsequent forecast by minimizing the following scalar cost function: where x and x b are vectors of analysis and forecast, respectively, which are the abovementioned prognostic aerosol variables over all model grids for this study; B and R are the background and observation error covariance matrices, which determine the relative contributions of the model forecast and observation terms to the final analysis; the nonlinear observation operator H maps the predicted aerosol mass mixing ratios to the satellite observed AODs for this study. The GSI 3DVAR adopts an incremental implementation of the cost function J(x) as following formulas where H is the nonlinear observation operator and H is the linearity of the observation operator in the vicinity of x b . By defining z = B −1 δx and assuming δx = z = 0, the minimum of the cost function is found iteratively with a preconditioned conjugate gradient algorithm to solve both δx and z [61]. In practical application, the recursive filters in GSI avoid inverting the huge matrix B and require only the standard deviation and horizontal and vertical length scales of the background error of the analysis variables [62]. In this study, the model named GENerate the Background Errors (GEN_BE) is applied to simulate the B with the differences of 12 and 24 h forecasts valid at the same time [63].

Himawari-8
Himawari-8 (H8) is a next-generation geostationary satellite located at 140.7 • E over the equator, which takes an Advanced Himawari Imager (AHI) to observe aerosol and clouds over the Asia-Oceania region from 80 • E to 160 • W between 60 • S to 60 • N [64]. The AHI has six bands in visible to near-infrared wavelength ranges, providing an operational full-disk near-real time aerosol optical properties every 10 min [65,66]. The aerosol optical properties in the latest version 031 are retrieved with a newly developed algorithm with the assimilated aerosol forecast as an a priori estimate [66], leading to more accurate aerosol retrievals. In addition, a scheme based on spatiotemporal aerosol variability is developed to improve merged hourly estimates of AOD with strict cloud screenings and provide the maximum number of available AOD retrievals [67]. Its products have been evaluated with observed station results by earlier studies and show good performance [68][69][70]. Taking advantage of the H8, we assimilate the AHI level 3 hourly merged AODs at 500 nm to improve the aerosol and subsequent SSR forecast by correcting the initial aerosol conditions in this study.
The officially public GSI 3DVAR aerosol assimilation system has the capacity to assimilate the Moderate Resolution Imaging Spectroradiometer (MODIS) AODs at 550 nm [61]. We further extend the assimilation system to include the H8 hourly merged AODs at 500 nm, and the efficient aerosol optical module specifically developed for the GOCART based on the look-up table in the Community Radiative Transfer Model (CRTM) is severed as the nonlinear observation operator H and its tangent linear operator H.

AERONET
The Aerosol Robotic Network (AERONET) is a globally distributed ground-based network of automatic tracking sun and sky scanning radiometers [71]. The directly spectral solar extinction measurements are used to retrieve the AODs at wavelengths of 340, 380, 440, 500, 675, 870, 1020, and 1640 nm with high accuracy. The estimated uncertainties of the retrieved AODs are about 0.01-0.021, with slightly higher errors in the ultraviolet wavelengths. The cloud-screened and quality-controlled AERONET AODs with the latest version 3 retrieval algorithm are used to independently validate the improvement of the aerosol forecast by assimilating the H8 AODs [72].
As WRF-Chem only outputs the AODs at 300, 400, 550, and 999 nm calculated with the aerosol optical module based on Mie scattering theory, the AERONET AODs at 440 and 675 nm are logarithmically interpolated to compare with the modeled ones at 550 nm. It is notable that the official version of WRF-Chem calculates the aerosol optical properties by assuming bi-modally lognormal size distributions of GOCART dry bulk sulfate and carbonaceous aerosols, and the mean diameters (geometric standard deviations) of the aitken and accumulation modes are 0.01 (1.7) and 0.07 (2.0), respectively. To be consistent with the size distributions assumed in the CRTM aerosol optical module, we modified the size distributions of sulfate and carbonaceous aerosols as mono-modally lognormal distributions. The mean diameters of the sulfate, BC, and OC are assigned as 0.139, 0.0236, and 0.0424 [73]. The simulated AODs are calculated with the WRF-Chem "aerosol chemistry to aerosol optical properties" module, which are the integral results of 8 sizes with drydiameter ranges from 0.039 to 10 um. The details of the particle dry-diameter range for 8 size bins can be found in an earlier study [74]. Totally, 94 stations are used in this study.

SONET
The Sun-Sky Radiometer Observation Network (SONET) is a ground-based aerosol observation network providing columnar atmospheric aerosol properties in the long run in China [75]. A total of 20 long-term observation stations are providing key aerosol parameters, including AOD, asymmetry factor, single scattering albedo, and Ångström exponent of aerosol particles within the entire atmospheric column. The multi-wavelength polarized sunsky radiometers CE318-DP are used to measure aerosol and water vapor at nine channels with center wavelengths of 340, 380, 440, 0.50, 670, 870, 1020, and 1640 nm at every 15 min [75,76]. The average differences in AOD between the SONET and AERONET are only 0.002 [75]. Similar to AERONET AODs, the SONET AODs at 440 and 675 nm are logarithmically interpolated to 550 nm. After quality control, 18 stations are used in this study.

SSR Stations
The hourly ground-based SSR data across China were available from the China Meteorological Administration (CMA), which contains data under the all-sky and clear-sky scenarios. After strict data quality control, such as the spike value test and homogenization [77,78], 138 sites are used to verify the performance of the simulated clear-sky SSR.

Experiments
As summarized in Table 2, three numerical experiments are conducted in this study. A base experiment (OR) with the default WRF-Solar model (which has not been coupled online with chemistry) configuration is run from 8 to 22 March 2021 as a reference to evaluate the performance of the original WRF-Solar version in China and surrounding areas. A control experiment (FR) with a full chemistry aerosol model GOCART coupled with WRF-Solar is conducted to evaluate the effects of the time-varying aerosols on SSR simulation. Similar to the FR experiment, an aerosol assimilation experiment (DA) is performed to correct the aerosol initial conditions every 24 h for evaluation of the effects of the Himawari-8 AOD assimilation on the hourly SSR simulation. All the experiments are conducted with a cold start at 06:00 UTC every day to ingest the meteorological fields, whereas the aerosol initial conditions are from the 24 h forecasts on the previous day in the FR experiment and the optimal ones by assimilating the H8 AODs in the DA experiment. The FR and DA experiments start from a clean atmosphere with a spin-up period of a week as the aerosol lifetime is generally less than one week.

√ √ √
Statistical metrics, including the mean error (BIAS), the root mean square error (RMSE), the correlation coefficient (CORR), and the index of agreement (IOA), are calculated to reveal the model performance [79]. IOA is a dimensionless statistical measure of model performance as following: where the O i means the observations, the M i means the model results, and the O gives the average values of observations.

DA Impacts on Aerosol Spatial Column Distributions
Due to the lack of aerosol simulation in OR experiment, only aerosol spatial column distributions in FR and DA experiments are analyzed here. Figure 2  becomes smaller and moves northwest, while a low-to-high shift (from about 0.70 to nearly 1.1) can be found in areas in South and Southeast Asia (Figure 2b). Similar results have been detected by earlier studies [80,81], which reveal that the observed AODs higher than 0.7 are consistently examined over the heavily industrialized and densely populated regions of India, Bangladesh, and China in spring. The relatively non-obvious change (with an absolute value below 0.04) may result from the absence of observations of Himawari-8 to be assimilated over the west of 80 • E, leading to a lack of assimilation data on the periphery of the measurement range.
AOD values are found in Xinjiang province, central China, and South Asia, where human activities are frequent. At the same time, low AOD values below 0.05 are found over the Tibetan Plateau and neighboring areas where anthropogenic aerosol emissions are limited. As shown in Figures 2b and c, after DA, it is apparent that DA enlarges the AOD values over the human-active areas except for Northwest China, part of Southwest China, and part of Central and North Asia. There are significantly increasing AODs in South Asia, Southeast Asia, the North China Plain (NCP), and the southwest and southeast coast of China on land, with average values passing 0.32. In the DA experiment, the range of low-value centers becomes smaller and moves northwest, while a low-to-high shift (from about 0.70 to nearly 1.1) can be found in areas in South and Southeast Asia (Figure 2b). Similar results have been detected by earlier studies [80,81], which reveal that the observed AODs higher than 0.7 are consistently examined over the heavily industrialized and densely populated regions of India, Bangladesh, and China in spring. The relatively non-obvious change (with an absolute value below 0.04) may result from the absence of observations of Himawari-8 to be assimilated over the west of 80°E, leading to a lack of assimilation data on the periphery of the measurement range. Then, the influences of DA on the spatial distributions of aerosol component burdens during a severe spring Eastern Asian dust event is analyzed. Figure 3 shows the spatial distributions of the daily mean values of the main individual aerosol component burdens (containing dust, sulfate aerosol, organic carbon (OC), and black carbon (BC)) on March 16th in the FR and DA experiments, as well as their differences on that day. With obvious increasing dust from Northwest China, DA has contributed to an increasing dust value beyond 400 mg/m 2 in NCP and Northwest China (Figure 3a-c). The noticeable increasing dust transport from South and Southeast Asia to South China can also be detected in DA Then, the influences of DA on the spatial distributions of aerosol component burdens during a severe spring Eastern Asian dust event is analyzed. Figure 3 shows the spatial distributions of the daily mean values of the main individual aerosol component burdens (containing dust, sulfate aerosol, organic carbon (OC), and black carbon (BC)) on March 16th in the FR and DA experiments, as well as their differences on that day. With obvious increasing dust from Northwest China, DA has contributed to an increasing dust value beyond 400 mg/m 2 in NCP and Northwest China (Figure 3a-c). The noticeable increasing dust transport from South and Southeast Asia to South China can also be detected in DA experiment. Concurrently, after DA, there are increased sulfate aerosol, OC, and BC over most parts of South and Southeast Asia and the adjacent seas, and reduced ones in most areas in East China. It should be mentioned that the reductions over East China in the DA experiment are dominated by the sulfate aerosol because of the significant reductions of the SO 2 due to China's clean air action [82], with regional values up to 16 mg/m 2 , which is consistent with earlier work [49]. experiment. Concurrently, after DA, there are increased sulfate aerosol, OC, and BC over most parts of South and Southeast Asia and the adjacent seas, and reduced ones in most areas in East China. It should be mentioned that the reductions over East China in the DA experiment are dominated by the sulfate aerosol because of the significant reductions of the SO2 due to China's clean air action [82], with regional values up to 16 mg/m 2 , which is consistent with earlier work [49]. The above results demonstrate that the assimilation of Himawari-8 AODs can significantly alter the spatial distribution of AODs and the aerosol column burdens of relevant individual components, thus changing the spatial pattern of SSR.

Verification of DA on Aerosol ICs and Forecasts
Then, the improvements of aerosol initial conditions and aerosol forecast by assimilating Himawari-8 AODs are verified. We conducted a comparative analysis of the Himawari-8 observations, simulated AODs in the FR experiment (FR AODs), the first guess, and analyses (FG and DA AODs) in the DA experiment, the results are shown in Figure 4. The first guess AODs are the forecasts for another 24 h from a prior DA cycle a day ago, and the DA AODs are the results of the assimilation of Himawari-8 data. It should be noted that the comparisons are conducted between the observed and simulated AODs at 06:00 UTC every day. Such evaluations are considered internal checks rather The above results demonstrate that the assimilation of Himawari-8 AODs can significantly alter the spatial distribution of AODs and the aerosol column burdens of relevant individual components, thus changing the spatial pattern of SSR.

Verification of DA on Aerosol ICs and Forecasts
Then, the improvements of aerosol initial conditions and aerosol forecast by assimilating Himawari-8 AODs are verified. We conducted a comparative analysis of the Himawari-8 observations, simulated AODs in the FR experiment (FR AODs), the first guess, and analyses (FG and DA AODs) in the DA experiment, the results are shown in Figure 4. The first guess AODs are the forecasts for another 24 h from a prior DA cycle a day ago, and the DA AODs are the results of the assimilation of Himawari-8 data. It should be noted that the comparisons are conducted between the observed and simulated AODs at 06:00 UTC every day. Such evaluations are considered internal checks rather than independent validations and can reveal the influence of aerosol assimilation on the model simulations [43].
Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 20 than independent validations and can reveal the influence of aerosol assimilation on the model simulations [43]. The spatial distributions of the mean AODs observed by the Himawari-8, simulated in the FR experiment and DA experiment, and their differences with Himawari-8 at daily 06:00 UTC are shown in Figure 4. Clearly, the spatial distributions of the FG and DA AODs are much more comparable to the Himawari-8 observations than those in the FR experiment, especially over South and Southeast Asia and North China. It proves that the Himawari-8 AODs are successfully assimilated into the model to correct the aerosol initial conditions, leaving much more accurate ICs. The FR AODs miss most of the high values (beyond 1.1) over South and Southeast Asia and the adjacent Indian Ocean, with mean biases up to −0.50 (Figure 4c). The DA approach reproduces a more reasonable spatial distribution and corrects the generally negative initial biases in South and Southeast Asia The spatial distributions of the mean AODs observed by the Himawari-8, simulated in the FR experiment and DA experiment, and their differences with Himawari-8 at daily 06:00 UTC are shown in Figure 4. Clearly, the spatial distributions of the FG and DA AODs are much more comparable to the Himawari-8 observations than those in the FR experiment, especially over South and Southeast Asia and North China. It proves that the Himawari-8 AODs are successfully assimilated into the model to correct the aerosol initial conditions, leaving much more accurate ICs. The FR AODs miss most of the high values (beyond 1.1) over South and Southeast Asia and the adjacent Indian Ocean, with mean biases up to −0.50 (Figure 4c). The DA approach reproduces a more reasonable spatial distribution and corrects the generally negative initial biases in South and Southeast Asia (Figure 4e), with lower mean biases below 0.30 (absolute value). The FG AODs further confirm that the aerosol forecasts benefit from correcting aerosol initial conditions by assimilation of the Himawari-8 AODs.
Figure 5a-c further show the comparisons of the probability density of the FR, FG, and DA AODs with the Himawari-8 observations at 06:00 UTC every day. Obviously, the FR AODs tend to be underestimated compared with the Himawari-8 observations (Figure 5a), which will lead to large biases when serving as the ICs for the forecasting of the next days. With assimilation, the FG AODs are more concentrated, and more points align with the 1:1 line (Figure 5b). The BIAS values are −0.340, −0.135, and −0.040 for the FR, FG, and DA AODs, respectively, showing a significant bias reduction by DA. With the DA approach, the RMSE value of FG AODs is lower (0.290) and decreased by 37.90% compared to the FR ones. Meanwhile, DA improves the AOD ICs with a CORR value over 0.937 for DA AODs (Figure 5c), followed by a CORR value of 0.726 for the FG AODs, the latter increases by 38.02% compared to the FR AODs. For the FR, FG, and DA AODs, their IOA values are showing an increasing trend, from 0.554 to 0.957. Figure 5d further shows frequency distributions of AOD deviations (PDF). Tending to underestimate the AODs over South and Southeast Asia (Figure 4), the distribution of the deviations for FR AODs is positively biased, with 26.11% of deviations dropping into the departures greater than 0.50. After successfully assimilating the AODs into the WRF-Chem-Solar model, the distribution of the deviations for DA AODs essentially reduces bias and obeys normal distribution. The distribution of FG AODs shows reduced bias, with 43.85% of deviations dropping into the absolute departures less than 0.10, indicating that the assimilation exhibits improved consistency with observations compared with the FR AODs.
Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 20 ( Figure 4e), with lower mean biases below 0.30 (absolute value). The FG AODs further confirm that the aerosol forecasts benefit from correcting aerosol initial conditions by assimilation of the Himawari-8 AODs.
Figures 5a-c further show the comparisons of the probability density of the FR, FG, and DA AODs with the Himawari-8 observations at 06:00 UTC every day. Obviously, the FR AODs tend to be underestimated compared with the Himawari-8 observations (Figure  5a), which will lead to large biases when serving as the ICs for the forecasting of the next days. With assimilation, the FG AODs are more concentrated, and more points align with the 1:1 line (Figure 5b). The BIAS values are −0.340, −0.135, and −0.040 for the FR, FG, and DA AODs, respectively, showing a significant bias reduction by DA. With the DA approach, the RMSE value of FG AODs is lower (0.290) and decreased by 37.90% compared to the FR ones. Meanwhile, DA improves the AOD ICs with a CORR value over 0.937 for DA AODs (Figure 5c), followed by a CORR value of 0.726 for the FG AODs, the latter increases by 38.02% compared to the FR AODs. For the FR, FG, and DA AODs, their IOA values are showing an increasing trend, from 0.554 to 0.957. Figure 5d further shows frequency distributions of AOD deviations (PDF). Tending to underestimate the AODs over South and Southeast Asia (Figure 4), the distribution of the deviations for FR AODs is positively biased, with 26.11% of deviations dropping into the departures greater than 0.50. After successfully assimilating the AODs into the WRF-Chem-Solar model, the distribution of the deviations for DA AODs essentially reduces bias and obeys normal distribution. The distribution of FG AODs shows reduced bias, with 43.85% of deviations dropping into the absolute departures less than 0.10, indicating that the assimilation exhibits improved consistency with observations compared with the FR AODs.  The simulated hourly AODs throughout [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] March 2021 are also evaluated with independent observations from AERONET and SONET. As shown in Figure S1, the high values of the observed mean AODs are located in stations near South (beyond 0.70) and Southeast Asia (beyond 1.10), NCP (beyond 1.10), and the coastal regions of Southeast China (beyond 0.70), which are in consistency with results in Figure 3b. Six typical sites are chosen to verify the AOD time series. Both stations in China and nearby regions are included to explore the DA effect, and only three SONET stations in China (Taiwan, Sanya and Jiaozuo stations) are chosen due to data missing. The comparisons of observed and simulated hourly AOD are given in Figure 6. The high AODs on that day is successfully reproduced by the DA experiment, and it is consistent with the DA results in Figure 3. DA helps to decrease its RMSE by 54.95%. At Hankuk University (Hankuk_UFS station), the sudden increase of AODs on 16 March is also detected by the DA experiment. In Doi Ang Khang station, the observed AODs remain at relatively high values due to serious anthropogenic emissions in Southeast Asia on March, which is consistent with DA results in Figure 3. Its RMSE value reaches 1.235 in the FR experiment, with a poor CORR value of -0.150. A significant RMSE reduction and CORR elevation can be found by DA, with relative improvement up to 53.20% for RMSE. All in all, with higher CORR values, all the stations correctly modify the AODs variations to match the observations better. Figure 7a and b show the comparison of the probability density of the observed hourly AODs versus the ones from the simulated results with and without the assimilation. Compared to the independent observations, apparently, the simulated AODs in the FR experiment led to obvious underestimations (Figure 7a). AODs assimilation essentially reduced the tendency and is superior to those in the FR experiment as indicated by the BIAS, RMSE, CORR, and IOA values (Figure 7b) Figure 7c shows the frequency distributions of AOD deviations. In the FR experiment, it reveals a positively bias, whereas in the DA The high AODs on that day is successfully reproduced by the DA experiment, and it is consistent with the DA results in Figure 3. DA helps to decrease its RMSE by 54.95%. At Hankuk University (Hankuk_UFS station), the sudden increase of AODs on 16 March is also detected by the DA experiment. In Doi Ang Khang station, the observed AODs remain at relatively high values due to serious anthropogenic emissions in Southeast Asia on March, which is consistent with DA results in Figure 3. Its RMSE value reaches 1.235 in the FR experiment, with a poor CORR value of -0.150. A significant RMSE reduction and CORR elevation can be found by DA, with relative improvement up to 53.20% for RMSE. All in all, with higher CORR values, all the stations correctly modify the AODs variations to match the observations better. Figure 7a and b show the comparison of the probability density of the observed hourly AODs versus the ones from the simulated results with and without the assimilation. Compared to the independent observations, apparently, the simulated AODs in the FR experiment led to obvious underestimations (Figure 7a). AODs assimilation essentially reduced the tendency and is superior to those in the FR experiment as indicated by the BIAS, RMSE, CORR, and IOA values (Figure 7b) Figure 7c shows the frequency distributions of AOD deviations. In the FR experiment, it reveals a positively bias, whereas in the DA experiment, the results not only show reduced bias but show a peak closer to 0. Compared with FR ones, about 36.00% of AOD deviations drop into the absolute departures less than 0.10 for the DA experiment, and 43.48% of AOD deviations are reduced in data dropping into the departures greater than 0.50. We also find that the DA is more efficient in correcting highly overrated AODs with departures greater than 0.50 compared with underestimated values with departures less than −0.50. Similar results have been observed during earlier works on assimilation AODs [49,83]. The above findings indicate that AOD assimilation can significantly improve the model performance.
ote Sens. 2022, 14, x FOR PEER REVIEW 12 of 20 experiment, the results not only show reduced bias but show a peak closer to 0. Compared with FR ones, about 36.00% of AOD deviations drop into the absolute departures less than 0.10 for the DA experiment, and 43.48% of AOD deviations are reduced in data dropping into the departures greater than 0.50. We also find that the DA is more efficient in correcting highly overrated AODs with departures greater than 0.50 compared with underestimated values with departures less than −0.50. Similar results have been observed during earlier works on assimilation AODs [49,83]. The above findings indicate that AOD assimilation can significantly improve the model performance.

Improvement of the Clear-Sky Solar Power Prediction by DA
The variations of AODs will change the evaluations of clear-sky solar power prediction. To access the performance of clear-sky GHI simulation at an hourly scale in WRF-Chem-Solar with and without assimilation, the GHI simulated by the original WRF-Solar version is also represented. The GHI values simulated by all three experiments are performed in independent validations over the study period through a comparison with SSR stations in China. It should be mentioned that to avoid the influence of nearby sources, buildings, and heat islands, only rural and suburban stations are analyzed here, their locations are shown in Figure S2. Beyond that, we also evaluate simulation performance in the urban stations, which show consistent results (Table S1 and Figure S4).
As shown in Figure 8, the probability density and PDF plots of the observed hourly GHI versus the simulated ones under clear skies for the three experiments are presented.

Improvement of the Clear-Sky Solar Power Prediction by DA
The variations of AODs will change the evaluations of clear-sky solar power prediction. To access the performance of clear-sky GHI simulation at an hourly scale in WRF-Chem-Solar with and without assimilation, the GHI simulated by the original WRF-Solar version is also represented. The GHI values simulated by all three experiments are performed in independent validations over the study period through a comparison with SSR stations in China. It should be mentioned that to avoid the influence of nearby sources, buildings, and heat islands, only rural and suburban stations are analyzed here, their locations are shown in Figure S2. Beyond that, we also evaluate simulation performance in the urban stations, which show consistent results (Table S1 and Figure S4).
As shown in Figure 8, the probability density and PDF plots of the observed hourly GHI versus the simulated ones under clear skies for the three experiments are presented. Due to the lack of aerosols, GHI in the OR experiment is overestimated (Figure 8a), with BIAS beyond 72 W/m 2 and RMSE beyond 159 W/m 2 . Apart from the minor improvement of CORR and IOA values, the BIAS and RMSE valves in the FR and DA experiments are all superior to the OR ones and present best performances in DA experiment (Figure 8b,c). Their relative improvements (compared with the OR experiment) are 44.52% and 61.95% for BIAS, 9.75%, and 11.50% for RMSE, separately. Compared to the OR experiment, the PDF distributions for the other two experiments all show reduced bias, as indicated by higher frequencies of deviations within 20 (Figure 8d). Meanwhile, the weak negative biases of the OR experiment are improved by the experiments online coupled with chemistry, leading to peaks closer to 0. We also find that the DA is more efficient in correcting overly underestimated GHI with departures less than −50 W/m 2 , which is related to the better performance of AODs DA with departures greater than 0.5 in Figure 7b. All in all, 38 Due to the lack of aerosols, GHI in the OR experiment is overestimated (Figure 8a) Figure 8d). Meanwhile, the weak negative biases of the OR experiment are improved by the experiments online coupled with chemistry, leading to peaks closer to 0. We also find that the DA is more efficient in correcting overly underestimated GHI with departures less than −50 W/m 2 , which is related to the better performance of AODs DA with departures greater than 0.5 in Figure 7b. All in all, 38  We further explore the GHI simulation performance in the main seven electricity grids in China according to the grid distribution [84], including Northeastern China, Northern China, Eastern China, Southern China, Central China, Western China, and Tibetan Plateau (Tibetan), separately. The partition is shown in Figure 1b. The comparisons of observed and simulated mean GHI and their relative change at rural and suburban stations under clear skies in different regions are given in Table 3. Comparatively speaking, AOD assimilation plays a weak role in improving the CORR and IOA. It is clear that the major and most obvious improvements are represented in revisions of the large BIAS and RMSE values. The GHI values are obviously overestimated in the OR experiment for   Then seven typical rural and suburban stations from the main seven electricity regions are picked out, their information is given in Table 4, and their time series of the observed and simulated hourly GHI under clear skies are shown in Figure 9. Apparently, GHI in OR and FR experiments tend to overestimate the daily high values in most stations, especially in Huaian and Tengchong stations. After DA, the simulated GHI values are closer to the observed ones, and the overestimation trends are corrected. Huaian station shares the best corrections, the RMSE value decreases by 43.07% in the DA experiment compared with the OR one. Overall, DA help to reproduce better variations from hours to days.

Conclusions
Assimilating high-quality geostationary satellite aerosol optical depth (AOD) datasets is a useful approach to improve global horizontal irradiance (GHI) predictions. In this study, the chemistry aerosol model is firstly coupled online with the radiation model in WRF-Solar to calculate the surface solar radiation (SSR) with real-time changing aerosol. Based on the newly developed WRF-Chem-Solar model, a Three-Dimensional Variational (3DVAR) aerosol assimilation system is further adopted to assimilate the AOD observations from a new generation geostationary satellite (Himawari-8) every 24 h window from 8 to 22 March 2021. The performances of the simulated hourly output of AODs and SSR are evaluated by the assimilated Himawari-8 AOT observations, the independent AERONET/SONET AOD observations and SSR observations. The major findings from this work are summarized as follows:

Conclusions
Assimilating high-quality geostationary satellite aerosol optical depth (AOD) datasets is a useful approach to improve global horizontal irradiance (GHI) predictions. In this study, the chemistry aerosol model is firstly coupled online with the radiation model in WRF-Solar to calculate the surface solar radiation (SSR) with real-time changing aerosol. Based on the newly developed WRF-Chem-Solar model, a Three-Dimensional Variational (3DVAR) aerosol assimilation system is further adopted to assimilate the AOD observations from a new generation geostationary satellite (Himawari-8) every 24 h window from 8 to 22 March 2021. The performances of the simulated hourly output of AODs and SSR are evaluated by the assimilated Himawari-8 AOT observations, the independent AERONET/SONET AOD observations and SSR observations. The major findings from this work are summarized as follows: Firstly, the results show that the assimilation of Himawari-8 AODs can significantly alter the spatial distribution of AODs and the relevant individual aerosol components.
Secondly, the impact of AOD data assimilation (DA) on aerosol initial conditions (ICs) and forecasts are verified by Himawari-8 AOD data and independent observational sites. Compared with experiments without assimilation (FR), DA reproduces a more reasonable spatial distribution of AOD ICs and corrects the generally negative initial biases in South and Southeast Asia and the adjacent seas. For the first guess AODs, the CORR value increases by 38.02% compared with the FR ones. Generally speaking, DA can significantly improve model AOD ICs by mitigating the underestimations. After DA, the negative BIAS and RMSE of forecast in the FR experiment are greatly improved with 50% reductions, and the CORR and IOA biases are reduced by about61.29% and 40.61%, separately.
Finally, improvement of the clear-sky GHI by DA is quantified by independent observations. Overall, due to the lack of online coupling of aerosol processes, GHI in the original WRF-Solar version is overestimated. The FR and DA experiments present significant better performances with relative improvements of 44.52% and 61.95% for BIAS, 9.75% and 11.50% for RMSE, separately. We also find that the DA experiment is more efficient in correcting overly underestimated GHI than the overestimated ones.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/rs14194990/s1, Figure S1: Spatial distributions of the (a) observed mean AODs; Figure S2: The locations of the SSR stations and their types; Figure S3