Impacts of meteorological uncertainties on the haze formation in Beijing–Tianjin–Hebei (BTH) during wintertime: a case study

In the present study, a persistent heavy haze episode from 13 to 20 January 2014 in Beijing–Tianjin– Hebei (BTH) is simulated using the WRF-CHEM model through ensemble simulations to investigate impacts of meteorological uncertainties on the haze formation. Model results show that uncertainties in meteorological conditions substantially influence the aerosol constituent simulations at an observation site in Beijing, and the ratio of the ensemble spread to the ensemble mean (RESM) exceeds 50 %. The ensemble mean generally preforms well in reproducing the fine particles’ (PM2.5) temporal variations and spatial distributions against measurements in BTH. The meteorological uncertainties do not alter the PM2.5 distribution pattern in BTH principally or dominate the haze formation and development, but remarkably affect the simulated PM2.5 level, and the RESM for the simulated PM2.5 concentrations can be up to 30 % at the regional scale. In addition, the rather large RESM in PM2.5 simulations at the city scale also causes difficulties in evaluation of the control strategies. Therefore, our results suggest that the ensemble simulation is imperative to take into account the impact of the meteorological uncertainties on the haze prediction.


Introduction
Over the past 3 decades, rapid industrialization and urbanization have caused severe air pollution in China, particularly during wintertime heavy haze, with extremely high levels of fine particles (PM 2.5 ) frequently engulfing the north of China (e.g., Chan and Yao, 2008;Fang et al., 2009;Zhao et al., 2013;Huang et al., 2014;Guo et al., 2014;Wu et al., 2017;Li et al., 2017a).Elevated atmospheric aerosols or PM 2.5 not only influence the Earth's climate system, but also remarkably impair visibility and potentially cause severe health defects (e.g., Penner et al., 2001;Pope and Dockery, 2006;Zhang et al., 2007).
Meteorological condition is critical for understanding the formation, transformation, diffusion, transport, and removal of the pollutants in the atmosphere.Dabberdt et al. (2004) have listed the meteorological research needs for improving air quality forecasting, one of which is to provide the model uncertainty information through ensemble prediction capabilities and quantify uncertainties and feedbacks between meteorological and air quality modeling components.Numerous studies have been performed in China to explore the role of meteorological conditions in the air pollution formation (e.g., Gao et al., 2011;Zhang et al., 2012Zhang et al., , 2015;;Wu et al., 2013;Wang et al., 2014;Bei et al., 2016a, b).Most recently, Liu et al. (2017) have investigated the meteorological impacts on the PM 2.5 concentrations over Beijing-Tianjin-Hebei (BTH) in December 2015.Their results have demon-strated that the unfavorable meteorological conditions are the main reason for deterioration of the air quality in BTH, while the undertaken emission control measures have only mitigated the air pollution slightly.
Previous studies on the air quality forecasting sensitivity to meteorological uncertainties mainly include Monte Carlo simulations (e.g., Dabberdt and Miller, 2000;Beekmann and Derognat, 2003) and adjoint sensitivity studies (e.g., Menut, 2003).The ensemble approach has also been applied to photochemical and secondary organic aerosol (SOA) simulations in various numerical models (e.g., Galmarini et al., 2004;McKeen et al., 2005), photo-chemical reactions (e.g., Delle Monache and Stull, 2003), emission scenarios (e.g., Delle Monache et al., 2006), physical parameterizations (e.g., Mallet and Sportisse, 2006), and meteorological initial conditions (e.g., Zhang et al., 2007;Bei et al., 2012).The ensemble means have generally performed better than most of the individual models.Uncertainties in meteorological initial conditions have been shown to substantially influence both ozone (O 3 ) and SOA simulations, including the peak time concentrations, the horizontal distributions, and the temporal variations (Zhang et al., 2007;Bei et al., 2012).Recently, Sharma et al. (2016) have evaluated uncertainties in surface O 3 simulations over the South Asian region during the premonsoon season due to different emission inventories and different chemical mechanisms.They have suggested that the assessment of the tropospheric O 3 budget and its implications for public health and agricultural output should be conducted prudently considering the huge uncertainties caused by emission inventories and chemical mechanisms.Solazzo et al. (2017) have emphasized the high interdependencies among meteorological and chemical variables and the related errors, indicating that the evaluation of the air quality model performance needs to be confirmed by more complementary analysis of meteorological fields and chemical precursors.
The purpose of the present study is to explore impacts of the uncertainties in meteorological conditions on the PM 2.5 simulations or forecasts in BTH through ensemble simulations using the WRF-CHEM model.The methodology and model are presented in Sect. 2. The analyses, results, and discussions are included in Sect.3. The summary and conclusions are given in Sect. 4.

