Sensitivity of WRF-Chem model resolution in simulating particulate matter in South-East Asia

Frequent haze episodes commonly caused by biomass burning has been a concern in South East Asia. One of such 15 events was the June 2013 severe haze in the region. This study assessed the ability of WRF-Chem in capturing the spatial variability and concentrations of particulate emissions during this period. It analyzed the regional biomass burning emissions and its transport leading to higher particulate matter levels in the region. In order to analyze the effect of grid scale, the horizontal resolution of the simulation was varied between low-resolution (100 km) and high-resolution (20 km). Evaluations of the simulations were made against meteorological observations pertinent to emission and transport of particulate matter, 20 including surface and vertical air profile variables such as temperature, relative humidity, and wind speed and direction. Particulate matter (PM10 and PM2.5) levels were evaluated using ground measurements in Brunei and Singapore respectively. The meteorological parameters were adequately represented across the model simulations. Increasing the horizontal resolution of the simulations generally improved the representation of meteorology and air quality but some progronostic variables maintained similar or better performance with coarse resolution simulation. With the high-resolution simulation, PM10 25 concentration in Brunei had a correlation coefficient around 0.4, and the simulated PM2.5 level in Singapore had correlation coefficient around 0.9. Whereas, the low-resolution simulation had correlation coefficidents around 0.2 and 0.8 for PM10 and PM2.5 levels at Brunei and Singapore, respectively. Both simulations could not repeat aerosol optical depth (AOD) from reanalysis unless the biomass burning emissions were enhanced. An enhacement factor of 6 with high resolution simulation gave PM10 and PM2.5 correlations around 0.6 and 0.9 in Brunei and Singapore respectively. 30 https://doi.org/10.5194/acp-2019-692 Preprint. Discussion started: 23 March 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Biomass burning is very predominant in the South-East Asian region due to agricultural activities including frequent slash and 35 clearing of the residual crop in preparation for the new planting season (Oozeer et al., 2016). These activities lead to higher emission of different air pollutants and particulate matter across the region. Regional biomass burning could be observed using remote sensing such as two-dimensional images and thermal anomalies captured by satellites. The June 2013 haze in South East Asia was evident from satellite imageries (NASA, 2013). During this period, air quality data from the Department of Previous studies of the June 2013 haze episode mostly included statistical studies of the effect of particulate matter on visibility and air quality (Betha et al., 2014;Bhardwaj et al., 2016;Dotse et al., 2016Dotse et al., , 2018Gaveau et al., 2015;Lee et al., 2017;Lin et 45 al., 2014). Oozeer et al. (2016) numerically studied and analyzed convective mechanisms responsible for the uplift and transport of the particulate emissions from Sumatra over to Peninsular Malaysia during the 2013 event (Oozeer et al., 2016).

Model Chemistry and Physics
In all simulations, the 2-moment Morrison scheme (list references) was used to represent the resolved-scale microphysics. The radiative processes (both longwave and shortwave) were represented by the Rapid Radiative Transfer Model for general circulation models (GCMs) longwave and shortwave radiation schemes (RRTMG) (Mlawer et al., 1997;Pincus et al., 2003) respectively. These two schemes can interact with WRF-Chem to include effects of aerosols on radiation through absorbing 95 and scattering. Convective parameterization was performed by the Grell 3D cumulus scheme, which was an improved version of the Grell-Devenyi (GD) cumulus scheme (Grell and Dévényi, 2002). In addition, the Grell 3D cumulus scheme allows for feedback from parameterized convection to the radiation schemes. Mellor-Yamada-Janjić (MYJ) boundary layer scheme (Hong and Pan, 1996) was adopted, and the moisture and thermal fluxes from the surface was obtained using the Noah landsurface model (Chen and Dudhia, 2001). We selected the Regional Acid Deposition Model, 2nd generation (RADM2) 100 (Stockwell et al., 1990) gas-phase chemistry mechanism coupled with the Modal Aerosol Dynamics Model for Europe (MADE), which included the Secondary Organic Aerosol Model (SORGAM) (Schell et al., 2001), for the simulation.
The MADE/SORGAM module is a modal scheme that describes three log-normally distributed modes to simulate particle size distribution: the Aitken mode (<0.1 m diameter), the accumulation mode (0.1-2 m diameter), and the coarse mode (>2 m diameter), and predicts mass and number concentrations for each aerosol mode. In each mode, particles are assumed to have 105 the same chemical composition (internally mixed), while they are externally mixed among different modes (Zhao et al., 2010).
The aerosol direct effect was turned on by activating aerosol radiative feedback in the simulation.

