Downscaling Winter Daily PM2.5 Concentrations From Large-Scale Meteorological Dataset -a Case Study for a City in North China


 Winter air pollution in North China becomes a serious environmental problem in recent years, which has aroused a widespread concern. Estimating PM2.5 concentration is necessary for the government to take actions in leading times, or to reproduce the historical values. In this study, we attempt to construct statistical downscaling (SD) models based on large-scale meteorological variables, to estimate the PM2.5 in Jiaozuo, a city in the heavy-pollution area of North China. Predictors were screened from large-scale meteorological variables, by selecting the grid boxes with the values highly correlated to the PM2.5 in Jiaozuo. Correlation maps show that PM2.5 is usually related to comparatively low pressure and high relative humidity at the lower atmosphere and comparatively high pressure at the upper atmosphere, high air temperature in all the troposphere and weak winter monsoonal winds. After training the SD models with the winter samples during 2014-2017, the daily PM2.5 is simulated with correlation coefficients of 0.607 (mean) and 0.548 (maximum) to the observation, by the logarithmic transformed PM2.5 values. The SD models can roughly reflect the daily variation of PM2.5. Compared to the PM2.5 model based on the local meteorological data, larger correlations are obtained (0.57 versus 0.51-0.53, without logarithmic transformation). In the days of “APEC Blue Sky” and the “SCO Blue Sky” with emissions largely reduced by strict controls on industry and traffic in surrounding areas, low PM2.5 is indeed observed, however, the SD models without considering any change of emissions also produced low PM2.5 simulations, implying that the large-scale circulation plays a main role in the daily variation of air pollution. Similar effects were obtained for the traffic-control month and the Chinese Spring festival with massive fireworks burning. Basing on the large-scale reanalysis variables and assuming the same emission level as that in the model training period, the downscaled winter PM2.5 has a significant decreasing trend during 1979-2017.