WRF-CHEM model
A specific version of the WRF-CHEM model is used to examine impacts of the uncertainties in meteorological conditions on the PM 2.5 simulations or the haze formation in BTH, which is developed by Li et al. (2010Li et al. ( , 2011aLi et al. ( , b, 2012) ) (Binkowski and Roselle, 2003).The inorganic aerosols are predicted using ISORROPIA version 1.7 (Nenes et al., 1998).The SOA formation is simulated using a non-traditional SOA module, including the volatility basis set (VBS) modeling method and the SOA contributions from glyoxal and methylglyoxal.A detailed description of the WRF-CHEM model can be found in Li et al. (2010Li et al. ( , 2011aLi et al. ( , b, 2012)).A persistent heavy haze pollution episode from 13 to 20 January 2014 in BTH is simulated.The model simulation domain is shown in Fig. 1, and detailed model configurations can be found in Table 1.

Ensemble initialization method
The ensemble initialization method used in the present study is called the "climatological ensemble initialization method" (Zhang et al., 2007;Bei et al., 2012).In the approach, dynamically consistent initial and boundary conditions are statistically sampled from a seasonal meteorological data set.In order to represent the wintertime climatological statistics, a data set during the period from 1 November 2013 to 28 February 2014 is generated using NCEP-FNL 1 • × 1 • reanalysis data.The perturbed variables include the horizontal wind components, potential temperature, perturbation pressure, and mixing ratio of water vapor.Other prognostic variables such as vertical velocity and mixing ratios of hydrometeors are not perturbed.In general, the perturbation in horizontal wind components constitutes the most important uncertainty in those variables (Bei et al., 2008(Bei et al., , 2010)).Thirty  Guenther et al. (2006) ensemble members are randomly chosen from this climatological data set.Similarly, boundary conditions for each ensemble member are generated from the same data set beginning at the randomly selected initial time of the given member, and extended for the same length of time as the simulated episode.Deviations of the initial and boundary condition data for each member from the climatological mean for the entire period are then scaled down to 20 % to reduce the ensemble spread to less than typical observation error magnitudes (Nielsen-Gammon et al., 2007) and added to the unperturbed initial and boundary conditions derived directly from the NCEP-FNL analyses valid at 12:00 UTC on 12 January 2014, which are used for the 6 km domain ensemble simulation.Figures 2a-d show the vertical distribution of the average initial ensemble spread which is calculated as the standard deviation of the perturbations imposed on each ensemble member's initial field.The average spread is 0.5-3.0m s −1 for horizontal winds (U and V components), 0.5-1.1 K for temperature, 0.02-0.48hPa for pressure, and 0-0.15 g kg −1 for the water vapor mass mixing ratio.The initial ensemble spreads of meteorological variables are generally less than their typical observation error magnitudes.It is worth noting that all the ensemble simulations used the same initial and boundary conditions for chemical fields, as well as the same anthropogenic emission inventory.

Pollutant measurements
The hourly near-surface CO, SO 2 , NO 2 , O 3 , and PM 2.5 mass concentrations in BTH are released by the China Ministry of Environmental Protection (China MEP) and can be downloaded from the website at http://www.aqistudy.cn/.The Aerodyne High Resolution Time-of-Flight Aerosol Mass Spectrometer (HR-ToF-AMS) with a novel PM 2.5 lens is used to measure the sulfate, nitrate, ammonium, and or-ganic aerosols (OA) from 9 to 26 January 2014 at the Institute of Remote Sensing and Digital Earth (IRSDE), Chinese Academy of Sciences (40.00 • N, 116.38 • E) in Beijing (Fig. 1) (Williams et al., 2013).The positive matrix factorization (PMF) technique is utilized with constraints implemented in SoFi (Canonaco et al., 2013) to analyze the sources of OA and five components are separated by their mass spectra and time series.The components include hydrocarbonlike OA (HOA), cooking OA (COA), biomass burning OA (BBOA), coal combustion OA (CCOA), and oxygenated OA (OOA).HOA, COA, BBOA, and CCOA are interpreted for surrogates of primary OA (POA), and OOA is a surrogate for SOA.Detailed information about the HR-ToF-AMS measurements and data analysis can be found in Elser et al. (2016).A lidar has also been deployed at IRSDE and the aerosol backscatter signal is used to retrieve the planetary boundary layer (PBL) height.
3 Results and discussions

Synoptic overview
Figure 3 shows temporal evolutions of the observed PM 2.5 mass concentrations averaged over 13 cities (see Fig. 1) in BTH during the severe haze episode from 13 to 21 January 2014.The observed PM 2.5 mass concentrations are frequently higher than 250 µg m −3 in the 13 cities during the episode, exceeding the standard of severe pollution (hourly PM 2.5 mass concentration exceeding 250 µg m −3 , Feng et al.,  2016).The haze in BTH was in the stage of development from 13 to 15 January, with the gradual increase in the PM 2.5 concentration.BTH was most polluted when the haze was in the maturity stage on 16 January, with the PM 2.5 concentration exceeding 400 µg m −3 in most of the cities.From 17 to 19 January, the PM 2.5 concentrations fluctuated considerably, which was primarily caused by the transition between different synoptic situations.During nighttime on 19 January, the haze in BTH rapidly dissipated, with the PM 2.5 concentration decrease of several hundreds of µg m −3 in 2 or 3 h.In addition, the diurnal cycles of the observed PM 2.5 mass concentrations were not clear, demonstrating the obvious regional pollution characteristics in BTH.For the four megacities in BTH, the PM 2.5 levels in Shijiazhuang and Baoding were much higher than Beijing and Tianjin, which is caused by the massive local emissions in Shijiazhuang and Baoding.NCEP-FNL reanalysis data are used to examine the effect of synoptic conditions on the air pollution during the haze episode in BTH.Figures S1-S3 in the Supplement show the synoptic conditions at the surface level, 850 and 500 hPa, respectively.On 13 January, BTH is to the north of a high pressure at the surface level, causing the southerly wind in the east of BTH, and sandwiched between the trough in the northeast of BTH and the high pressure in the southwest of BTH at 850 hPa, inducing the westerly surface wind in the west of BTH.At 500 hPa, BTH is situated in the rear of the trough, and the westerly airflow is dominant.The air pollutants in BTH are subject to transport to the east but hindered by the southerly wind, causing accumulation of air pollutants.On 14 January, the high-pressure system begins to control BTH at the surface level and 850 hPa, and the wind is varied and weak, favorable for the accumulation of air pollutants in BTH.On 15 January, BTH is still controlled by the high pressure at the surface level and 850 hPa, and the westerly wind prevails at 500 hPa.The weak surface wind, together with the stable stratification, further facilitates accumulation of air pollutants in BTH.On 16 January, a trough develops over BTH at 850 and 500 hPa, and BTH is situated near the trough line, in which the northerly and southerly winds occur at the same time.At the surface level, the northerly wind prevails in the north of BTH and the southerly wind prevails in the south of BTH, leading to evacuation of air pollutants in the north of BTH and the high level of air pollutants in the south of BTH.On 17 January, the trough at 850 hPa commences weakening and the controlling region of the trough at 500 hPa becomes narrow.The northwesterly wind is dominant over BTH, leading to divergence of the air pollutants in BTH.On 18 January, BTH is located near the ridgeline at 850 hPa and at the verge of the high pressure at the surface level.The controlling scope of the high-pressure system at the surface level is wide, inducing the varied wind over BTH, and is not conducive to the evacuation of air pollutants in BTH.On 19 January, the prevailing southerly wind in the south of BTH and the strong westerly wind in the west of BTH lead to the convergence of air pollutants at the surface level.At 850 and 500 hPa, BTH is situated in the southeast of the trough and southwesterly wind is prevalent.On 20 January, BTH is located in the southwest of the trough at 500 and 850 hPa, and the strong northwesterly wind prevails over BTH.At the surface level, BTH is situated between the high pressure in the west and the low pressure in the east, inducing the strong northwesterly wind over BTH.The cold clean air sweeps BTH and efficiently decreases the air pollutant concentrations in BTH.

Uncertainties in meteorological simulations
Figures 4a-d provide the temporal profiles of the ensemble simulations of the surface meteorological fields and the corresponding observations at the meteorological site in Beijing from 13 to 20 January 2014.The U component exhibits larger ensemble spread than the V component (Fig. S4), but the ensemble mean (ENSM) of the U component generally yields the observed diurnal variations.The ensemble prediction of the V component fails to reproduce the observed intensified southerly or northerly winds.The meteorological site is located in the north of the Yanshan Mountains, substantially influenced by the mountain-valley circulation (MVC).Apparently, the WRF-CHEM model lacks the ability to simulate the occurrence and development of MVC well, causing the considerable biases of the ensemble prediction of the V component.The ensemble prediction performs well in producing the diurnal variation of the surface temperature, but the underestimation or overestimation is still large when the V component prediction is biased.The relative humidity (RH) shows a rather large ensemble spread (Fig. S4d), and the ENSM reasonably tracks the observed diurnal variation, with high nighttime and low afternoon simulated RH.The RH simulation is sensitive to the simulated surface temperature.Generally, the overestimation of the surface temperature corresponds well to the underestimation of the RH, or vice versa.The ENSM considerably overestimates the PBL height during daytime on 13 and 14 January, and underestimates it on 15 January (Fig. 4e).In addition, most of the ensemble members frequently underestimate the observed PBL height during nighttime, and all ensemble members fail to produce the peak PBL height on 17 and 20 January.The PBL height is principally determined by the vertical shear of horizontal winds and the ground thermal condition.Therefore, uncertainties of wind and temperature field simulations cause large biases of the PBL height simulation.

Uncertainties in aerosol species simulations
Figure 5 shows the temporal profiles of the ensemble simulations of the aerosol species and the observations at IRSDE in Beijing.The ENSM reasonably produces the observed variations of the POA concentrations.However, all ensemble members fail to capture the peaks in the morning on 16 January and in the evening on 17 January, indicating that the underestimation might not be caused by the meteorological uncertainties, but by emission biases.The POA in the atmosphere is from multiple sources, including the direct emissions from vehicles, cooking, biomass, and coal combustion.Diurnal variations of those sources might constitute one of the major reasons for the biases of the POA simulations.The ENSM generally performs reasonably well in simulating the SOA concentration against the measured OOA.The ratio of the ensemble spread to the ensemble mean (RESM) for the SOA prediction is large compared to that of POA (Fig. S5a, b).Four SOA formation pathways are included in simulations: oxidations of anthropogenic and biogenic volatile organic compounds (VOCs), oxidation HOA semi-volatile vapors, and irreversible uptake of glyoxal and methylglyoxal on aerosol surfaces.Therefore, uncertainties in meteorological fields influence not only the transport of the SOA precursors, but also the SOA formation processes in the atmosphere, causing the rather large RESM of SOA simulations.The ENSM generally reproduces the observed variations of sulfate, nitrate, and ammonium (SNA), but the RESM of SNA is also considerably large (Fig. S5c-d).During haze days, sulfate is primarily formed through heterogeneous reactions of SO 2 on aerosol surfaces, which is highly dependent on the relative humidity (Li et al., 2017b).Nitrate formation is determined by the HNO 3 and N 2 O 5 that originated from the NO 2 oxidation, is sensitive to the temperature and relative humidity, and is also influenced by the level of sulfate in the particle phase and ammonia in the atmosphere.The ammonium aerosol is formed through neutralization of sulfate and nitrate aerosols by NH 3 .Additionally, in the present study, ISORROPIA (version 1.7) is used to calculate the thermodynamic equilibrium between the sulfatenitrate-ammonium-water aerosols and their gas-phase precursors H 2 SO 4 -HNO 3 -NH 3 -water vapor.Therefore, uncertainties in meteorological fields propagate to the transport, atmospheric oxidation, and thermal dynamic processes, which all have contributions to the large RESM of the SNA simulations.Apparently, uncertainties in meteorological conditions substantially affect the aerosol species simulations at a single observation site, which is consistent with the previous studies (Bei et al., 2012).

Uncertainties in PM 2.5 simulations in BTH
Heavy haze with high levels of PM 2.5 frequently constitutes a regional pollution event.Figure 6 shows the temporal profiles of the ensemble simulations and observations of air pollutants averaged at the monitoring sites in BTH from 13 to 20 January 2014.The RESM of the average air pollutants is much less than those of aerosol species at the single observation site (Fig. S6).For the primary air pollutants, SO 2 and CO, the ENSM generally tracks reasonably the observed variations.However, sometimes all the ensemble members underestimate or overestimate the observation.There are two possible reasons for the biases of ensemble simulations of SO 2 and CO: uncertainties of emissions and systematic errors of meteorological fields.In the evening on 15 January, the ensemble prediction substantially overestimates the observed SO 2 concentration, but CO overestimation is not large.In contrast, in the morning on 16 January, the ensemble prediction slightly underestimates the SO 2 observation but noticeably underestimates the CO concentration.Therefore, the overestimation of SO 2 on 15 January and underestimation of CO on 16 January might be primarily contributed by the emission uncertainties.In the morning on 18 January, the ensemble prediction significantly underestimates both SO 2 and CO observations, indicating the plausible uncertainties caused by the systematic errors of meteorological fields.
The ENSM of the average surface O 3 and NO 2 over the monitoring sites in BTH is in good agreement with observations.The ensemble prediction is prone to underestimating the O 3 observation during nighttime, but is very consistent with the NO 2 observation.Considering the massive NO x emission and the titration of NO, the nighttime O 3 concentrations are generally very low, particularly during wintertime when the daytime O 3 concentrations are not high.Hence, the underestimation of nighttime O 3 concentrations is perhaps caused by the observation uncertainties, such as the setting of a lower detection limit.In addition, the ENSM does not reproduce the high O 3 level during nighttime on 19 January when the northwesterly wind is intensified to evacuate the air pollutants in BTH.Rapid increase in the observed O 3 concentrations during nighttime shows the substantial contribution of the background O 3 transport.Therefore, the background O 3 uncertainties constitute the major reason for the O 3 underestimation on 19 January.
The ENSM also performs well in replicating the observed PM 2.5 observation, except for the underestimation on 16 and 18 January.However, the RESM of the PM 2.5 simulations is larger than those of O 3 , NO 2 , SO 2 , and CO (Fig. S6).The average ENSM of the PM 2.5 concentration over the monitoring sites during the simulation period is 189.5 µg m −3 , close to the observed 197.6 µg m −3 .In addition, the ensemble member of 16 and 30 respectively) produces the highest and lowest PM 2.5 levels, with average PM 2.5 concentrations of 231.5 and 167.3 µg m −3 , respectively.The PM 2.5 mainly includes the primary aerosols which are determined by direct emissions, and the secondary aerosols which are determined by their precursor emissions and the homogeneous and heterogeneous oxidation process in the atmosphere.Therefore, the large RESM of SOA and SNA simulations enhances the ensemble spread of the PM 2.5 simulations.
Figure 7 presents the spatial distributions of ENSM and observations of the daily average near-surface PM 2.5 mass concentrations during the haze episode, along with the simulated wind fields.The ENSM predicted PM 2.5 spatial patterns are generally in good agreement with the observations at the ambient monitoring sites in BTH.The ENSM successfully reproduces the haze development and maturity stages from 13 to 16 January 2014.From 17 to 18 January, the northeasterly wind develops and decreases the PM 2.5 level in BTH, but not strongly enough to evacuate the air pollutants.The PM 2.5 pattern of ENSM is very consistent with observations, but on 18 January, the PM 2.5 concentrations are remarkably underestimated in four cities in BTH.On 19 January, the westerly wind prevails in BTH, causing the divergence of the PM 2.5 .On 20 January, the intensified northwesterly wind begins to empty the PM 2.5 in BTH.However, apparently, the occurrence of the intensification of the northwesterly wind is early, causing considerable underestimation of the PM 2.5 concentration in the ENSM.
The uncertainties in meteorological fields are generally less than observational and analysis errors, but the ensemble simulations still exhibit considerable spreads.In order to contrast the PM 2.5 simulations of different ensemble members, we have selected two members: EN-16 and EN-30, representing the highest and lowest PM 2.5 simulations in BTH, respectively.Figures 8 and 9 provide the horizontal distributions of the daily average surface PM 2.5 concentrations along with surface winds during the episode in respectively. Similar PM 2.5 distribution patterns are simulated in  showing that the meteorological uncertainties do not dominate the haze formation and development principally.The PM 2.5 level in EN-16 is much higher than that in EN-30 in BTH, which is mainly caused by the considerable discrepancies in the surface winds between the two members.The simulated southerly wind in EN-16 is generally more intense than that in EN-30, but the northerly wind in EN-16 is weak compared to EN-30, which is more favorable for the air pollutant accumulation in EN-16 than in EN-30.On 13 and 14 January, the winds in EN-30 are weak or calm in BTH and the PM 2.5 is mainly attributed to the local production.However, in EN-16, the prevailing southerly winds also deliver the air pollutants from the southern areas to BTH, substantially enhancing the PM 2.5 level.On 15 January, although EN-16 and EN-30 both produce the prevailing southerly wind in BTH, the westerly wind in EN-30 is intense compared to EN-16, considerably decreasing the PM 2.5 level in EN-30.On 16 January, the northeasterly wind in EN-30 is intensified and evacuates the PM 2.5 in the north of BTH.However, in EN-16, the simulated northeasterly wind is weak and the PM 2.5 level in the north of BTH still remains high.On 17 January, the simulated northerly wind in EN-16 is weak compared to that in EN-30, causing a higher PM 2.5 concentration in EN-16 than EN-30 in BTH.On 18 January, the intensified southerly wind in EN-16 considerably increases the PM 2.5 level in BTH compared to EN- 30.On 19 January, the westerly wind is prevalent in EN-30 and the PM 2.5 level begins to decrease, but in EN-16, the southwesterly wind still causes high PM 2.5 concentrations in BTH.On 20 January, the stronger northeasterly wind in EN-30 more efficiently evacuates the PM 2.5 than that in EN-16.
3.5 Uncertainties in PM 2.5 simulations in mega-cities EN-16 and EN-30 both predict the haze occurrence and development in BTH during the episode, although the difference in the PM 2.5 level between those two members is considerable, showing that the meteorological uncertainties do not dominate the regional haze formation.Previous studies have shown that the meteorological uncertainties substantially impact the air quality simulations at the city scale (Bei et al., 2012).Figure 10 presents the temporal variation of the ensemble simulations and observations averaged at four mega-cities in BTH during the episode.The ENSM of the PM 2.5 concentrations in Beijing, Tianjin, and Baoding is in good agreement with the observation.However, the ENSM remarkably underestimates the observed PM 2.5 concentration in Shijiazhuang from 16 to 19 January, which is hardly interpreted by the emission biases.The ENSM performs well in simulating the PM 2.5 variations from 13 to 15 January, and overestimates the observation on 20 January in Shijiazhuang.One of the possible reasons for the underestimation in Shijiazhuang is that the westerly wind is systematically overestimated from 16 to 19 January along the foothills of the Taihang Mountains, causing the haze plume to move eastward.
Although the ENSM produces reasonably well the PM 2.5 variations in the four mega-cities against the measurement, the meteorological uncertainties still cause large uncertainties in the PM 2.5 concentration (Fig. S7).During the first 3 days of the episode, the ENSM is very consistent with the observations in the four mega-cities, but the PM 2.5 level discrepancy between the members with the highest and lowest PM 2.5 concentrations is rather large, causing troubles for the assessment of the control strategies.For example, in Shijiazhuang, the average PM 2.5 concentrations during the first 3 days in the members with the highest and lowest PM 2.5 concentrations are 403.5 and 213.8 µg m −3 , respectively, and the difference is about 190 µg m −3 .In Beijing, the average PM 2.5 concentrations in the two members are 103.9 and 196.3 µg m −3 .It is worth noting that, according to the Chinese air quality standard released in 2012, the PM 2.5 concentration of 103.9 µg m −3 is defined as a "lightly polluted condition", but that 196.3 µg m −3 is defined as a "heavily polluted condition".If the heavy air pollution occurs, the control strategies will be implemented.Therefore, it is necessary to use the ensemble simulation to avoid the impact of the meteorological uncertainties the haze prediction.

Summary and conclusions
In the present study, the uncertainties in simulating haze formation due to meteorological uncertainties are investigated using the WRF-CHEM model through ensemble simulations.A persistent heavy haze episode that occurred in BTH from 13 to 20 January 2014 is simulated.A climatological ensemble initialization approach is used to produce initial and boundary conditions for each ensemble member.
The ENSMs of the aerosol constituents are generally in good agreement with the observations at an observation site in Beijing, including the sharp buildup of the aerosol constituents in the evening on 15 January and rapid falloff in the morning on 20 January.However, the ENSM considerably underestimates the observed primary aerosols in the evening on 17 January.The ensemble spread is rather large for the aerosol constituent simulations, and the RESM exceeds 50 %, respectively.
The ENSM performs well in simulating the temporal variations of the average surface CO, SO 2 , NO 2 , O 3 , and PM 2.5 mass concentrations over the monitoring sites in BTH, and the RESM of the air pollutants is generally less than 30 %.The RESM of PM 2.5 simulations is larger than the other air pollutants, which is due to the complicated composition of PM 2.5 , including the contributions of primary and secondary aerosols.The meteorological uncertainties do not principally dominate the haze formation and development, but considerably alter the simulated PM 2.5 level.The average PM 2.5 difference during the episode exceeds 60 µg m −3 between the two members with the highest and lowest PM 2.5 simulations.
Although the meteorological uncertainties do not dominate the regional haze formation, they still substantially influence the PM 2.5 simulations at city scale.The ENSM predicts the PM 2.5 variations in the four mega-cities against the measurements reasonably well, including Beijing, Tianjin, Baoding, and Shijiazhuang, but the RESM of the PM 2.5 simulations is rather large, causing troubles for the evaluation of the control strategies.Therefore, the ensemble simulation is needed to take into consideration the impact of the meteorological uncertainties on the haze prediction.It is worth noting that aside from meteorological fields, uncertainties in emissions or various chemistry/aerosol schemes also considerably influence the WRF-CHEM model simulations.The extended response surface modeling (ERSM) technique can be used to quantify the relative importance of each uncertainty source in the WRF-CHEM model (Zhao et al., 2017).
Data availability.Data used in the present study can be provided by Guohui Li (ligh@ieecas.cn).
Competing interests.The authors declare that they have no conflict of interest.Special issue statement.This article is part of the special issue "Regional transport and transformation of air pollution in eastern China".It is not associated with a conference.
at the Molina Center for Energy and the Environment.The model includes a new flexible gas-phase chemical module and the CMAQ/Models-3 aerosol module developed by the US EPA

Figure 1 .
Figure 1.WRF-CHEM simulation domain.The filled red (in BTH) and blue (outside of BTH) circles represent centers of cities with ambient monitoring sites.The size of the circle denotes the number of ambient monitoring sites of cities.The filled black triangle and rectangle denote the deployment location of the HR-ToF-AMS and the surface meteorological site in Beijing, respectively.

Figure 2 .
Figure 2. Vertical distribution of the mean of initial ensemble spreads and the standard deviation for (a) horizontal winds (U and V components), (b) temperature, (c) pressure, and (d) water vapor mixing ratio.

Figure 3 .
Figure 3. Observed hourly PM 2.5 concentrations averaged in (a) four mega-cities (Beijing, Tianjin, Baoding, and Shijiazhuang) and (b) nine non-mega-cities of BTH during the period from 13 to 20 January 2014.

Figure 4 .
Figure 4. Temporal evolution of the surface (a) U component, (b) V component, (c) temperature, and (d) relative humidity at the meteorological site, and (e) the PBL height at IRSDE in Beijing from each ensemble member (thin green lines), the ensemble mean (bold black line), and observations (black dots) from 13 to 20 January 2014.

Figure 5 .
Figure 5. Temporal evolution of the (a) POA, (b) SOA, (c) sulfate, (d) nitrate, and (e) ammonium mass concentrations at IRSDE in Beijing from each ensemble member (thin green lines), the ensemble mean (bold black line), and observations (black dots) from 13 to 20 January 2014.

Figure 6 .
Figure 6.Temporal evolution of the (a) PM 2.5 , (b) O 3 , (c) NO 2 , (d) SO 2 , and (e) CO mass concentrations averaged over monitoring sites in BTH from each ensemble member (thin green lines), the ensemble mean (bold black line), and observations (black dots) from 13 to 20 January 2014.

Figure 7 .
Figure 7. ENSM of the daily average surface PM 2.5 concentration distributions (color contour) along with the ENSM of the daily average surface winds (black arrows) from 13 to 20 January 2014.The colored circles denote the PM 2.5 measurements in cities.

Figure 8 .
Figure 8. Same as Fig. 7 but for the ensemble member of 16 with the highest simulated PM 2.5 concentration.

Figure 9 .
Figure 9. Same as Fig. 7 but for the ensemble member of 30 with the lowest simulated PM 2.5 concentration.

Figure 10 .
Figure 10.Temporal evolution of the PM 2.5 mass concentrations averaged in (a) Beijing, (b) Tianjin, (c) Baoding, and (d) Shijiazhuang from each ensemble member (thin green lines), the ensemble mean (bold black line), and observations (black dots) during the period from 13 to 20 January 2014.The red and blue lines represent the simulations in the members with the highest and lowest PM 2.5 concentrations, respectively.