A Study on the Influence of Air Pollution on Temperature Forecast Skill Based on Operational Weather Forecast in BTH Region

: Surface air temperature is a comprehensive function of aerosols in the atmosphere and various weather factors. However, there is no real-time aerosol concentration feedback in most operational numerical weather prediction (NWP) models. This raises a scientific question of how abnormal changes in air pollutants in a short period of time will affect the temperature prediction skill of NWP models. Thus, the study was carried out to investigate the possible influence of air pollution on the temperature forecast skill based on the operational NWP model over the Beijing–Tianjin–Hebei (BTH) region during January–February 2020. The results show that the average concentrations of PM 2.5 , SO 2 , NO 2 and CO over the BTH region in February were smaller than those in January by 38.5%, 35.1%, 48.0% and 33.1%, respectively. Simultaneously, the forecast skills for surface temperature in February from both regional (RMAPS, Rapid-refresh Multi-scale Analysis and Prediction System) and global (ECMWF, European Centre for Medium-Range Weather Forecasts) operational NWP models improved markedly compared with that in January. In both models, the underestimation of maximum temperature and the overestimation of minimum temperature in most cities over the BTH region in February were significantly reduced. With the 24 h (24 h) forecast lead time, the RMSE (root mean square error) of BTH daily mean, maximum and minimum temperature prediction in February based on RMAPS were 17.3%, 9.8% and 21.6% lower than that in January, respectively. These are generally consistent with the other statistical indices such as deviation and regression coefficient. As the forecast lead time extended to 48 h and 72 h forecast, the phenomena still existed and were also evident in the ECMWF model. The improvement of temperature forecast skill of NWP models may be attributed to the unexpected dramatical reduction of air pollutants. Less aerosols during the daytime allow more solar radiation reaching the surface and cause a warming in the near-surface temperature, while less aerosols during the nighttime favor the outgoing long-wave radiation and then lead to a cooling near the ground.


Introduction
The Beijing-Tianjin-Hebei (BTH) region is located in the north of the North China Plain with approximately 110 million residents and 216,000 km 2 in size. With the rapid urbanization and industrial development over the past three decades, the BTH region has become one of China's most economically developed regions and the third economic engine in China. However, with the heavy industrialization and rapid urbanization, as well as high-population density and special terrain and meteorological conditions, the BTH region has also become one of the most polluted hotspots in China [1][2][3]. Although the unprecedented and nationwide environmental governances have been executed since 2013 to reduce air pollutants emission and the air quality in most Chinese cities have been greatly improved, the BTH region is still one of the most polluted areas in China [4][5][6][7]. Air pollution or haze events in the BTH region, characterized by low visibility and high concentration of aerosol particles in aerodynamic diameter less than 2.5 µm (PM 2.5 ), greatly threaten human health, traffic safety and ecosystems [8][9][10][11].
Meteorological conditions play important roles in the formation and evolution of air pollution. At the same time, they are also affected by aerosol pollution through direct or indirect, physical or chemical effects. Therein, atmospheric aerosols are liquid or solid particles suspended in the air, which are multi-facet agents that play great roles in visibility, air quality, weather and climate and others through various mechanisms, such as atmospheric radiation, cloud physics, and hydrologic cycle processes [12][13][14][15][16][17]. An important and interesting question is whether atmospheric aerosol pollution would affect the forecast performance of the real-time operational numerical weather prediction (NWP) model systems, especially when the pollution status changes greatly and unexpectedly in a short period of time. To date, this is still an open scientific problem that needs more exploration.
Because of the high cost of running a fully integrated air pollution system at the high resolution of current NWP systems, the real-time coupling effects between atmosphere and chemistry have not been taken into account in the mainstream operational NWP models. They usually use the climatologically radiative transfer parameters based on the outdated aerosol datasets instead of the realistic aerosol-radiation interactions [18,19]. Some studies, mainly in North America and Europe, suggested that the numerical simulations with aerosol feedback mechanism can generally improve the prediction skills of temperature, relative humidity and other factors, especially near the large pollution sources [20][21][22]. Although some recent studies pointed out that, whether the feedback effect of real-time air pollution is considered in numerical weather prediction models, may have a nonnegligible impact on the forecast bias of temperature and precipitation, and many weather systems such as Madden-Julian oscillation (MJO), intertropical convergence zone (ITCZ) and monsoon systems [23][24][25][26], there are fewer studies for understanding the specific impact of air pollution on the regional or local weather forecast performance from the perspective of real-time and operational NWP models, especially for the seriously polluted BTH region.
In a short period of time, an operational NWP system should be generally stable in its forecasting skills. Theoretically speaking, without the real-time meteorology-chemistry coupling effects, the NWP model is unlikely to adjust promptly to their forecasting capabilities with the different air pollution levels in time. Due to the spread of coronavirus disease , the strict industrial production and transportation restrictions and quarantines were executed, approximately, since the end of January in most of China and then many other countries, which led to an obvious reduction in pollutant emissions and improvement of air quality [27,28]. The spaceborne NO 2 column observations revealed unprecedented decreases of NO 2 over China, South Korea, Western Europe and the U.S. as a result of public health measures enforced to contain the coronavirus disease outbreak in January-April 2020 [29]. The concentrations of other pollutants such PM 2.5 , SO 2 and CO also decreased significantly during this period in China, especially in BTH, which will be presented below based on the hourly station observations. It is scientifically important to understand whether the unexpected significant decrease in pollutants (especially the aerosols, i.e., PM 2.5 ) would cause a discernable impact on weather forecasts over the BTH region during this period. A systematic understanding of the forecast bias driven by pollution is beneficial to improve the temperature prediction skill, and even further to serve the accurate prediction of pollutant concentration based on diverse forecasting methods, such as machine learning or AI (artificial intelligence) models [30][31][32]. Therefore, the purpose of this study is to investigate the possible influence of unexpected significant changes in pollutant concentration in a short period of time on the temperature prediction skill in BTH based on the operational NWP models. The possible associations between the meteorological condition, air pollutant and temperature forecast are also preliminary explored in the study.
The rest of this paper is organized as follows. Section 2 describes the data and methods used and the case selection. The pollutant concentration changes and the comparison of temperature prediction performance based on different NWP models are presented in Section 3. Section 4 discusses the possible physical links involved in the meteorological Atmosphere 2023, 14, 1229 3 of 18 condition and temperature prediction skills. Finally, the main conclusions are summarized in Section 5.