Introduction
The impact of the changing climate on environment is a public concern around the world today. A great amount of future climate predictions for the globe were produced by general climate models (GCM) in many organizations. However, these predictions are in rough resolutions which cannot meet the needs of local impact studies in different areas. Downscaling is useful to resolve such scale discrepancy. There are two categories of downscaling techniques can be used: dynamical downscaling and statistical downscaling(SD). Dynamical downscaling refers to the using of regional climate models to simulate the climate in the interested areas in higher resolutions, driven by the datasets of large-scale climate prediction which is simulated by GCMs. Dynamical downscaling can generate reasonable downscaled result but it is computationally expensive for downscaling long-term climate scenarios. SD refers to the technique that derives local variables from large-scale variables (LSV) in climate scenarios by constructing empirical relation between the large-scale variables and the local variables (Maraun et al., 2010). Compared to dynamical downscaling, SD is computationally very cheap. The Lucio-Eceiza et al., 2019). Nevertheless, there is no widely applicable SD models which can be directly used for anywhere. In other words, SD models need to be tailored for different sites/areas respectively. The techniques for constructing SD models, including screening predictors, reducing the dimensionality of predictors and selecting mathematical structure of models, were presented in many previous studies ( There have been many SD models for downscaling surface air temperature and precipitation now, but very few studies have been presented for one public concern in China: the concentration of PM2.5 (suspended particulate matter (PM2.5) with aerodynamic diameter less than 2.5 mm). Epidemiological studies have reported that exposure to PM2.5 has adverse effects on human health (Kunzli et Zhang et al., 2020). That means, under a given emission level of pollutants, the relative concentration level of PM2.5 can be estimated/predicted from the meteorological factors. Actually, at a certain area, the human activities generally keep stable during different years. The emission levels can be reduced year by year through the gradual innovation of technologies, but after the reduction reaches a certain limit, a minimum emission will be maintained. Thus, at a presumed constant emission level of pollutants, the climatic variation or trend of PM2.5 concentration can be estimated from pure meteorological factors.
Estimation of PM2.5 from pure meteorological factors seems practically meaningless at the daily scale, since there have been many accurate and operational methods for estimating daily PM2.5 concentration, such as numerical models or machine learning models (Analitis et al., 2020;Chen et al., 2019;Zheng et al., 2016). Most of these advanced high-accuracy models are highly relied on many localized datasets (such as the data of local meteorology, remote sensing and pollutant emissions) (Zhang et al., 2019) and can be used to reconstruct the daily/hourly PM2.5 concentration in the period with all these available datasets, however, such models is irrelevant to the climate science, and cannot be used to re ect/predict the relative PM2.5 levels in the periods with some of the datasets unavailable. The downscaling of PM2.5 concentration from large-scale meteorological variables is useful in certain aspects. For example, people may want to know the relative concentration level of PM2.5 which is caused by the meteorological conditions in the last ve decades in which most of the time there was no observational PM2.5 data. It is worth noting that there was no PM2.5 records in China before 2013, and the remote sensed AOD data is also only available for the last two decades. People may also want to know how the future PM2.5 concentration will change under the background of climate change. In such a sense, downscaling the relative level of PM2.5 concentration has the meaning in climate science just as that for the downscaling of precipitation, wind and air temperature. In weather prediction, the PM2.5 concentration can be generally well predicted by numerical models, but predicting PM2.5 through statistical downscaling based on large-scale circulation is also valuable, since this method is very simpler and computationally cheaper than numerical models. Till now, few studies have addressed such downscaling problem.
In this study, we made an attempt to construct SD models for downscaling daily PM2.5 concentration in winter at a site in North China.
The aim of this study is to investigate following related questions: whether the large-scale circulation representing variables can be used to get satisfying predictors; what large-scale circulations are related to the site PM2.5; what predictors can be identi ed from different large-scale variables; how the performance can be obtained by such SD models when compared with the models constructed from local ground based meteorological variables; what performance can be obtained in reproducing the daily PM2.5 concentration in the winters of recent years. As a research on downscaling based on purely meteorological large-scale variables, we do not expect the downscaled outputs be very accurate performance as that achieved by many studies on modelling PM2.5.

Study site
As a case study, we conduct this statistical downscaling of the PM2.5 in Jiaozuo, a city in North China Plain (NCP). The area in NCP, including the Beijing-Tianjin-Hebei region (BTH) and the North of Henan province, is one of the most populous regions in China, and has been the region having the most severe air pollution before 2014 (Zhang et al., 2012). To cope with the severe air pollution problem, the State Council of China established the toughest-ever air pollution prevention and Control Action Plan in 2013. By taking such effective measures, the air quality is improved year by year (Xiao et al., 2020). However, in order to maintain the sustainable economic development to guarantee the necessary employment, most of the industries have to be preserved. Thus, a least emission level will have to be last for a long time in the future before important technological innovations being available, and the air pollution problem cannot be completely resolved. Heavy air pollution is still frequent in winter times now, although the most serious pollution situation has been improved compared to previous years.
As a case study, we use the PM2.5 concentration in the urban area of Jiaozuo (Fig. 1). Jiaozuo is a medium-sized city (according the standard for the classi cation of cities in China) in the south part of NCP and belongs to the northwest of Henan Province. Jiaozuo is a neighbor city of Zhengzhou, the metropolis of Henan Province. The administration region of Jiaozuo has an area of 4071 km 2 and a population of 3.6 million. The urban area in the region includes the center urban area of Jiaozuo and other six county-level cities. The whole Jiaozuo administration area and its surrounding areas in the north of Henan province are highly industrialized. The large industry production in this area gives rise to large pollutant emissions, because heavy chemical industry accounts for a large proportion of gross domestic production (GDP) in the north of Henan Province. Due to the shielding effect of Taihang Mountain on atmospheric circulation and the high industrial activities in North China, a long belt of heavy air pollution stretches along the East and south side of Taihang Mountain in the NCP (Fig. 1). The area of Jiaozuo is just in the south of this belt, and Jiaozuo is also famous for being one of the most polluted cities in China in the last two decades. Due to the fact that the central urban area of Jiaozuo just lies at the foot of Taihang mountain, the pollution diffusion condition of the area is worse than the surrounding areas.

Data
In this study, we obtained the hourly PM2.5 concentration data recorded by three air quality monitoring stations (Huanbaoju, Gaoxinqu and Yingshicheng) in the city of Jiaozuo, which is managed by the Environmental Protection Bureau of Jiaozuo. Considering very few missing values exists in the PM2.5 observation at station Huangbaoju, we only used the data at this station to construct SD models. The data for stations Gaoxinqu and Yingshicheng were used as reference to compare or validate. The dataset for Huanbaoju covers the period 2014-2018, only the values in the rst several months of 2014 are missing. The mean and the maximum of each day were extracted from this hourly dataset. Both the daily mean PM2.5 and maximum PM2.5 were used as predictands for the statistical downscaling.
The ERA-Interim Reanalysis dataset (ERAI) from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Dee et al., 2011) were downloaded. The variables in ERAI with multiple pressure levels (500hpa, 700hpa, 850hpa and 1000hpa) were used to screen predictors: relative humidity (R), speci c humidity(Q), latitudinal direction wind (V), longitudinal direction wind (U), air temperature (T) and geopotential height (Z). Here, the ERA-Interim data forecasted at 0h UTC but for 12h UTC in every day were used to represent the meteorological status of that day. The selected dataset covers the period 1979-2018. This dataset has the domain of 60°E-150°E and 10°N-65°N, which covers the areas surrounding China. To facilitate the data processing, the dataset was degraded into a 1°×1°resolution grid when downloaded. Only the datasets from the winter months (November, December, January and Februry) were used. In this study, in order to preserve the temporal continuity of the samples in each winter, we take the samples in November and December of each year and those in January and February of the next year into one winter sample set. Therefore, the winter sample set of 2014 should also include the samples in January and February of 2015.
Considering there is no published ground-based meteorological data available in the administrative area of Jiaozuo, we also collected surface meteorological data (including air temperature, wind speed, air pressure, relative humidity, pan evaporation and sun shine duration) observed by Zhengzhou meteorological station, which is one of the nearest site (within 100km) to Jiaozuo, and both Jiaozuo and Zhengzhou shares generally the same weather. Accordingly, the PM2.5 concentration data observed by the Center Air Quality Monitoring Station of Zhengzhou during 2014-2018 was also collected, which is used to construct a simple model for estimating PM2.5 concentration based on ground based meteorological data, in order to compare the performance with other different models. Here, we do not plan to construct the downscaling model based on the large-scale circulation for Zhengzhou, considering the models for these cities should have theoretically little differences. The correlation maps can actually re ect the relation of the LSV's anomalies to the site PM2.5, and cannot re ect the true circulation over the area. In order to compare different circulations at high PM2.5 days and low PM2.5 days, we divide the winter daily samples into a high-concentration group and a low-concentration group by the threshold of PM2.5=38.5mg/m 3 , the medium PM2.5 value of the samples.
The geopotential height, U-wind and V-wind in each group were averaged over all the samples in the group, respectively.