Data for Model Set-up and Emission Input 110
This study used meteorological initial and boundary condition from the National Center for Environmental Prediction FiNal reanalysis (NCEP-FNL) for the initialization of meteorological variables (NOAA/NCEP, 2000). The NCEP-FNL has a horizontal resolution of 1-degree by 1-degree prepared operationally every six hours. The dataset is available on the surface and at 26 mandatory (and other pressure) levels from 1000 hPa to 10 hPa. Parameters include surface pressure, sea level pressure, geopotential height, temperature, sea surface temperature, soil values, ice cover, relative humidity, u-and v-winds, 115 vertical motion, vorticity, and ozone (NOAA/NCEP, 2000).
The trace gas and aerosol emissions fields were preprocessed using an adapted version of PREP-CHEM-SRC (Freitas et al., 2011). PREP-CHEM-SRC integrates emissions inventories and forest fire hotspots onto a specified domain. WRF-Chem needs three emissions input inventories to simulate the air pollution: (1) anthropogenic emissions; (2) biomass burning emissions and; (3) biogenic emissions. The biomass burning emissions were estimated using the Brazilian Biomass Burning Emissions 120 https://doi.org/10.5194/acp-2019-692 Preprint. Discussion started: 23 March 2020 c Author(s) 2020. CC BY 4.0 License. Model (3BEM); , which were injected into the atmosphere using the plume rise model described by Freitas and co-workers ( (Freitas et al., 2006(Freitas et al., , 2007. Daily fires detected from different satellites could be utilized and processed by the 3BEM. The forest fire hotspots database integrated into the model typically includes a combination of the MODIS C6 fire location database and the Wildfire Automated Biomass Burning Algorithm (WFABBA) database (University of Wisconsin-Madison, 2013). Aqua and Terra products from MODIS with 1000m resolution (at nadir) have been adopted in 125 many cases but small fires that are common in the Borneo and South-East Asia region are less likely to be captured in this product. We have leveraged on the combination of similar product with higher resolution -Visible Infrared Imaging Radiometer Suite (VIIRS) carried on the Suomi-NPP satellite with a 375m resolution (at nadir), to augment and enhance the fire detection database fed into the 3BEM. The 3BEM filters all the fire products and removes fires within a kilometer of each other to avoid double counting between fire datasets. The anthropogenic emissions were obtained from two global inventory 130 datasets -REanalysis of TROpospheric chemical composition over the past 40 years (RETRO) and Emissions Database for Global Atmospheric Research (EDGAR).
The biogenic emissions data is provided by the Model of Emissions of Gases and Aerosols from Nature version 2.1 (MEGAN2.1). MEGAN2.1 is a model framework for a resolution of 1 km. This model estimates natural emissions that occur in a terrestrial ecosystem (Guenther et al., 2012). 135 Chemical initial and boundary conditions are needed to account for initial concentrations and inflow/background concentrations. WRF-Chem uses idealized chemical profile generated from the NALROM simulation (based upon northern hemispheric, mid-latitude, clean environment conditions). In this study, the initial chemical boundary condition was set in WRF-Chem using 'mozbc' tool (MOZART GCM output, reference).
To ensure that the model did not violate the Courant-Friedrichs-Levy (CFL) stability criterion (Courant et al., 1928) the time 140 step of the simulation was set at 600 and 120 seconds for the low-resolution (100km) and high-resolution (20km) simulations respectivelyfollowing the recommended maximum of 6 × grid spacing. The WRF-Chem simulations were initialized at 00UTC on 15 June 2013 and ended at 00UTC on 30 June 2013. The emission inventories were updated for each day of simulation and the model outputs were generated every hour. The simulations were restarted at the beginning of each day with initial values from the previous run. 145