Research Area and Period
This study is focused on the BTH region ( Figure 1) for mainly two reasons. One is that BTH is one of the regions with the most serious air pollution in China and thus with typical representativeness. The second is that BTH is the main service area of the regional NWP model (RMAPS) operated by the author's team, whose data of real-time forecast and observations are readily accessible. condition and temperature prediction skills. Finally, the main conclusions are summarized in Section 5.

Research Area and Period
This study is focused on the BTH region ( Figure 1) for mainly two reasons. One is that BTH is one of the regions with the most serious air pollution in China and thus with typical representativeness. The second is that BTH is the main service area of the regional NWP model (RMAPS) operated by the author's team, whose data of real-time forecast and observations are readily accessible.
In view of the main research objective, it is important that the research period or case should be featured with the dramatically and unexpectedly changes (increase or decrease) of atmospheric pollutant concentration in a short time. Due to the strict industrial production and transportation restrictions and quarantines amid the sudden outbreak of coronavirus disease , the air pollutant concentrations in BTH significantly dropped in February than those in January 2020 [28,29]. Therefore, the comparative analysis of this study is concentrated in January and February 2020.

NWP Forecast Data and Observations
RMAPS is an operational regional NWP model and has been running on the highperformance computer of Beijing Meteorological Bureau since 2015. It was developed by the Institute of Urban Meteorology in the China Meteorological Administration in collaboration with the National Center for Atmospheric Research (NCAR). RMAPS is built based on the advanced research WRF (ARW) model and WRF data assimilation (WRFDA) of version 3.8.1 [33,34]. In addition to the temperature forecast data at different lead times from the RMAPS operations, the daily updated temperature forecast data from European Centre for Medium-Range Weather Forecasts (ECMWF, also use EC for short in the text) were also examined for comparison.
The hourly temperature observation data and other meteorological variables (e.g., relative humidity, wind speed and cloud cover) during January-February 2020 from 44 synoptic meteorological stations were obtained from China Meteorological Administration. Moreover, the hourly concentrations of air pollutants (including PM2.5, SO2, NO2, and CO) at 76 environmental stations operated by the Ministry of Environmental Protection In view of the main research objective, it is important that the research period or case should be featured with the dramatically and unexpectedly changes (increase or decrease) of atmospheric pollutant concentration in a short time. Due to the strict industrial production and transportation restrictions and quarantines amid the sudden outbreak of coronavirus disease (COVID-19), the air pollutant concentrations in BTH significantly dropped in February than those in January 2020 [28,29]. Therefore, the comparative analysis of this study is concentrated in January and February 2020.