Statistical downscaling Method
The statistical downscaling comprises three steps: screening predictors, training statistical models and assessing the model performance.
Predictor screening is vital to SD models, and the most commonly used LSVs for downscaling precipitation are related to air pressure, moisture and air circulation. For PM2.5, theoretically, the above variables are also applicable and were also used in this study. The pressure levels and the grid boxes of the LSVs were considered. In this study, we use a correlation map based predictor screening method. In this method, the correlation maps should be created by relating each candidate variables with the predictand, then the grid boxes with The LASSO based GLM tting was used to estimate the GLM parameters. The predictors screened from the multiple variables are also usually interior correlated, thus using all the predictors can inevitably cause over tting problem. However, there is no strict rules to use or discard each predictor by considering the correlation coe cient values among different predictors. In a standard over tting problem, the regression coe cient/parameter of a predictor has a wrong sign which is different to the sign of the correlation coe cient of this predictor to the predictand. For example, one predictor screened from the air pressure has positive correlation to the predictand, however, during the model tting process, this predictor gets a negative regression coe cient. This effect is caused by the high correlation of this air pressure based predictor to another predictor extracted from the air temperature (or air humidity): the signal of the air pressure based predictor is substituted by the signal of the air temperature based predictor, while a meaningless regression coe cient with wrong sign is obtained; the regression coe cients for both these two predictors can give rise to a wrong output when new samples other than the training samples are fed to the model. To avoid this over tting problem, the predictor based on air pressure should be discarded, and the functionality of the discarded air pressure predictor is played by the air temperature predictor. The Least Absolute Shrinkage and Selection Operator (LASSO) algorithm is such a method (Hammami et al., 2012). By using a threshold, LASSO can automatically determine which predictor to be discarded in the model tting process. If over tting problem still exists after a LASSO based model tting, then the threshold in LASSO should be increased to a larger value and re t the model.

Model validation
One-winter-leave-out cross validation was used to validate the models. There are only four years (2014-2017) of observed PM2.5 available, therefore, in each training-validation step, three years were used for training, and one year were used for validation. Rooted mean squared error (RMSE) and the Spearman's correlation coe cients were used to evaluate the model performance.
During the winter days of 2014-2017, there are some speci c periods in which signi cant emission reduction or increase was caused by some social activities. These speci c periods can be regarded as ideal scienti c experiments to study the mechanism of air pollution. The detailed information for these speci c periods is introduced here: 2) "SCO Blue Sky" (Dec. [14][15]2015). An important leaders meeting of the Shanghai Cooperation Organization (SCO) was held in Zhengzhou, Henan Province, during December 14-15, 2015. To ensure the good air quality, strong restrictions were executed in Henan province. As a city neighboring to Zhengzhou, Jiaozuo also got several days of good air quality.
3) Tra c control with even and odd-numbered license plates in multiple cities of Henan province. To avoid the happening of bad air quality in the last month of 2017, in each day, half of the cars were restricted to travel in the main urban areas of multiple cities (Jiaozuo, Zhengzhou, Kaifeng and Anyang) in December 4-31, 2017.
4) The Chinese Spring Festival. The Spring Festival (the rst day of each year in the Chinese Lunar Calendar) is well known as the very important traditional festival in China, and burning reworks are always a necessary activity in the New Year's Eve (the night before the Spring Festival). In recent years, more reworks are burning than in the past because the residents become richer than ever. This activity is bound to create severe pollution and every one can smell the strong sul de in air during this festival. Nevertheless, since the Spring Festival in 2016, burning reworks were banned in the main urban area of Jiaozuo to reduce air pollution, but no restriction was required in county-level urban areas and the wide rural areas. Similar restrictions were also executed in other cities in the surrounding areas of Jiaozuo.

