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

Abstract

Although the Weather Research and Forecasting model with solar extensions (WRF-Solar) is tailed for solar energy applications, its official 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 significantly eliminate the AOD underestimations and reasonably reproduce the AOD temporal distributions, improving 51.63% for biases and 61.29% for correlation coefficients. Compared with the original WRF-Solar version, coupled online with an advanced aerosol module minifies the bias value of global horizontal irradiance (GHI) up to 44.52%, and AOD data assimilation contributes to a further reduction of 17.43%.

1. 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 carbon-neutral scenario to mitigate climate change in China affordably. To achieve the carbon-neutral 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 satellite-based 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 evaluated [15]. For example, the Global Environmental Multiscale (GEM) model [18], the European Center for Medium-Range Weather Forecasts [14], the Advanced Multiscale Regional Prediction System (ARPS) model [19], several versions of the Weather Research and Forecasting (WRF) model [20], and so on. Among them, by downscaling to a finer physical resolution realistically, mesoscale models have significant advantages in local forecast accuracy compared with global models [4,15,21,22].

The Weather Research and Forecasting model with solar extensions (WRF-Solar) is the first NWP model explicitly developed for solar power applications [23]. After improving the representation of the cloud–aerosol–radiation system [16,24] and incorporating a fast-radiative transfer algorithm [25], it is proved that GHI is improved by 46% on clear-sky assessment compared with the standard WRF at all the surface radiation budget network sites.

However, results based on NWP such as WRF-Solar still have significant uncertainties [26], which are caused by random and systematic errors and uncertainties when describing atmospheric conditions [6,27,28]. In addition, the prediction errors may come from a bad estimate of the initial and boundary conditions [29,30]. Boundary conditions can determine simulated aerosol effects on convective cloud fields [31] and the sensitivity of boundary layer dynamics in the model is highly related to initial meteorological conditions [32]. In addition, the official version of WRF-Solar uses a prescribed climatology-based aerosol parameterization in the radiation-related calculation [33,34], which is based on a model dataset [35]. Considering the rapid change in emissions, its 1975-average level may be far away from reality. Under clear skies, the aerosol-radiation interaction is the largest source of uncertainty in solar irradiance [23], especially in highly polluted areas [36] such as East and South Asia. Consequently, improving the AOD forecast under clear skies is essential for an effective and accurate solar prediction [37,38].

Aerosol data assimilation (DA) is a statistically optimal method to reduce the uncertainty in the initial aerosol field [39,40]. With the development of large observing systems and remote sensing techniques, abundant data are becoming available these years. This provides promising prospects for improving model initial and boundary conditions through DA [29]. Among all the DA methods, three-dimensional variational (3DVAR) is one of the most prevailing approaches applied in AOD assimilation [41,42]. Meanwhile, remote sensing optical properties show a significant advantage in covering large space spans and vertical detailed aerosol information [43,44,45]. AOD data from Moderate Resolution Imaging Spectroradiometer [43,45,46,47], Geostationary Ocean Color Imager and Visible Infrared Imaging Radiometer Suite [41], Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) satellite [48] and Himawari-8 [42,49] have been assimilated to WRF-Chem to improve AOD analyses and forecasts. At the same time, only radiances from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) aboard the Meteosat Second Generation (MSG) geostationary satellite have been assimilated into WRF Solar by 3DVAR to reduce the mean bias error [2]. There is still lacking research focusing on the effects of AOD DA in Asia using geostationary instruments such as Himawari-8 in WRF-Solar.

In this study, the chemistry aerosol module is firstly coupled online with the radiation module in WRF-Solar (i.e., WRF-Chem-Solar) to calculate the surface solar radiation (SSR) with time changing aerosols. Then we implement a 3DVAR aerosol DA system with the WRF-Chem-Solar model using the Himawari-8 geostationary satellite AOD data. Our major objectives are to quantitate the sensitivity of SSR to the aerosol initial conditions (ICs) and then investigate its forecasting performance under clear skies in China. The paper is organized as follows: Section 2 describes the WRF-Solar and its configurations as well as the data assimilation system. Section 3 reveals the satellite and observation datasets and the experiment designs. Section 4 provides the results and the statistical analysis. Section 5 draws conclusions.

2. Methodology

2.1. Forecast Model