NWP Forecast Data and Observations
RMAPS is an operational regional NWP model and has been running on the highperformance computer of Beijing Meteorological Bureau since 2015. It was developed by the Institute of Urban Meteorology in the China Meteorological Administration in collaboration with the National Center for Atmospheric Research (NCAR). RMAPS is built based on the advanced research WRF (ARW) model and WRF data assimilation (WRFDA) of version 3.8.1 [33,34]. In addition to the temperature forecast data at different lead times from the RMAPS operations, the daily updated temperature forecast data from European Centre for Medium-Range Weather Forecasts (ECMWF, also use EC for short in the text) were also examined for comparison.
The hourly temperature observation data and other meteorological variables (e.g., relative humidity, wind speed and cloud cover) during January-February 2020 from 44 synoptic meteorological stations were obtained from China Meteorological Administration. Moreover, the hourly concentrations of air pollutants (including PM 2.5 , SO 2 , NO 2 , and CO) at 76 environmental stations operated by the Ministry of Environmental Protection of the People's Republic of China were also used in the study. Most of these environmental and meteorological stations are located in the urban areas of 13 cities in the BTH region ( Figure 1

Method
Commonly used statistical methods, such as contrastive analysis and composite analyses with a two-tailed Student's t-test were applied in this research. In order to simply characterize the deviation or bias of the prediction from the observed value, we defined a PMO index, namely PMO = T prediction − T observation (1) Here, T predicton and T observation indicate the predicted and observed temperatures, respectively. The significance of the PMO index is that, positive PMO indicates the predicted temperature is higher than the observed, namely the NWP model forecast overestimates the temperature. On the contrary, negative PMO means that the predicted temperature is lower than the observed, namely the NWP model forecast underestimates the temperature.
And then a dPMO index is designed to denote the difference of the mean PMO in February and January, dPMO = PMO February − PMO January (2) Here, PMO January and PMO February indicate the mean PMO in January and February, respectively. The dPMO index reflects the difference of the model's temperature forecast bias between February and January. When dPMO > 0, the forecast deviation in February is greater than that in January. On the contrary, if dPMO < 0, the forecast deviation in February is less than that in January.
In order to understand the change of forecast performance more intuitively, we further calculated the root mean square error (RMSE) between the forecast temperature and observation for the daily mean, maximum and minimum temperature at different lead times, respectively. The calculation formula of RMSE is where T observation and T prediction denote the observed and predicted temperatures, respectively. The difference of RMSE can comprehensively reflect the variation characteristics of temperature prediction skill in January and February in the BTH region, without being affected by positive or negative phase of the deviation. In addition, a linear fitting equation based on the least square method is used to reflect the response function between the predicted temperature and the observation. The equation is as follows: where β is a constant, and α the regression coefficient. Generally, the closer to 1 α is, the better the prediction effect is.