Schemes for comparison
In order to analyze whether the GLM outputs are sensitive to different choices on using predictors, two different schemes for comparison are considered. The rst scheme is to analyze whether the predictors extracted from wind eld has large in uence on the output. The design of this scheme is due to the fact that multiple predictors are usually obtained from the wind eld, including U-wind and V-wind, which makes the process for screening predictors complex. The second comparison scheme is to use the single-grid box value at the high correlation center or the average value in multiple grid boxes of the high correlation area to get predictor values. Here, in order to obtain multiple grid box values, we determine each window through manual inspection on each correlation map. The sizes for the windows are different for different high-correlation areas.
Four combinations by the above two schemes were designed for comparison: 1) use the value of the single grid box at each highcorrelation center, and no wind based predictors used; 2) use the average value of multiple grid boxes in a certain window, and no wind based predictors used; 3) use the single grid box value at each high-correlation center, with wind based predictors used; 4) use the average value over the multiple grid boxes, with wind based predictor used.

Variance in ation
In statistical downscaling, the models usually produce outputs with smaller variance than the observed, because the random variability is ignored by the models. Such underestimation of variance is not suitable for risk warning. Variance in ation techniques (Maraun et al., 2010) can be used to enlarge the variance. In this study, we use a simple variance in ation: 2.4 Estimating PM2.5 based on local meteorological data Estimating PM2.5 from local ground based meteorological data is a common way. In this study, we also constructed PM2.5 models based on the local meteorological data. In following paragraphs, such models are represented as the 'local models', while the above models based on LSVs in ERAI are represented as the 'downscaling models'. Local models were constructed here is to gure out which model has the advantage. Local models based on GLMs were constructed to estimate the PM2.5 concentrations both at Jiaozuo (only Huanbaoju station) and Zhengzhou, based on the same meteorological data observed in Zhengzhou. The parameters were estimated by the LASSO-GLM method. Similarly, the one-winter-leave-out cross validation was also used here to obtain the Pearson's correlation coe cient (CC) between the model's output and the observed PM2.5, as a metric of the model's performance.

Statistics of PM2.5 concentration in Jiaozuo
Histograms for the daily mean and maximum PM2.5 in Jiaozuo show skewed distributions (Fig), indicating the Pearson's coe cient and common linear regression are not suitable to be used for PM2.5.Transformation is usually used to convert such skewed-distribution variables into more Gaussian distribution variables (Qian et al., 2015). Here, we use a logarithmic transformation.
The logarithmic transformed data seem to be more Gaussian like. All the following modelling process is based on these logarithmic transformed data. Although common linear regression can be used here, but we still prefer GLM considering that GLM is also suitable for the data of normal distribution.

Correlation maps
The CCs were calculated between the mean/max daily PM2.5 of Jiaozuo and the time series extracted at each grid box from the LSVs, then the correlation maps were produced. According to these maps, the high-correlation areas are not necessarily close to Jiaozuo or in its surrounding area.

Geopotential height
The correlation maps between the PM2.5 and the geopotential height show that heavily polluted days correspond to the negativecorrelation area at the Mongolian Plateau at the 500 hPa layer and the positive-correlation area stretching from the east of North China to the Korean Peninsula (Fig. 3). The negative-correlation area covering the south part of China and the positive-correlation area covering the Northwest Paci c and the Indian peninsula on the 850 hPa and 1000 hPa layers. This shows that the main circulation in heavy air pollution days of Jiaozuo is characterized by the relative low-pressure troposphere on the Mongolian Plateau and the relative high pressure from the east of North China to the Korean Peninsula, indicating that the cold air is weakly transported from the Mongolian Plateau to the south. The low pressure in the South Central China appears in the lower troposphere. At the location of Jiaozuo, the relative high pressure appears at the top of troposphere and the relative low pressure appears at the bottom of troposphere. The existence of relative high pressure over the Indian subcontinent and the Paci c Ocean indicates that the low-layer atmosphere tends to convergent circulation, thus, under such conditions, the pollution diffusion in eastern China is weak.
According to the correlation maps, there are one negative-correlation center (Mongolia-Siberia) and one positive-correlation center (North China-Korean Peninsula) at the 500hPa layer, while there is one negative correlation center (the south-eastern part of China) and two positive-correlation centers (Indian subcontinent and the Northwest Paci c). The geopotential at these centers were used as predictors to construct SD models. The correlation centers at 700hPa layer and 850hPa layer were not used here because these two layers can be regarded as the transition from the 500hPa layer to the 1000hPa layer and are not typical enough.