Air Quality and Meteorological Data for Validation
Chemical transport models such as WRF-Chem serve as a tool to quantify and establish the deterministic relationship between the amount of emissions and measured concentrations. The performance of the model must be measured and quantified in order to establish the confidence of such relationship (Kumar et al., 2015). Thus, ground measurements were used to evaluate the model performance of the simulations in this study. Ground observations for PM10 were obtained from Department of  Radiosonde data provided by the Department of Atmospheric Science at the University of Wyoming were obtained for Brunei Airport weather station (4.93° N, 114.93° E) and Singapore station (1.36º N, 103.98º E) to evaluate the WRF-Chem simulated atmospheric temperatures (º C), relative humidity (%), wind speed (m/s) and wind direction (degree) (University of Wyoming, 2013). Based on the geographical location of Singapore and Brunei Darussalam in South East Asia, evaluating the simulated convective mechanism in both region is very suitable in the assessment of the June 2013 haze event in the region. Data 160 availability can be seen in Table 1.

Grid Staggering and Meteorological Evaluation
Aside considering model resolutions, performance of model simulation is influenced by the arrangements of prognostic 165 variables on staggered grids relative to observation. It is also influenced by the interpolation method (e.g. nearest neighbor, bilinear interpolation) adopted in calculating variables from simulations.
In WRF and WRF-Chem model, grid cells are in rectangular distribution and according to Arakawa and Lamb C grid staggering (see Figure 1 top left) introduced by Arakawa and Lamb (Arakawa and Lamb, 1977). The height (h) is evaluated at the intersection of grid cells while components of wind variables (u, v) are evaluated at mid-points between grid cells. Scalar 170 variables such as thermodynamic variables (e.g. pressure and temperature) are evaluated at the theta points -not staggered and at the grid center (Collins et al., 2013). Ground measurements at the Brunei Airport station were evaluated against surface parameters from simulations such as the temperature at 2 m, surface relative humidity, wind speed and direction at 10 m. Also, the vertical air profile from the simulations were evaluated against the sounding (radiosonde) data for Singapore and Brunei Airport provided by the 190 Department of Atmospheric Science at the University of Wyoming (University of Wyoming, 2013). Vertical profiles of temperature (º C), relative humidity (%), wind speed (m/s and wind direction (degree) were evaluated. The wind speed (m/s), wind direction (degree), relative humidity (%), atmospheric temperature (º C), pressure (hPa) and the geopotential height (m) corresponding to Singapore station (1.36º N, 103.98º E) and Brunei Airport weather station (4.93° N, 114.93° E) were extracted between 15th to 30th June 2013, and compared to the sounding data. The statistical metrics used to evaluate the model 195 performance are the Pearson correlation coefficient (r) (Eq. 1), root mean square error (RMSE) (Eq. 2), normalized root mean square error (NRMSE) (Eq. 3) and normalized mean bias factor (NMBF) (Eqs. 4a. and 4b.). The metrics are defined as follows,   Figure 2 shows the time series comparison of temperature at 2 m (º C) of the observed hourly temperature to the WRF-Chem simulations at Brunei station. In Brunei station, surface temperature simulated by both high and low horizontal resolution simulation had very good correlation (above 0.9) with the observation. However, it is noticeable from 215 the comparison of time-avergaed surface temperature distrubtion across the region between reanalysis and simulations that higher resolution simulations resolved terrain with temperature gradient better at regions in Sumatra and Borneo with higher altitudes (Figure 3). Although WRF-Chem_20km had lower error and bias in Brunei station using nearest neighbor calculation method with RMSE, NRMSE and NMBF of 1.760 º C, 0.154 and -0.039 respectively, WRF-Chem_100km simulation had slightly better correlation (0.924 against 0.908). Also, the bilinear interpolation calculation yielded better results in WRF-220 Chem_100km than WRF-Chem_20km for surface temperature in Brunei station (see Table 2).
Considering air temperature in Brunei, both low and high-resolution simulations maintained similar and very good results (correlation coefficient of 0.999) with barely any difference between nearest neighbor and bilinear interpolation estimation methods. Overall, the bilinear interpolation in the high-resolution simulation gave the least error (RMSE and NRMSE of 1.759 º C and 0.016 respectively), while nearest neighbor estimation gave the least bias value of -0.011 (Table 3). 225 In Singapore, the air temperature simulation results were similar to Brunei. Both low and high-resolution simulations had the same correlation coefficient of 0.999 with nearest neighbor and bilinear interpolation techniques. While low-resolution simulation estimated using bilinear interpolation has the least error (RMSE and NRMSE of 1.753 º C and 0.016) but the highest bias (NMBF of 0.013), increasing the resolution reduced the bias greatly (NMBF lowered to 0.007) with no significant increase in error (Table 4).. 230

Relative humidity
The surface relative humidity in Brunei simulated using high-resolution and low-resolution simulations had a very good correlation coefficient (above 0.8). Although WRF_Chem_20km suffered slightly in correlation (by 0.017) compared to WRF-Chem_100km using nearest neighbor estimation technique, the overall performance was better using 20 km horizontal https://doi.org/10.5194/acp-2019-692 Preprint. Discussion started: 23 March 2020 c Author(s) 2020. CC BY 4.0 License. resolution ( lower error and bias) with RMSE, NRMSE and NMBF of 7.544 %, 0.142 and -0.027 respectively. Also, the 235 bilinear interpolation method maintained almost the same correlation for both resolutions and the 20 km resolution simulation had lower error and bias with RMSE, NRMSE and NMBF of 7.293 %, 0.138 and -0.024 respectively (see Table 2). This confirms that higher horizontal grid resolution improves the representation of surface relative humidity in our WRF-Chem simulations.
In the case of air profile relative humidity in Brunei, increased simulation resolution worsened the correlation and gave higher 240 error especially with bilinear interpolation method (RMSE and NRMSE of 16.101 % and 0.186). Although high-resolution simulation with nearest neighbor estimation reduced the bias greatly (NMBF of 0.100), the same simulation with bilinear interpolation method gave the highest bias (NMBF of 0.125) (see Table 3).
Similarly, the air profile relative humidity in Singapore had its performance worsened with higher resolution simulation. The correlation coefficient reduce from 0.844 in the low resolution simulation to 0.823 and 0.829 with nearest neighor and bilinear 245 interpolation method respectively, in high-resolution simulation. Although bias was lowest with bilinear interpolation on highresolution simulation, the errors were higher especially with the nearest neighbor estimation (RMSE and NRMSE of 15.281 and 0.193) ( Table 4).

Wind speed
The surface wind speed in Brunei station underperformed with increased parametrization using high-resolution simulation. 250 Using nearest neighbor calculation, WRF-Chem_20km simulation had lower correlation coefficient relative to WRF-Chem_100km simulation (0.329 compared to 0.469). Also, the error and bias for WRF-Chem_100km simulation were lower than WRF-Chem_20km simulation with RMSE, NRMSE and NMBF of 1.165 m/s, 0.226 and 0.019 compared to 1.221 m/s, 0.237 and 0.101, respectively. Although, bilinear interpolation technique improved high-resolution simulation results but the performance (correlation, errors and bias) were lesser than the calculated results using nearest neighbor method with coarse 255 simulation (see Table 2).
Considering air profile wind speed in Brunei, all the simualtions had a good correlation coefficient (above 0.8). Increased resolution simulation results had worse correlation coefficients but the bilinear interpolation estimation improved the correlation slightly. Overall, the high-resolution simulation gave lower error and bias compared to low-resolution simulation (see Table 3). 260 In Singapore on the other hand, the air profile wind speed simulated using high resolution improved overall results. The correlation coefficient increased to 0.915 from 0.913 and 0.897 for bilinear interpolation and nearest neighbor techniques respectively. Also, the RMSE, NRMSE and NMBF reduced to 3.384 m/s, 0.114 and -0.093 respectively (see Table 4). https://doi.org/10.5194/acp-2019-692 Preprint. Discussion started: 23 March 2020 c Author(s) 2020. CC BY 4.0 License.

Wind direction 265
The transport of aerosol is mainly governed by wind and its direction. The wind direction at 10 meters above ground level using finer resolution had better representation than the low-resolution simulation. Using nearest neighbor method, The correlation coefficient for WRF-Chem_20km was slightly higher than WRF-Chem_100km (0.258 against 0.252). Also, the error and bias observed were lower in high-resolution simulation. WRF-Chem_20km had RMSE, NRMSE and NMBF of 85.082 degree, 0.243 and -0.051 compared to WRF-Chem_100km with 103.342 degree, 0.295 and -0.115, respectively. 270 Overall, bilinear interpolation calculation method yielded better results for high-resolution simulation with RMSE, NRMSE and NMBF of 78.669 degree, 0.225 and -0.039 respectively (see Table 2).
In Brunei, air profile wind direction simulated using high-resolution simulation had worsened correlation but lowered bias compared to low-resolution simulation with either nearest neighbor or bilinear interpolation estimation. The bilinear interpolation calculation greatly improved the results especially for the low-resolution simulation with correlation coefficient 275 of 0.775, and RMSE and NRMSE of 51.222 degree and 0.180 respectively (see Table 3).
Whereas in Singapore, wind direction profile simulated using high-resolution simulation yielded better correlation coefficent (0.681) with either bilinear interpolation or nearest neighbor methods. The bias was also reduced to -0.046 and -0.055 while estimating using nearest neighbor and bilinear interpolation techniques respectively. While the error was lowered with increased resolution simulation, the least error (RMSE and NRMSE of 66.541 degree and 0.208) was obtained when using 280 bilinear interpolation calculation on coarse resolution simulation (see Table 4).
Overall, the model in this study simulated air profile up to 20 km altitude (in the upper troposphere and lower stratosphere region). The air temperature and wind speed profile in Brunei and Singapore had an overall improvement (lower errors and bias/higher correlation) with finer resolution simulation. Wind direction profile using high resolution simulation suffered a 285 slight increase in errors, but with higher correlation in Singapore and lower bias in both Brunei and Singapore.
The high-resolution WRF-Chem simulations performed better in meteorology representation, though the low-resolution simulations results were also very good. Thus, the model developed and presented in this work proved to be an effective tool in simulating atmospheric conditions and deep convection in the South-East Asia region.

Air Quality Evaluation 290
The simulated aerosols were evaluated by comparing the results to the observed mean daily particulate matter concentrations in Singapore and Brunei Darussalam. in Singapore and PM10 in Brunei, respectively (Table 5) (Figures 6 and 7). The parameterization involved in this coarse simulation is relatively minimal and over-approximation of topography could have induced higher inconsistencies in the spatial distribution of particulate matter across the region. 300 In the case of the WRF-Chem_20km simulation, the simulated PM2.5 in Singapore had a correlation coefficient of 0.8864 and the simulated PM10 levels in Brunei had 0.3469 correlation coefficient (Table 5) (Figures 6 and 7). Despite that particulate matter levels were generally underestimated across the region in 20 km grid scale than 100 km grid scale simulation (larger bias, RMSE and NRMSE), the improvement in correlation showed that the higher horizontal resolution with 20 km grid scale had a better representation of PM. 305 Based on the current knowledge of processes (physical and chemical) and observations in estimating biomass burning emissions, many uncertainties still existed, for which enhancements were required and adopted in several models to scale the bottom-up inventories to the top-down constraints (Kaiser et al., 2012). Wu et al (2011) (Kaiser et al., 2012;Wu et al., 2011). In order to study this effect, the fire emissions were elevated 310 by changing the default biomass burning emission enhancement factor from 1.3 to 6 in simulatng the June 2013 haze episode in the region, and the new results were compared with default coarse and finer resolution simulations (see Figures 4 and 5).
After enhancement of biomass burning emissions, the simulation (WRF-Chem_20kmX) gave a very good representation of particulate matter distribution across the South East Asia region. The PM10 levels in Brunei were adequately represented in this simulation with a correlation coefficient of 0.5974. Although PM2.5 levels in Singapore were slightly overestimated, the 315 correlation coefficient was maintained at 0.8847 (Table 5) (Figures 6 and 7).
Moreover, it is evident from the comparison of the time-averaged total column aerosol optical depth at 550 nm of the simulations against reanalysis from MERRA-2 (M2T1NXAER v5.12.4) that the default factor of 1.3 for biomass burning emission in the model underestimated the aerosols and could not repeat the AOD reanalysis in both lower and higher horizontal grid scale simulations. On the other hand, the enhacement factor of 6 used in high resolution simulation fairly repeated AOD 320 reanalysis with positive bias and overestimation of the aerosols especially across the highly polluted zones in the South East Asia region (see Figure 8). From this analysis,it could be suggested that an optimum enhancement factor for fire emissions to simulating aerosols in this region should be above 1.3 and below 6. https://doi.org/10.5194/acp-2019-692 Preprint. Discussion started: 23 March 2020 c Author(s) 2020. CC BY 4.0 License.
The biomass burning emissions from Sumatra was very intense in June 2013 and it was mainly responsible for the haze and increase in PM concentrations in Singapore. It was observed that the fire emissions elevated around the eastern part of Sumatra 325 on the 19th of June 2013 which is in correlation with the increased frequency of hotspots captured by MODIS. The winds blowing in the northeast and easterly direction over Sumatra and Peninsular Malaysia through the South China Sea is mainly responsible for the long-range transport of particulate matter to Singapore and Borneo region from the biomass burning sources in Sumatra. Some sparsely spread fires were also detected in the Sarawak region in Borneo and they had more influence on air pollution levels in Brunei during this period. This spatial distribution was clearly seen and evident from the simulations with 330 20 km horizontal resolution than 100 km (Figures 4 and 5).The elevated fire emissions and air pollution level in this region persisted until the 23rd of June when it peaked and subsequently started its decline. The wind direction gradually shifted around 23rd June from northeasterly direction to the north over the South China Sea, while westerly winds started blowing from eastern Borneo and weak winds prevailed over the Borneo region. The decline in fire emissions and change in wind direction could be responsible for the decline in particulate matter levels observed in Brunei Darussalam after this period. 335

Summary and Conclusion
We have simulated the particulate matter distribution from biomass burning emissions across South East Asia during the June 2013 haze event. We synergistically leveraged on hot-spot datasets across the region from MODIS and VIIRS and assessed the effect of horizontal grid scales and enhancement factor for biomass burning emissions for particulate matter simulations. 340 The result of the statistical metrics for the model simulations against the observation in Brunei and Singapore revealed that the model simulated the meteorology (surface and vertical air profile for variables such as temperature, relative humidity, wind direction and wind speed) adequately across the region.
It was also confirmed from the analysis that higher resolution simulations gave a better representation of some prognostic meteorological variable while some underperformed. Interpolation method for obtaining simulation results was also varied 345 between nearest neighbor and bilinear interpolation. Considering surface meterology evaluations using nearest neighbor method, higher resolution yielded lower bias and errors in all meteorological variables assessed except wind speed. On the other hand, increased resolution simulation resulted in lesser correlation in all variables except wind direction.
On the other hand, bilinear interpolation performed better than nearest neighbor technique for surface meteorology at Brunei in higher and lower model resolutions simulations except wind speed. It proved to be an efficient method for estimating station 350 variables such as surface temperature on low resolution simulation and it gave overall best performance in obtaining surface relative humidity and wind direction in Brunei station. However, nearest neighbor extraction seemed to yield better result for wind speed in Brunei. Therefore, we can conclude that the model succeeded in capturing the particulate matter emissions across the region during the haze episode but had a large underrepresentation in some parts such as the Borneo region. The discrepancies observed may 365 be as a result of errors present in the biomass burning emission inventories. Despite leveraging on VIIRS to augment MODIS hot-spot data, it was very likely that some burning hot-spots were not captured due to cloud cover and vegetation canopies around the Borneo region.
In future research, improvement in the simulation could be made by employing finer emission inventory and ground measurement of biomass burning emissions for WRF-Chem model emissions input to overcome the undetected emissions 370 from satellites. Also, to improve the particulate matter representation, the emissions in the model require enhancement. A biomass burning enhancement factor above 1.3 and below 6 could be included to bring the bottom up inventories to the observed values. Furthermore, discrepancy caused by the underrepresentation of the model resolutions should be avoided by adopting a finer horizontal resolution in resolving the small-scale circulation and the cloud convection while considering prognostic variables staggering on the grid. These could improve significantly the understanding of particulate matter 375 distribution from biomass burning emissions across the region.

Author Contributions
ARA, with support from LD, conducted the simulations, processed the data, and prepared the data visualization with contributions from all authors. LD, MIP, and LCD sourced for all observation datasets required. ARA prepared the paper and