Characteristics of Pollutant Concentration Changes over the BTH in January and February 2020
Based on the hourly pollutant records at the 76 environmental stations, the monthly mean PM 2.5 , SO 2 , NO 2, and CO concentrations in BTH in January and February 2020 were calculated. The results show that the mean PM 2.5 concentrations in BTH were 99.0 µg/m 3 (varying from 32.7 µg/m 3 to 150.8 µg/m 3 among the 13 cities) and 60.9 µg/m 3 (varying from 39.3 µg/m 3 to 79.6 µg/m 3 ) in January and February, respectively. The mean PM 2.5 concentration of these cities in February was 38.5% lower than that in January. There were few differences in the change of pollutant concentration between daytime and nighttime. The variation characteristics of SO 2 and CO concentration were basically the same as that of PM 2.5 concentration ( Table 1). The mean SO 2 and CO concentrations of these cities in February were 35.1% and 33.1% lower than those in January, respectively. The reduction in NO 2 concentration was more dramatical, with the mean concentration of the 13 cities dropped by 48.0%. Moreover, the decrease of NO 2 concentration in daytime (−50.2%) was significantly greater than that in nighttime (−46.2%). This is likely because the NO 2 concentration depends more on vehicle exhaust emissions, which is more concentrated in the daytime due to people's daily life. To understand the general variation characteristics of the pollutant concentration in the BTH region during the period, we first examined the concentrations of PM 2.5 , SO 2 , NO 2 and CO in the 13 cities during the period January to February (Figure 2), respectively. Firstly, combined with the spatial distribution of the 13 cities in BTH, it is shown that the high concentrations were mainly concentrated in central and southern BTH in both January and February. Secondly, from the perspective of day-by-day variations, most of the cities in BTH region had a significant decrease in pollutant concentrations in February compared with January, especially the concentration of PM 2.5 ( Figure 2a) and NO 2 ( Figure 2b). Among them, the concentration of NO 2 decreased the most, reaching −48.0%. As far as daytime and nighttime are concerned, the percentage or magnitude of pollutant concentration changes in February is roughly the same as that in January, although there are some differences or exceptions (Table 2). One exception, for example, is that the PM 2.5 concentration in Beijing in February 2020 was slightly higher than that in January. In addition, both of PM 2.5 and CO concentrations in Zhangjiakou in February 2020 were higher than those in January. These might be attributed to the heavy pollution caused by adverse weather conditions, which are conducive to the accumulation of pollutants and the formation of secondary aerosols and then causing an increase in PM 2.5 .    Generally, despite some differences existing in the spatial distribution and magnitude range of change, the mean PM 2.5 and other pollutant concentrations in both daytime and nighttime dropped dramatically in February due to the strict quarantine or lockdown amid the outbreak of COVID-19. This is consistent with the results of other studies [27,28]. Therefore, the significant decrease of air pollutants in the BTH region in February 2020 compared with January can be considered as an unexpected change.
For an operational NWP system, the great difference of aerosol concentration in a short time constitutes the different background conditions. Next, we will examine whether it can cause some variation of the temperature prediction skills based on the regional NWP system of RMAPS and the global NWP system of ECMWF. Figure 3 examines the average PMO of the daily mean temperature (Tmean), maximum temperature (Tmax) and minimum temperature (Tmin) with the forecast lead of 24 h (24 h for short) at the 13 cities in BTH in January and February. In January, the spatial distribution of Tmean PMO showed large spatial variability ( Figure 3a); Tmax was generally underestimated except at the two southernmost cities (Figure 3b) while Tmin was overestimated at all 13 cities (Figure 3c). Concretely, the PMO of Tmax in January varied from −2.17 to 0.23 • C with an average of −1.15 • C. In contrast, the PMO of Tmin in January varied from 0.06 to 3.36 • C with an average of 1.46 • C.   Table 2 presents the root mean square error (RMSE) between the predicted temperature and the observed temperature to eliminate the effect of bias signs. The RMSE of Tmean, Tmax and Tmin in January were 1.31 • C, 2.50 • C and 2.37 • C, respectively, while they were 1.08 • C, 2.26 • C and 1.86 • C in February. That is to say, the RMSE of the daily mean, maximum and minimum temperature prediction based on RMAPS model in BTH in February 2020 were 17.3%, 9.8% and 21.6% lower than that in January. Generally, the absolute deviation and the RMSE of the maximum and minimum temperature prediction in February were smaller than those in January. Moreover, for the forecasts with longer lead time (i.e., 48 h and 72 h forecast in advance), this feature is still evident ( Table 2).

Differences in Temperature Forecast Skills between January and February 2020
The PMO of Tmean, Tmax and Tmin in the 24 h forecast at the 13 cities in BTH in January and February based on the EC model was also examined. As shown in Figure 4, the systematic forecast bias of Tmax in January based on the EC model is quite different from that based on the RMAPS model. However, the overestimation of Tmin in January is in agreement with RMAPS, except for the underestimation of several cities in the south. Statistically, the RMSEs of Tmean, Tmax and Tmin based on the EC model were 1.08 • C, 2.02 • C and 2.09 • C in January, respectively, and then decreased to 0.88 • C, 1.73 • C and 1.65 • C in February (i.e., a reduction by 18.5%, 14.5% and 20.9%), respectively. As the forecast time extended to 48 and 72 h, this phenomenon was also still evident, although there are slight differences in degree ( Table 2). That is, the distinct comparison characteristics of the bias and RMSE in the two months are generally consistent with those in the RMAPS model.   In order to better understand the NWP performance change of temperature forecast skills in February 2022 relative to that in January and its physical significance, we further examined the dPMO of the Tmean, Tmax and Tmin with the forecast lead of 24 h (Figure 5), 48 h (Figure 6), and 72 h (Figure 7) based on both the RMAPS and EC models. First of all, it can be seen that all cities are dominated by negative dPMO of Tmin in the BTH region, and most of them are significant at 95% confidence level (Figure 5c). However, most of the cities in BTH are covered by the positive dPMO of Tmax (Figure 5b). For the mean temperature, almost all the cities are featured by negative dPMO (Figure 5a), but the absolute values are generally lower than that in the minimum temperature. The negative (positive) dPMO indicates that the PMO in February is smaller (larger) than that in January, namely the RMAPS-based temperature forecast deviation was decreased (increased) in February.
value of dPMO in Figure 5c indicates that, compared with January, the overestimation of the Tmin by RMAPS was decreased significantly in February. That is, the prediction skill of daily minimum temperature in February based on RMAPS has been improved. Meanwhile, the prediction of Tmax based on RMAPS in BTH in January was basically underestimated ( Figure 3b) and was improved generally in February, especially in the central part of BTH region (Figure 5b). The spatial distribution characteristics and the comparisons of PMO and dPMO are basically similar when the forecast lead increases to 48 h ( Figure 6) and 72 h (Figure 7).      Figure 5).