Relative humidity and speci c humidity
The correlation maps calculated from relative humidity and speci c humidity at the four pressure levels are generally similar and show a broken structure of high-correlation areas which can pass the signi cant tests at the 0.05 level. There are many such small patches of high or low correlation areas. Nevertheless, the broken structure obtained from speci c humidity is not apparent as that from relative humidity.
From the correlation map of relative humidity at 500hPa, there is a negative-correlation belt stretching from the Indian subcontinent to Southwestern China and a tiny high-correlation region in the lower Yellow River Basin. At the 700hPa level, there are negative-correlation areas in India, Northeastern China, Northwest Paci c, a large positive-correlation area at the Mongolia-Siberia region and a small positivecorrelation area between the lower reaches of the Huai River and the Yangtze River (the JiangHuai region). The speci c humidity at 700hPa has some similar patterns as the relative humidity at 700hPa, except that there is no apparent positive-correlation area at the Mongolia-Siberia region, and no negative-correlation area in Northeastern China. At the 700hPa level, the positive-correlation region obtained by speci c humidity at the Jiang-Huai region has a larger area than that obtained by the relative humidity. At the 850hPa/1000hPa level, the positive-correlation area covering the Jiang-Huai region becomes larger and stretches to the south of China.
The negative correlation center is also apparent in the west Paci c Ocean. The above patterns indicate that heavy air pollution in Jiaozuo is related to high humidity in the southeast of China at the bottom of the Troposphere. From lower levels to upper levels of the troposphere, the high-correlation region surrounding Jiaozuo gradually moves northward and becomes smaller.
We select the three high-correlation centers from the relative humidity at 700hPa as model predictors: the negative-correlation area at India subcontinent, the negative-correlation area in Northeastern China and the positive-correlation area in the Jiang-Huai region. Two correlation centers from the relative humidity at 1000hPa are also selected: the negative-correlation center at the Northwest Paci c, and the positive-correlation center at the eastern part of China.

Temperature
The correlation maps obtained from air temperature show a high-positive center in the east of China and a negative region in the Mongolia-Siberia region, indicating that heavy PM2.5 in Jiaozuo is related to the relatively warmer air temperature in the surrounding area of Jiaozuo and the relatively colder air temperature in the Mongolia-Siberia region, implying a stationary and weak diffusion condition for air pollution. Usually, in winter, when cold waves move through China, the air temperature in eastern China becomes colder, and the heavy air pollution can be blown away in a relatively short time. This circulation for cold-wave days is contrary to the above patterns of air temperature. From lower levels to upper levels in the troposphere, the high-correlation region surrounding Jiaozuo gradually moves northward, indicating the warm air mass moves to Eastern China earlier in the upper-level troposphere than in the lower-level troposphere. The positive-area of air temperature at 1000hPa in Eastern China seems broken or becomes weaker in a belt along the Taihang Mountain, and this weaker belt stretches southward to the Yangtze River basin. Considering all the above patterns, we selected both the centers in Mongolia-Siberia region and the eastern China as model predictors.

Wind
From the correlation maps from U-wind, a long belt of negative correlation stretches along the longitudinal direction, from North China to Japan, in 500hPa and 700hPa levels, indicating a relatively westward wind anomaly. In lower level troposphere (850hPa and 1000hPa), the negative area is connected with a large region in the Northwest Paci c, implying that the heavy air pollution in Jiaozuo is related to the westward wind anomaly widely existing in East Asia. Meanwhile, in the upper levels (500hPa and 700hPa), there is also a large and long positive-correlation area in the north of Asia continent, covering all the Central Asia, Mongolia, Northeastern China and Siberia, indicating a strong west wind anomaly in the upper-level troposphere. This pattern is not shown in lower levels, implying that it is related to the stable westerly jet.
In the upper levels of the U-wind, there is a large positive-correlation area from the whole area of East China to the Mongolia Plateau and Siberia, indicating a strong northward wind anomaly. Meanwhile a negative-correlation area exists in the Northwest Paci c, indicating a southward wind anomaly there. In the lower levels, this positive-correlation area becomes smaller and weaker, and the negative-correlation area seems moves westward and covers Northeastern China.
All the above patterns show the heavy air pollution is related to a low-pressure system with weakened wind blowing from the Mongolia-Siberia region to the surrounding areas in Eastern China. Under this atmospheric condition, the strong westerly in the subarctic area can prevent cold air mass from moving to mid-latitude East Asia. For wind variables, we select the negative center (South China) and the positive center (North China) of U850, V500 and V850 as predictors.