In this study, the latest version of Weather Research and Forecasting model coupled with Chemistry (WRF-Chem v4.4) is applied to predict the SSR [50,51]. The model blends the advantages of the Multi-Sensor Advection Diffusion Nowcast (MADCast) [52,53] and WRF-Solar (i.e., MAD-WRF) [23] to achieve improved end-to-end solar irradiance forecasting. WRF-Solar is a specific configuration and augmentation of the WRF model. Compared with the original WRF version, the WRF-Solar focuses on improving the representation of the cloud–aerosol–radiation system by providing new parameterizations. It also outputs the direct and diffuse components of the irradiance at high temporal resolution while keeping the computational cost low [25]. Then, the feedbacks between sub-grid scale clouds and shortwave irradiance are implemented in a shallow cumulus parameterization [54]. All mentioned above contribute to the solar augmentations for solar forecasting in the WRF model.

In the frame of WRF-Solar, the effects of aerosol on the radiation simulation have been added only to the Goddard [55] and the Rapid Radiative Transfer Model for Global Circulation Models (RRTMG) [56] shortwave radiation schemes. However, the Goddard parameterization cannot output the clear-sky radiation at present; therefore, RRTMG is selected to simulate the SSR considering the aerosol-induced reduction with the predicted aerosols. The Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) is adopted to predict the major tropospheric aerosol components and their direct effects on SSR (i.e., WRF-Chem-Solar). In GOCART aerosol model, fifteen prognostic aerosol variables are used to predict the mass mixing ratios of the bulk sulfate, hydrophobic and hydrophilic organic carbon and black carbon (OC₁, OC₂, BC₁, and BC₂), other fine particulate matter (PM_2.5), dust in five bins and sea salt in four bins. An improved bulk microphysics [57] combined with a scale-aware cumulus parameterization [58] is employed to explicitly predict the cloud water, rain, cloud ice, snow, and graupel. The major configurations are shown in

Table 1. Configurations of the model used in the study.

Model Setting	Description	Reference
Microphysics	Thompson	[57]
Radiation	RRTMG scheme for SW and LW	[10,11]
Land Surface	Noah Land Surface Model	[12]
Cumulus Parameterization	Grell–Freitas ensemble scheme	[15]
Aerosol module	GOCART	[60]
Dust emission	GOCART dust emissions	[21]

, and the other configurations are as same as the earlier work [59]. The model is configured with a domain covering Eastern Asia with a 25 km horizontal resolution and 45 vertical levels. The top of the model is located at 50 hPa. The domain is shown in Figure 1a. In the current study, only the anthropogenic sources of aerosols are considered by using the HTAPv3 mosaic emission inventory for the year 2018 [60]. Considering the need for operational forecasts of air quality and SSR, the initial and lateral boundary meteorological fields are generated from the NCEP Global Forecast System (GFS).

2.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:

$$J\left(x\right)=\frac{1}{2}{\left(x-x^b\right)}^T{B}^{-1}\left(x-x^b\right)+\frac{1}{2}{\left[H\left(x\right)-y\right]}^T{R}^{-1}\left[H\left(x\right)-y\right],$$

(1)

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

$$J\left(\delta x\right)=\frac{1}{2}\delta x^T{B}^{-1}\delta x+\frac{1}{2}{\left[\mathrm{H}\delta x-d\right]}^T{R}^{-1}\left[\mathrm{H}\delta x-d\right],$$

(2)

$$\delta x=x-x^b,$$

(3)

$$d=y-H\left(x^b\right),$$

(4)

where $H$ is the nonlinear observation operator and $\mathrm{H}$ is the linearity of the observation operator in the vicinity of $x^b$. By defining $z=B^{-1}\delta x$ and assuming $\delta x=z=0$, the minimum of the cost function is found iteratively with a preconditioned conjugate gradient algorithm to solve both $\delta 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].

3. Observations and Experimental Design

3.1. 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 $\mathrm{H}$.

3.2. 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 dry-diameter 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.

3.3. 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 sun–sky 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.

3.4. 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.

3.5. Experiments

As summarized in

Table 2. Experimental design for the sensitivity tests in this study.

	Version	Chemistry	Online Couple	Data Assimilation
OR	WRF-Solar	-	-	-
FR	WRF-Chem-Solar	√	√	-
DA	WRF-Chem-Solar-DA	√	√	√

, 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:

$$\mathrm{IOA}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{O}}_{\mathrm{i}}{-\mathrm{M}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left(\left|{\mathrm{O}}_{\mathrm{i}}- \overline{\mathrm{O}}\right|+\left|{\mathrm{M}}_{\mathrm{i}}-\overline{\mathrm{O}}\right|\right)}^2}$$

(5)

where the ${\mathrm{O}}_{\mathrm{i}}$ means the observations, the ${\mathrm{M}}_{\mathrm{i}}$ means the model results, and the $\overline{\mathrm{O}}$ gives the average values of observations.

4. Results

4.1. 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 shows the horizontal distributions of the simulated time-averaged AODs in the FR and DA experiments and their differences (DA experiment minus FR experiment). For the FR experiment, high 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 Figure 2b,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² 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₂ due to China’s clean air action [82], with regional values up to 16 mg/m², 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.

4.2. 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].

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.

The simulated hourly AODs throughout 8–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.

Compared with the observations, results in the FR experiment underestimates the AODs, especially in relatively high-value stations. Their large RMSE values of the FR experiments are significantly removed by the assimilation of the Himawari-8 AODs. For example, observed AODs in Sanya show obvious daily variations, especially on 16 March 2021. 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). For example, the BIAS values are −0.430 in FR and −0.208 in DA; the latter reduces by about 51.62% compared to the former. The RMSE values show similar results. Meanwhile, the CORR and IOA values increase by about 61.29% and 40.61% compared to the FR experiment, indicating that the forecast fields are significantly improved by DA. 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.

4.3. 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² and RMSE beyond 159 W/m². 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², which is related to the better performance of AODs DA with departures greater than 0.5 in Figure 7b. All in all, 38.61% of stations have significant improvements (both for RMSE and CORR values), 59.51% of stations have lower RMSE values, and 53.46% of stations have higher CORR values in rural and suburban areas.

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. The comparison of the hourly observational GHI with simulated ones in OR, FR, and DA experiments in seven electricity regions, only rural and suburban stations are analyzed.

Regions	N	Experiments	BIAS	RMSE	CORR	IOA
North Eastern	1084	OR	49.980	93.664	0.920	0.943
		FR	30.918	84.651	0.921	0.953
		DA	23.240	83.873	0.917	0.954
Northern	731	OR	101.039	143.234	0.865	0.871
		FR	63.619	114.327	0.881	0.913
		DA	46.346	104.749	0.886	0.927
Eastern	342	OR	88.366	120.377	0.928	0.923
		FR	47.814	94.106	0.929	0.951
		DA	27.412	90.098	0.921	0.955
Southern	562	OR	111.054	176.936	0.840	0.887
		FR	71.735	151.104	0.846	0.915
		DA	41.776	133.711	0.854	0.932
Central	951	OR	72.450	208.284	0.737	0.841
		FR	19.840	185.317	0.766	0.870
		DA	17.020	184.346	0.768	0.871
Western	1328	OR	41.392	156.566	0.842	0.910
		FR	14.320	149.571	0.846	0.916
		DA	12.102	151.194	0.842	0.914
Tibetan	273	OR	146.070	230.096	0.818	0.851
		FR	137.654	226.330	0.815	0.854
		DA	135.003	226.468	0.810	0.854

N shows number of data used.

. 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 all seven regions. Moreover, the DA experiment shows the best performance among the three experiments for the BIAS values, with an average relative improvement beyond 53%. Meanwhile, East, Central, and West China show the most remarkable decreasing BIAS, with relative values beyond 68.98%, 76.50%, and 70.76% compared with the OR experiment. Except for West and Tibetan of China, the RMSE values in other regions are the lowest in the DA experiment. The inconspicuous effect of DA in the two regions may come from the absence of observations of Himawari-8 to be assimilated, as we talked about in Section 3.1.

Then seven typical rural and suburban stations from the main seven electricity regions are picked out, their information is given in

Table 4. The information of the 7 typical stations.

Station Name	Regions	Latitude	Longitude	Station Tsype
Erlianhaote	Northern	43.63°N	111.94°E	Prairie
Jiamusi	North Eastern	40.08°N	113.41°E	Plain
Panzhihua	Central	26.58°N	101.72°E	Mountainside
Cuona	Tibetan	27.98°N	91.95°E	Plateau
Huaian	Eastern	33.64°N	118.93°E	Suburban
Tengchong	Southern	24.98°N	98.51°E	Suburban
Xining	Western	36.66°N	101.73°E	Suburban

, 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.

5. 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.