Figure 7. The dPMO of the Tmean, Tmax and Tmin derived from RMAPS (a-c) and EC (d-f) in 72
h forecast lead time in advance (same as Figure 5).

Figure 7. The dPMO of the Tmean, Tmax and Tmin derived from RMAPS (a-c) and EC (d-f) in 72 h
forecast lead time in advance (same as Figure 5).
On the other hand, the obvious difference between Figures 3 and 5 is that the PMO of Tmin in the whole BTH region is positive (Figure 3c), while the dPMO of Tmin is negative ( Figure 5c). As mentioned above, the positive of PMO in Figure 3 suggests that the RMAPS model overestimated the minimum temperature in January. Thus, the negative value of dPMO in Figure 5c indicates that, compared with January, the overestimation of the Tmin by RMAPS was decreased significantly in February. That is, the prediction skill of daily minimum temperature in February based on RMAPS has been improved. Meanwhile, the prediction of Tmax based on RMAPS in BTH in January was basically underestimated ( Figure 3b) and was improved generally in February, especially in the central part of BTH region (Figure 5b). The spatial distribution characteristics and the comparisons of PMO and dPMO are basically similar when the forecast lead increases to 48 h ( Figure 6) and 72 h ( Figure 7). Figures 8-10 further examine the difference of forecast skills of Tmax and Tmin in BTH in January and February based on RMAPS using scatter diagrams and fitting equations. It can be seen that, in both January and February, the dispersion of forecast and observation increased with the extension of the forecast lead time. That is to say, with the extension of the forecast lead time, the forecast skills of the model for both Tmax and Tmin were generally decreased. However, it can be quantitatively known from the fitting equations whether Tmax or Tmin and the regression coefficients in February are closer to 1 and greater than those of January. Taking Tmin as an example, for the forecast lead of 24 h, 48 h and 72 h, the regression coefficients of the fitting equation in February are 0.83, 0.77 and 0.68, respectively, significantly larger than the counterparts in January (0.77, 0.67 and 0.53). This undoubtedly shows that the forecast skill of the RMAPS model for the daily maximum and minimum temperatures in the region in February was better than that in January. h, 48 h and 72 h, the regression coefficients of the fitting equation in February are 0.83, 0.77 and 0.68, respectively, significantly larger than the counterparts in January (0.77, 0.67 and 0.53). This undoubtedly shows that the forecast skill of the RMAPS model for the daily maximum and minimum temperatures in the region in February was better than that in January.     We also examine the model's forecast performance of the temperature diurnal variation in January and February. Firstly, the observed and RMAPS-predicted temperature diurnal variations averaged over the 13 cities in BTH in January and February are presented for different forecast lead times in Figure 11a, Figure 11b and Figure 11c, respectively. Shown in Figure 11a, an obvious feature is that the RMAPS forecast temperature in January was significantly higher at 0-10 a.m. (local time) and lower at 1-7 p.m. than the observation. In contrast, the forecast and observed temperature diurnal variations in February are intuitively consistent, and the deviations at all times are much smaller than those in January. With the extension of the forecast lead (Figure 11b,c), although there is a slight difference in the magnitude of the specific changes, this phenomenon was still evident. We also examine the model's forecast performance of the temperature diurnal variation in January and February. Firstly, the observed and RMAPS-predicted temperature diurnal variations averaged over the 13 cities in BTH in January and February are presented for different forecast lead times in Figure 11a, Figure 11b and Figure 11c, respectively. Shown in Figure 11a, an obvious feature is that the RMAPS forecast temperature in January was significantly higher at 0-10 a.m. (local time) and lower at 1-7 p.m. than the observation. In contrast, the forecast and observed temperature diurnal variations in February are intuitively consistent, and the deviations at all times are much smaller than those in January. With the extension of the forecast lead (Figure 11b,c), although there is a slight difference in the magnitude of the specific changes, this phenomenon was still evident. We also examine the model's forecast performance of the temperature diurnal variation in January and February. Firstly, the observed and RMAPS-predicted temperature diurnal variations averaged over the 13 cities in BTH in January and February are presented for different forecast lead times in Figure 11a, Figure 11b and Figure 11c, respectively. Shown in Figure 11a, an obvious feature is that the RMAPS forecast temperature in January was significantly higher at 0-10 a.m. (local time) and lower at 1-7 p.m. than the observation. In contrast, the forecast and observed temperature diurnal variations in February are intuitively consistent, and the deviations at all times are much smaller than those in January. With the extension of the forecast lead (Figure 11b,c), although there is a slight difference in the magnitude of the specific changes, this phenomenon was still evident.
The second row in Figure 11 shows the diurnal variation of PMO (T prediction − T observation ) in January and February at the different forecast lead times. It is clear that the RMAPS model overestimated the nighttime temperature and underestimated the daytime temperature in January. The turning points were around 11 o'clock in the morning and 20 o'clock in the night, respectively. In February, the RMAPS model mainly underestimated the average temperature of the 13 cities throughout the day, with only a few positive deviations around 9 o'clock in the morning. The deviation mode in February is significantly different from that in January. It is worth noting that the absolute deviation of the temperature forecast in February is much smaller than that in January, especially at night.
Meanwhile, it is interesting to note that, regardless of whether the forecast period is 24 h, 48 h, or 72 h, the positive dPMO only appeared from 13 to 17 or 18 o'clock in the afternoon, and the rest was dominated by negative dPMO. Here, the positive values reflect the improvement of the underestimation of the daytime temperature forecast, while negative values indicate the improvement of the overestimation of the nighttime temperature forecast. That is, the diurnal variation characteristics of dPMO also suggested that the forecast skill of the RMAPS model for temperature over the BTH region in February was significantly improved than that in January. These results are basically consistent with those based on the EC model. afternoon, and the rest was dominated by negative dPMO. Here, the positive values flect the improvement of the underestimation of the daytime temperature forecast, wh negative values indicate the improvement of the overestimation of the nighttime temp ature forecast. That is, the diurnal variation characteristics of dPMO also suggested th the forecast skill of the RMAPS model for temperature over the BTH region in Februa was significantly improved than that in January. These results are basically consiste with those based on the EC model.  Based on the differences in bias, RMSE, spatial distribution and diurnal variation characteristics of the temperature forecast in January and February over the BTH region, it can be concluded that the temperature prediction skills in the BTH region, based on both RMAPS and EC models, were generally improved in February 2020 compared with those in January, especially for the minimum temperature prediction in nighttime, although the changes are model-and city-dependent to some extent. In order to understand the possible association between the improved temperature prediction skills and the reduced aerosol pollution, the changes and climatology characteristics of temperature, relative humidity, wind speed, cloud cover in January and February in BTH are discussed in the following.