Main circulations in high/low PM2.5 days
By comparing the geopotential height at 500hPa in low-PM2.5 days and high-PM2.5 days, low-PM2.5 days re ect comparatively low pressure and strong northwestern wind over the southeast of China. In both high PM2.5 and low PM2.5 days, west wind dominates the area north to 35°N, indicating the stable westerly jet. Comparatively, in low-PM2.5 days, the wind in Southeastern China become more southerly.
At 850hpa, both in high-PM2.5 and low-PM2.5 days, there is a high pressure belt stretching from southwest to northeast in China, covering the Sichuan basin, Loess Plateau, North China and Northeastern China. A low-pressure area also exists in the Northwest Paci c. In the high-pressure belt, the area over North China in high-PM2.5 days poses a weak and divergent air ow to the surrounding areas.
Comparatively, in high-PM2.5 days, the pressure in this belt is not as high as that in low-PM2.5 days and the pressure in the Northwest Paci c is also not as low as that in low-PM2.5 days. Therefore, the land-sea pressure gradient during low-PM2.5 days are larger than that in high-PM2.5 days, producing a strong wind which is conducive to the dissipation of PM2.5 in the eastern part China.
Combining the geopotential heights at the upper level (500hpa) and the lower level(850hpa), it is obvious that high PM2.5 concentration usually corresponds to comparatively high pressure in upper air and low pressure in lower air over North China, which implies a static strati cation of atmosphere. In low PM2.5 days, on the contrary, comparatively low pressure in upper air and high pressure in lower air imply a strong downward air ow, which is usually related to the cold wave coming from the north pole. 3.3 Performance of the downscaling models

Parameter estimation and the nal models
The LASSO algorithm was used during the GLM training (parameter estimation), and some of the predictor parameters/coe cients were replaced by zeroes. The parameters were estimated based on four different sample sets, respectively and the RMSE and the CCs (based on the logarithmic transformed PM2.5) (Fig. 2-3). The predictor set including wind (15 predictor values) get better model performance than the predictor set without using wind variable (21 predictor values). Generally, the predictors extracted from averaged multiple grid boxes has poorer model performance than the single grid box based predictors. The differences between the simulated PM2.5 under the four combinations of the comparison schemes (the single grid box combining 15 predictors, the average of multiple grid boxes combining 15-predictor, the single grid box combining 21 predictors, the average of multiple grid boxes combining 21-predictor) are generally small and share the similar variations (Fig. 6).  Here, the predictors X 1− X 19 are listed in Table 1.
We use these models based on samples in 2015-2017 to produce results for further analysis, since relatively good model performance was obtained based on these samples. Two additional PM2.5 series observed at other two sites (Gaoxinqu and Yingshicheng) in Jiaozuo were also used to validate the models.
For the models for daily mean PM2.5, the CCs between the observations and the modeled results are 0.52 for both the sites (Gaoxinqu and Yingshicheng). The average of daily meanPM2.5 over the three sites in Jiaozuo has a CC of 0.55 to the modeled result. All the above CCs are lower than the CC obtained at Huanbaoju (0.56).
For the maximum PM2.5, the CCs obtained are 0.50, 0.48, 0.52 and 0.54, for Huanbaoju, Gaoxinqu and Yingshicheng and the average over the three stations, respectively. It is very interesting that the model was trained based on the data of station Huanbaoju, theoretically, the highest CC should be obtained for this station, but smaller CC was obtained instead than that for the three-station average and the station Yingshicheng. This demonstrates that the single station of Huanbaoju has some local variabilities which is related to some local conditions.
The simulated daily PM2.5 mean/maximum values has variations roughly corresponding to the observed values (Fig. 6). The simulated peaks and valleys usually have one day leading or lagging behind the observed ones. This inconsistency maybe related to the local atmospheric condition which affects the pollution transportation. It is common that the heavy PM2.5 episodes usually lasts for several days, and the simulated PM2.5 in this study can re ect such multi-day variations.