Discussion
Why did the NWP performance in the temperature forecast over the BTH region differ greatly between January and February 2020? Generally speaking, in a relatively short period, without major upgrades or changes in an operational NWP system, the ability to predict a certain meteorological element should be basically the same. The significant change in the forecast skill might be related to the abnormal external influence. Therefore, it is necessary to explore whether this phenomenon is related to the response of meteorological conditions to the unexpected and significant changes in air pollutant concentrations. Figure 12 examines the mean air temperature at 2 m (T2) of January and February in 2020 and 2010-2019 and their differences are based on the hourly observations at the 44 synoptic meteorological stations located in the 13 cities. Whether in all day time (A), daytime (D) or nighttime (N), the climatological temperature in February and the temperature in February 2020 are higher than that of January. For the climatological temperature, the mean daytime (nighttime) T2 in January is −1.6 (−5.5) • C with the standard deviation of 1.9 (1.7) • C; while the mean daytime (nighttime) T2 in February is 1.9 (−2.4) • C with the standard deviation of 1.4 (1.3) • C, respectively. For the temperature in 2020, the mean daytime and nighttime T2 were 0.2 and −4.0 • C in January, and 3.8 and −0.2 • C in February, respectively. The increase of the monthly mean temperature in February relative to January is a climatic feature due to the seasonal variations of the solar radiation. On the other hand, both daytime and nighttime temperatures in January and February 2020 are significantly higher than the climatic means, which is likely to be attributed to the background of global and regional climate warming. CLOCov in February were higher than those in January (gray bars in Figure 13).
As for the actual observed cloud cover in 2020, the cloud cover was slightly larger during the daytime and smaller during the nighttime in February than that in January (red dots in Figure 13). However, from the comparison of monthly anomalies (blue squares in Figure 13), the daytime clouds and night clouds in February 2022 were obviously fewer than those in January 2020. The possible mechanism of the cloud on the surface temperature is likely consistent with that of aerosols-temperature relations, although there are a lot of uncertainties in aerosol and cloud feedbacks [35,36]. Moreover, the increases of northerlies in nighttime and southerlies in daytime can cause a cooling in nighttime and warming in daytime, respectively (Figure 13b). And a reduction of RH in nighttime can be also beneficial to the increase of OLR and the surface cooling ( Figure  13c).
Generally, the anomalous warming in daytime and cooling in nighttime in February 2020 are well consistent with changes of aerosol and meteorological factors in BTH, although the causal relationship between temperature variability and meteorological conditions is complex and highly uncertain. In theory, the effects of air pollutant concentrations on temperature can be achieved by an online-coupling meteorology-chemistry model, such as the WRF-Chem. Therefore, based on the online-coupling meteorology-chemistry model and the near authentic pollutant emission data, and with the assimilation of aerosol or not, a series of sensitivity numerical simulations are required to uncover the possible mechanism behind this phenomenon and improve the temperature forecasting skills. Figure 12. Then, mean T2 in January (a) and February (b) in 2020 (red dots) and their climatological values (grey histograms), and the differences of T2 in January and February (c) (A, D and N denote All day long, Daytime and Nighttime, respectively; Jan and Feb are abbreviations of January and February, respectively; Feb-Jan means the temperature difference between February and January; the red dots represent the temperature in 2020; the blue squares indicate the temperature anomaly of 2020 relative to the climatic state; the error bars of the gray histogram represent temperature variance of the 13 cities in BTH). Figure 12. Then, mean T2 in January (a) and February (b) in 2020 (red dots) and their climatological values (grey histograms), and the differences of T2 in January and February (c) (A, D and N denote All day long, Daytime and Nighttime, respectively; Jan and Feb are abbreviations of January and February, respectively; Feb-Jan means the temperature difference between February and January; the red dots represent the temperature in 2020; the blue squares indicate the temperature anomaly of 2020 relative to the climatic state; the error bars of the gray histogram represent temperature variance of the 13 cities in BTH). Therefore, it is necessary to exclude the effects of seasonal variations and climate warming when examining the differences of the temperature anomalies in January and February. Concretely, the anomaly daytime T2 of January and February 2020 were 1.63 and 2.30 • C, respectively, while the anomaly nighttime T2 of January and February 2020 were 1.54 and 1.50 • C, respectively. The anomaly daytime T2 deviation from February to January in 2020 was 0.67 • C, and the anomaly nighttime T2 deviation from February to January in 2020 was −0.04 • C.
The variations of the temperature anomaly in February and January are probably related to the changes in synoptic-scale weather factors, including relative humidity (RH), wind speed (WS), cloud cover (CLOCov), as well as chemical composition of the atmosphere. In addition, the variation of surface ground temperature in January and February 2020 also showed the similar comparative characteristics. The relatively high daytime temperature anomaly in February 2020 is likely related to the significant reduction of air pollutant concentration caused by epidemic controls [27,29]. A possible mechanism is that, the improvement of air quality largely means the decrease of aerosol concentration (i.e., PM 2.5 ) in the atmosphere and thus less scattering and reflection of solar shortwave radiation, which is conducive to the increase of near-surface temperature. Although the temperature anomaly at night in February changed little compared with that in January, only −0.04 • C, it seemed to be consistent with the interaction mechanism of air pollution and temperature. Due to the decrease of aerosols, the blocking and absorption effects of the near-surface atmosphere on the outgoing longwave radiation (OLR) from ground at night will be weakened, thus reducing the ground temperature at night.
The anomaly of temperature may also be related to relative humidity, wind speed, cloud cover, and other factors. Here, the variations of RH, WS, and CLOCov in daytime and nighttime in January and February are also examined, respectively. In terms of the climatological statistics, the daytime RH in February in BTH is lower than that in January, and the change at night is not obvious; and both the daytime and nighttime WS and CLOCov in February were higher than those in January (gray bars in Figure 13).

Conclusions
Due to the strict industrial production and transportation restrictions and quarantines amid the outbreak of coronavirus disease (COVID-19), the air pollutant concentrations in BTH significantly dropped in February 2020 compared with that in January. This study is to explore whether this sudden significant background change will affect the NWP performance prediction. Based on a large number of air quality and meteorological observations, as well as the forecast data derived from a regional (RMAPS) and global (ECMWF) real-time operational NWP models, a series of quantitative and qualitative analyses were performed. The main conclusions are summarized as follows.
(1) The air pollutant concentrations over the BTH in February 2020 were dramatically lower than that in January. For example, the concentrations of PM2.5 and NO2 were decreased by 38.5% and 48.0%, respectively. Simultaneously, after removing the climatological signal, the BTH region experienced an anomalous warming in daytime and an anomalous cooling in nighttime in February, respectively. The warming and cooling trends are generally consistent with changes of cloud cover, wind speed and relative humidity.
(2) The forecast performance of NWP models for daily maximum (Tmax) and minimum (Tmin) temperature over the BTH in February 2020 was significantly improved compared with that in January. In 24 h forecast lead time in advance, the root mean square errors (RMSE) of Tmax and Tmin prediction based on RMAPS in February were 9.8% and 21.6% lower than those in January, respectively. This phenomenon basically existed when the forecast lead time extended to 48 h and 72 h. And a similar result was also presented in the operational global ECMWF model.
(3) It is suggested that the temperature forecasting skill based on NWP models in February 2020 was generally improved, and this might be partly attributed to the significant reduction of aerosol pollutant concentrations in a short period of time, although the specific effect of aerosols on temperature prediction is likely to be regional-and modeldependent. The possible links deduce that less aerosols in daytime are beneficial to more solar radiation and a warming in the near-surface temperature. Contrarily, less aerosols in nighttime are not conducive to stop the outgoing long-wave radiation and lead to a cooling near the ground. The response mechanism of these meteorological conditions to the changes of aerosol concentrations is likely to be the reason for the improvement of NWP model's temperature prediction skills in February. Therefore, a further study with a series of sensitivity numerical simulations based on an online-coupling meteorologychemistry model is warranted for better understanding of the detailed mechanism. As for the actual observed cloud cover in 2020, the cloud cover was slightly larger during the daytime and smaller during the nighttime in February than that in January (red dots in Figure 13). However, from the comparison of monthly anomalies (blue squares in Figure 13), the daytime clouds and night clouds in February 2022 were obviously fewer than those in January 2020. The possible mechanism of the cloud on the surface temperature is likely consistent with that of aerosols-temperature relations, although there are a lot of uncertainties in aerosol and cloud feedbacks [35,36]. Moreover, the increases of northerlies in nighttime and southerlies in daytime can cause a cooling in nighttime and warming in daytime, respectively (Figure 13b). And a reduction of RH in nighttime can be also beneficial to the increase of OLR and the surface cooling ( Figure 13c).
Generally, the anomalous warming in daytime and cooling in nighttime in February 2020 are well consistent with changes of aerosol and meteorological factors in BTH, although the causal relationship between temperature variability and meteorological conditions is complex and highly uncertain. In theory, the effects of air pollutant concentrations on temperature can be achieved by an online-coupling meteorology-chemistry model, such as the WRF-Chem. Therefore, based on the online-coupling meteorology-chemistry model and the near authentic pollutant emission data, and with the assimilation of aerosol or not, a series of sensitivity numerical simulations are required to uncover the possible mechanism behind this phenomenon and improve the temperature forecasting skills.