Variance in ation
The simulation underestimated the daily variances of both mean PM2.5 and maximum PM2.5 (Fig. 7). Some speci c peak values were not produced by the models. In order to match the simulated variances with the observed, variance in ation was used. The simulated standard deviation is 0.43 of the observed standard deviation. Considering the emission levels are different for each year, we use a different in ation factor for each year.
After this in ating, the magnitude of the variations can match the observed, however, very few peaks are produced by the simulation. The mean PM2.5 are simulated with better performance than the maximum PM2.5, and the underestimation of maximum PM2.5 is also more signi cant than that of the mean PM2.5. In contrast, the models perform the poorest in winter 2015 and perform the best in winter 2014. Similar good performance is obtained in 2017 as that in 2014.
After in ating the variance, the increase of RMSE is smaller than 0.01 µg/m 3 , while the relative biases are reallocated for the winter of each year, and the average absolute biases are reduced from 0.64 to 0.59 and from 0.72 to 0.57 respectively, for the mean PM2.5 and maximum PM2.5. The variance in ation does not change CCs, and causes very small increase of biases, but can improve the risk warning on peak PM2.5 days.

Outputs in special periods
In the special periods of the APEC Blue sky and SCO Blue sky, the simulated PM2.5 at these days are also comparatively low, indicating that these blue-sky days were mainly caused by good atmospheric circulation (Fig. 8). Therefore, as for whether the emission control measures have a signi cant impact on reducing PM2.5 concentrations in the two special periods, this study cannot give a conclusion. According to the experiments through numerical transport simulation, the emission control measures explains a PM 2.5 reduction of -21.8% in Beijing (Gao et al., 2017).
Similarly, in the December 2017, the even-odd-number tra c control in multiple urban areas shows no signi cant impact on PM2.5, since the observed PM2.5 also has several peak days corresponding to the downscaled results.
In During the Spring Festival of 2016, there is actually a peak observed, but it is lower than the corresponding peak simulated by the downscaling model ( Fig. 9). During the Spring Festival of 2017, although the reworks were banned in the main urban areas, however, burning reworks in other wide areas were still allowed and the emissions certainly has important effects on PM2.5. This effect is also not signi cant by comparing with the simulated PM2.5.

Downscaled annual PM2.5 during 1979-2017
In this study, the historical simulation for the winters of 1979-2017 was produced purely based on large-scale atmospheric predictors without considering the annual changes of emissions. The mean daily PM2.5 and maximum PM2.5 are shown Fig. 10. Mann-Kendall tests (Kendall, 1955;Mann, 1945) shown that both the simulated average PM2.5 mean and maximum have decreasing trends which can pass test at the 95% signi cance level, during 1979-2017. This result implies that under the same assumed anthropogenic and natural emission level, due to the large-scale atmospheric factors, the PM2.5 concentrations is generally dropping. This result is contrary to the recent observed trend which is largely caused by the increased pollutant emission. The dropping trend is mainly resulted from the lower values after 2005.
Theoretically, the relation between pan evaporation and PM2.5 is very complex. Nevertheless, the pan evaporation in high PM2.5 days is comparatively low because in hazy days, the high relative humidity and weak air ow can reduce water evaporation. In this study, the CC between the observed EVP (represented by the data observed at Zhengzhou station) and observed PM2.5 during winters of 2014-2017 is -0.44, while that between the observed EVP at Zhengzhou and the simulated PM2.5 is -0.42 during the same period. (R=-0.39 during 1979-2018). 3.3.5 Comparison between the downscaling models and the local meteorological data based models The local meteorological data based models for daily mean PM2.5 concentrations were trained by the LASSO-GLM algorithm. In all cases based on different sample sets of the cross validation, the parameters for air temperature and pan evaporation were set to zeros by the LASSO algorithm, indicating that the useful signals carried by these two variables can be represented by the signals in other variables.  Table 4). The values of R and RMSE obtained by these local meteorological based models are close to that obtained by the downscaling models.
The outputs of the downscaling models and the local meteorological data based models were shown in Fig. 11(a). It is clear that the observed daily PM2.5 values in Jiaozuo and Zhengzhou are generally consistent to each other, indicating that the air quality in these two cities are highly controlled by the same weather system and that using the meteorological data observed in Zhengzhou to represent the weather in Jiaozuo is reasonable. The outputs of the local model and the downscaling model produces similar variations although they are not completely the same, and both of them underestimated the observed values. Based on the comparison between the modeled outputs and the observations, the downscaling models for Jiaozuo get a larger R value than the PM2.5 model based on local meteorological data, and the R value obtained by the local model for Jiaozuo is slightly smaller than that for Zhengzhou.