Conclusions
Due to the strict industrial production and transportation restrictions and quarantines amid the outbreak of coronavirus disease (COVID-19), the air pollutant concentrations in BTH significantly dropped in February 2020 compared with that in January. This study is to explore whether this sudden significant background change will affect the NWP performance prediction. Based on a large number of air quality and meteorological observations, as well as the forecast data derived from a regional (RMAPS) and global (ECMWF) real-time operational NWP models, a series of quantitative and qualitative analyses were performed. The main conclusions are summarized as follows.
(1) The air pollutant concentrations over the BTH in February 2020 were dramatically lower than that in January. For example, the concentrations of PM 2.5 and NO 2 were decreased by 38.5% and 48.0%, respectively. Simultaneously, after removing the climatological signal, the BTH region experienced an anomalous warming in daytime and an anomalous cooling in nighttime in February, respectively. The warming and cooling trends are generally consistent with changes of cloud cover, wind speed and relative humidity.
(2) The forecast performance of NWP models for daily maximum (Tmax) and minimum (Tmin) temperature over the BTH in February 2020 was significantly improved compared with that in January. In 24 h forecast lead time in advance, the root mean square errors (RMSE) of Tmax and Tmin prediction based on RMAPS in February were 9.8% and 21.6% lower than those in January, respectively. This phenomenon basically existed when the forecast lead time extended to 48 h and 72 h. And a similar result was also presented in the operational global ECMWF model.
(3) It is suggested that the temperature forecasting skill based on NWP models in February 2020 was generally improved, and this might be partly attributed to the significant reduction of aerosol pollutant concentrations in a short period of time, although the specific effect of aerosols on temperature prediction is likely to be regional-and model-dependent. The possible links deduce that less aerosols in daytime are beneficial to more solar radiation and a warming in the near-surface temperature. Contrarily, less aerosols in nighttime are not conducive to stop the outgoing long-wave radiation and lead to a cooling near the ground. The response mechanism of these meteorological conditions to the changes of aerosol concentrations is likely to be the reason for the improvement of NWP model's temperature prediction skills in February. Therefore, a further study with a series of sensitivity numerical simulations based on an online-coupling meteorology-chemistry model is warranted for better understanding of the detailed mechanism.