Conclusion
In this study, a statistical downscaling model of the winter daily PM2.5 in a city (Jiaozuo, Henan Province) in North China was constructed based on large-scale meteorological variables. The predictors were screened based on the correlation maps between the large-scale variables in ERA-Interim reanalysis and the observed daily mean and maximum PM2.5 concentrations. Just as some statistical downscaling studies for precipitation, the high-correlation areas re ecting the large-scale circulations around the site is usually not coincident with the location of the site. Within the same large-scale variable, the high-correlation areas are different among multiple pressure levels. The high-correlation areas covering the city Jiaozuo indicate that high PM2.5 concentration is usually related to comparatively low pressure in the lower atmosphere and high pressure in the upper atmosphere, high relative humidity in the lower atmosphere, high temperature in all the troposphere and weak winter monsoon. Besides this, some areas far from the Jiaozuo location also have high correlations to the PM2.5 of Jiaozuo: the low geopotential height and low air temperature in the upper troposphere in the west Mongolia and Siberia, the high geopotential height in the lower troposphere around Japan and the Northwest Paci c, the low humidity at the India subcontinent and the Northwest Paci c, the west wind anomaly in the subarctic areas and the southward wind anomaly in the Northwest Paci c. This demonstrate that predictors for estimating PM2.5 can be selected from some positions far from the location of the observed predictand.
Among the multiple high-correlation areas, we selected some values of the large-scale variables as predictors, then trained GLMs using the LASSO-GLM algorithm to eliminate the over tting problem which is caused by the collinearity in the multiple predictors. According the cross validation, the predictor extracted from the multiple grid boxes has no advantage over the predictor extracted from the single grid box, especially for the daily maximum PM2.5. Meanwhile, the incorporation of predictors from wind variables can give rise to better performance than the predictors without using wind variables. Nevertheless, such different choices result in very small difference in the model outputs.
The standard deviation of the modeled PM2.5 explains 0.43 of the observed counterpart. The CCs between the model outputs and the observed PM2.5 (on the logarithmic transformed version) is 0.607 and 0.548, for the mean daily PM2.5 and the maximum daily PM2.5, respectively. Since the modeled outputs underestimated the variance, variance in ation was used to the outputs, then the relative bias of the modeled outputs was reduced signi cantly with a very small increasing of RMSE. Compared to the PM2.5 estimating models based on the local ground observed meteorological data, the downscaling models based on large scale circulations have a better performance.
During the periods of "APEC Blue Sky" and "SCO Blue Sky", the statistical downscaling model also reproduced the low PM2.5 levels, implying that the large-scale circulation plays a main role in these blue-sky days. The effect of emission control measures during these periods is di cult to quantify by comparing the observed PM2.5 with the simulation of such statistical models. Similarly, the restriction on tra c control in December, 2017 and the emission from reworks for celebrating Chinese Spring Festival also produce no signi cant effects on daily PM2.5 variations when compared with the effects of atmosphere circulation conditions. In some cases, theoretically, burning reworks can aggravate air pollution, but in the Spring Festivals of 2015 and 2016, the observed air quality is very good while the statistical model simulated medium peaks.
Without considering the changes of pollutant emissions, the simulated winter PM2.5 with the statistical downscaling models driven by the large-scale circulations during 1979-2017 has a slight but statistically signi cant decreasing trend, which is contrary to the observed increasing trend or that felt by the public in recent years due to the increased pollutant emissions. Similar to the observed daily PM2.5, the simulated daily PM2.5 values in 1979-2017 have a weak correlation to pan evaporation.
Overall, the CCs obtained by the statistical downscaling models are not ideally large, but it de nitely demonstrates that the SD models can re ect the important impacts of the large-scale circulation on the PM2.5. It is promising that there is some space to improve the model performance by combining more predictors of local meteorological factors, local emissions or pollutant transportation with these largescale predictors. Meanwhile, in this study, not all the useful signals in the high-correlation areas were used and only the signal at some of the high-correlation centers were used, therefore, the model performance perhaps may be improved by utilizing more signals, which need more e cient data-mining methods. Moreover, such SD models can be used to derive future projections of PM2.5 for local sites, based on the rough-resolution outputs of global climate models. Although only the PM2.5 at a site of Jiaozuo is used to construct the downscaling models, the high correlation areas of each large-scale variable to different sites in surrounding areas may be similar, since the PM2.5 concentration values among the neighboring sites are highly correlated to each other.

Declarations
Declaration of Competing Interest   Simulated annual average of daily Mean (a) and daily maximum (b), and the relation between daily PM2.5 and the observed pan evaporation: (c) the simulated daily mean PM2.5 versus pan evaporation; (d) the observed daily mean PM2.5 versus pan evaporation.