Improving prediction of trans-boundary biomass burning plume dispersion: from northern peninsular Southeast Asia to downwind western North Pacific Ocean

Plumes from the boreal spring biomass burning (BB) in northern peninsular Southeast Asia (nPSEA) are lifted into the subtropical jet stream and transported and deposited across nPSEA, South China, Taiwan and even the western North Pacific Ocean. This paper as part of the Seven SouthEast Asian Studies (7-SEAS) project effort attempts to improve the chemical weather prediction capability of the Weather Research and Forecasting coupled with the Community Multiscale for Air Quality (WRF–CMAQ) model over a vast region, from the mountainous near-source burning sites at nPSEA to its downwind region. Several sensitivity analyses of plume rise are compared in the paper, and it is discovered that the initial vertical allocation profile of BB plumes and the plume rise module (PLMRIM) are the main reasons causing the inaccuracies of the WRF–CMAQ simulations. The smoldering emission from the Western Regional Air Partnership (WRAP) empirical algorithm included has improved the accuracies of PM₁₀, O₃ and CO at the source. The best performance at the downwind sites is achieved with the inline PLMRIM, which accounts for the atmospheric stratification at the mountainous source region with the FINN burning emission dataset. Such a setup greatly improves not only the BB aerosol concentration prediction over near-source and receptor ground-based measurement sites but also the aerosol vertical distribution and column aerosol optical depth of the BB aerosol along the transport route. The BB aerosols from nPSEA are carried by the subtropical westerlies in the free troposphere to the western North Pacific, while BB aerosol has been found to interact with the local pollutants in the Taiwan region through three conditions: (a) overpassing western Taiwan and entering the central mountain area, (b) mixing down to western Taiwan, (c) transport of local pollutants upwards and mixing with a BB plume on higher ground. The second condition, which involves the prevailing high-pressure system from Asian cold surge, is able to impact most of the population in Taiwan.

1 Introduction

Large quantities of gaseous and aerosol pollutants released from biomass burning affect regional air quality, radiative forcing, public health and economic burden, especially in Southeast Asia (Chen et al., 2017; Lee et al., 2017; Pani et al., 2018, 2020). The prolonged heat during the dry season (December to May) in peninsular Southeast Asia (PSEA) has led to the deterioration of biomass burning (BB) in northern PSEA (nPSEA) (Kim Oanh and Leelasakultum, 2011). The outflow of the BB smoke plumes from nPSEA usually occurs during the spring season (late February until mid-April), when the high-pressure system has retreated northwards back into the Asian continent. The mountainous structure over the north-central PSEA has lifted the BB plume into the subtropical Pacific High (700 to 800 hPa, ∼1–3 km) under prevailing southern wind (Dong and Fu, 2015b; Huang et al., 2020). The plume is then transported eastward to the West Pacific and frequently detected at the Lulin Atmospheric Background Station (LABS) in central Taiwan (Fu et al., 2012; Lee et al., 2011; Lin et al., 2017, 2014, 2013; Ou-Yang et al., 2014; S.-H. Wang et al., 2013). Moreover, there were several instances when the high-pressure system entered Taiwan and brought the upper-layer BB plumes down to populous southwestern Taiwan and altered the atmospheric chemistry and composition (Dong et al., 2018; Huang et al., 2016; Yen et al., 2013).

Spaceborne remote-sensing data from satellites and the high-resolution spatiotemporal data generated from the chemical weather prediction (CWP) model are often used for studying long-range transport of BB smoke across the region (e.g., Huang et al., 2020; Tsay et al., 2013). Previous studies have found that the numerical model is prone to overestimating the BB emissions including CO, PM_2.5 and PM₁₀ up to 3 times the measured amount at the major burning source in northern Thailand (Huang et al., 2013; Pimonsree et al., 2018). The exceedance of estimated emission at the near-source burning leads to the incorrect modeled signal at the downwind site (Fu et al., 2012). The modeled columnar aerosol optical depth (AOD) is found to be comparable with aerosol products of the Aerosol Robotic Network (AERONET) and the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor as well as columnar CO and NO₂ at the burning source over the nPSEA region, but great discrepancies are found for the spatial distribution of downwind plumes (Dong and Fu, 2015b; Fu et al., 2012). In those models, the vertical distribution percentage of BB emission was set to be constant throughout the case. However, there are many possible factors that govern the actual plume rise condition, including the fire size, vegetation cover, buoyancy heat flux, wind drag, boundary layer condition, etc. (Freitas et al., 2010; Kukkonen et al., 2014; Paugam et al., 2016; Val Martin et al., 2012). Furthermore, the accuracy of the model depends greatly on the plume rise condition.

As part of the local effort of the interdisciplinary Seven SouthEast Asian Studies (7-SEAS) project (Lin et al., 2013; Reid et al., 2013), this paper attempts to improve the modeling performance of the long-range transport of BB from the nPSEA region to the downwind region using the Weather Research and Forecasting coupled with the Community Multiscale Air Quality (WRF–CMAQ) model. The paper attempts to improve the ability of the CMAQ model and its plume rise module (PLMRIM) to predict the complexity of BB amount from its burning source in nPSEA to its downwind receptor LABS. With the availability of on-site and satellite lidar (light detection and ranging) measurement, the vertical plume rise profile can be better understood to ensure that BB plumes are distributed according to the actual conditions (Walter et al., 2016; S.-H. Wang et al., 2013). In this work, several factors, including the injection height, initial vertical distribution and smoldering fraction, are considered in the model. Knowing that the atmospheric circulation over nPSEA is also affected by terrain, the work now intends to incorporate the interaction of the atmospheric stratification and BB plumes into the PLMRIM. This research takes an approach mainly from the perspective of the vertical distribution profile of modeled BB aerosol concentration with the assistance of top-down and bottom-up vertical lidar profilers. The better-performing setting will be applied to test its applicability and to dissect the sources of high pollution at LABS and in western Taiwan.

The model experimental design (Sect. 2.1), model emission input (Sect. 2.2) and case study setup (Sect. 2.3) are explained in detail. The performance of the PLMRIM is then verified with ground-based measurement stations in Sect. 3.1 and vertical aerosol products from lidar sensors (MPLNET, CALIPSO) and MODIS columnar AOD (Sect. 3.2), where the reliability and accuracy of inline PLMRIM are discussed (Sect. 3.3). The resulting output is subsequently studied in Sect. 4 to answer the transport mechanism to the ground-based observation sites in western Taiwan. Conclusions to these findings are made in Sect. 5.

2 Methodology

The study focuses on the spring BB events in March 2013. With moderate burning occurring in nPSEA, this El Niño–Southern Oscillation (ENSO)-neutral year is chosen because the LABS mainly received the BB plumes with minimal influence from the Asian dust storm to Taiwan (NOAA-ESRL, 2020; TAQM, updated daily; Kong et al. 2021). The 7-SEAS spring campaigns carried out during the BB season supply an abundance of air quality and aerosol measurement data (http://rsm2.atm.ncu.edu.tw/, last access: 17 August 2021).

Figure 1 (a) Domain setup of model (domain 1–4) with terrain height information; (b) third domain covering Taiwan (d03) with information of terrain height (contour fill), AA' cross-section (dotted red line), and locations of Taiwan EPA air quality and CWB weather stations (black dots) and LABS receptor site (big red dot); (c) fourth domain covering part of nPSEA (d04) with terrain height (contour fill), BB' cross-section (dotted red line), and location of Thailand PCD ground air quality stations (black dots) and the DAK source site (big red dot). The MPLNET data collected are located at the DAK station (big red dot).

2.1 Model physics and experimental design

This work employs the Weather Research and Forecasting with the Advanced Research core (WRF–ARW v3.9.1) (Wang et al., 2017) model to hindcast the weather field and predict the corresponding air chemistry field with the chemical transport model CMAQ v5.2.1 (Byun and Schere, 2006). The model domain is dynamically downscaled through nesting from the majority of Asia (d01 resolution: 45 km) to cover the transport route from nPSEA to Taiwan (d02: 15 km), Taiwan only (d03: 5 km) and nPSEA only (d04: 5 km), as shown in Fig. 1. The weather input for the initial and lateral boundary condition is the 6-hourly $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ National Centers for Environmental Prediction (NCEP) Final Analyses (FNL) dataset (NCEP-ds083.2, 2000, updated daily). As an extension of the latter, data assimilation is applied for both grid and observation nudging. The weather data for observation nudging are obtained from the National Centers for Environmental Prediction (NCEP) Automated Data Processing (ADP) Global Surface (NCEP-ds461.0, 2004, updated daily) and Upper Air Observational Weather Data (NCEP-ds351.0, 2000, updated daily), with additional local sites operated by the Taiwan Central Weather Bureau (CWB) and the Thailand Pollution Control Department (PCD). The radii of influence (RIN) for both d03 and d04 are updated to 100 km based on the average distance between the observation stations (d03: 125 km; d04: 153 km) and the minimum distance between two stations (d03: 64 km; d04: 36 km). Wind speed and wind direction are substantially improved by observation nudging. A detailed discussion about meteorology performance is given in Sect. A in the Supplement. Other WRF–CMAQ settings and configurations are listed in Table 1.

On top of the ground-based measurement weather and air quality data, the lidar systems are also used to evaluate the performance of the model ability to estimate the vertical profile of BB aerosols. They are the bottom-up Micro-Pulse Lidar Network (MPLNET) and top-down Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) lidar sensors. The MPLNET is a federated network managed by NASA to measure the aerosol vertical structure (Welton et al., 2000). In line with the 2013 7-SEAS spring campaign conducted in nPSEA, the MPLNET device is located at the Doi Ang Khang (DAK) meteorology station to collect the near-source aerosol vertical distribution profile (L1.5a) data. The gridded extinction, diagnosed from the planetary boundary layer height and vertical aerosol extinction coefficient data collected, is used to verify the performance of the model output (S.-H. Wang et al., 2015). The CALIOP sensor mounted on the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) satellite is used to study the transport pattern with larger spatial coverage to complement the single-point cross-extinction profile provided by the MPLNET system. The diagnosed vertical feature mask (VFM) product is used to distinguish the aerosol types with consideration of observed backscatter strength and depolarization (Winker et al., 2010).

Table 1 WRF and CMAQ model settings.

2.2 Emission data

2.2.1 Anthropogenic and biogenic emission inventories

The anthropogenic emissions are re-gridded for the first, second and fourth domain (d01, d02 and d04 in Fig. 1) from the MIX dataset available at $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ for the year 2010 (Li et al., 2017; Zheng et al., 2018). The Model of Emissions of Gases and Aerosols from Nature (MEGAN v2.10) produces the biogenic emission input (Guenther et al., 2012) using the updated 8 d averaged leaf area index (LAI) (Yuan et al., 2011) and present-day plant functional types (PFTs) from the Community Land Model version 4.0 (CLM4.0) (Oleson et al., 2010). The third domain (d03) covering Taiwan uses the 2010 anthropogenic and biogenic emissions from the locally developed Taiwan national emission database (TEDSv8.1) (TEPA, 2017). Except the high quality of the East Asia national emission inventories (China, Taiwan, Japan and Korea), large uncertainties in Southeast Asia emission due to the scarce availability of region-specific emission factor are pointed out by the inventory developers (Kurokawa et al., 2013; Li et al., 2018; Ohara et al., 2007) and local modeling efforts (Dong and Fu, 2015a; Ooi et al., 2019). Such inaccuracies are likely to affect the performance of further modeling work in the area. Therefore, energy statistics based on the global anthropogenic emissions dataset Evaluating the Climate and Air Quality Impacts of Short-Lived Pollutants (ECLIPSE), developed by the International Energy Agency (IEA) (Klimont et al., 2017), is used in place of the MIX dataset for PSEA. The accuracy deviation between these two datasets in nPSEA is determined through the WRF–CMAQ model performance in Sect. 4. The detailed comparison of ECLIPSE and the MIX dataset in 2010 is shown in Sect. B in the Supplement.

Table 2 Case setup to evaluate PLMRIM performance.

2.2.2 Biomass burning emission inventory

The study region is composed of small fires with a small area burnt but has a rather substantial amount of fuel load and BB emissions due to the high woody compositions of the tropical and temperate forest covers. The global dataset Fire INventory from the NCAR (FINN v1.5; referred as “FINN” here onwards) has been applied in several previous works focusing on the region (Lin et al., 2014; Pimonsree and Vongruang, 2018) and is used as the input to the BB emission inventory into the model. A particular comparison work done for 2014 biomass burning episodes has shown that FINN, when used with NCEP FNL boundary condition, gives the greatest accuracy for PM₁₀ at the source region compared to the GFEDv4.1 fire emission dataset (Takami et al., 2020). Seeing that the temporal speciation is handled in this research work, the main difference between fire emission inventories is the total amount of emission produced (Liu et al., 2020); hence this paper will settle with the regionally more robustly tested FINN dataset for the subsequent studies. FINN is a 1 km×1 km resolution bottom-up daily emission dataset produced from the MODIS product of active fire, land-cover type and continuous vegetation field (Wiedinmyer et al., 2011). Each active fire is assumed for a 1 km² burnt area, and the emission factor is geographically and land-cover-dependent. The BB emission is processed with the fire_emis preprocessor to allocate to each grid and specify input into the WRF–CMAQ model at the hourly scale.

2.3 Case study setup

The plume rise module (PLMRIM) derives the initial plume top and bottom, plume rise, and its dispersion according to the atmospheric stability and its residual buoyancy flux (Kukkonen et al., 2014). Among a wide range of PLMRIM approaches, the simplest plume rise allocation method is the direct allocation of the initial plume top and bottom through prescribed height for all fires. This is the conventional method adopted in the case study region (Chuang et al., 2016b; Pimonsree et al., 2018). They can be determined on the basis of fixed height (J. Wang et al., 2013) and an empirical ratio of the plume height allocation (WRAP, 2004), and they can be adjusted with, for example, the stereo-height data from the space-based multi-angle imaging spectroradiometer (MISR) (Jian and Fu, 2014; Val Martin et al., 2012). The inline plume rise algorithm couples the interaction of BB plume dispersion with the basic weather dynamics to determine the effective plume rise height and subsequently the plume top and bottom. This inline PLMRIM is also able to resolve the fire on the sub-grid scale and feed the plume dynamics information back into the atmospheric dynamics (Gillani and Godowitch, 1999). However, the more complex the PLMRIM gets, the higher the quality and quantity of input data are required to be to ensure its reliability.

Figure 2 Initial CO emission rate (mol s⁻¹) profile at Mae Hong Son, Thailand, on 13 March 2013 (UTC) for each case setup in Table 2 with (a) F0, (b) F800, (c) F2000, (d) FWrp, (e) IDef and (f) IWrp/IWrp + EC.

In this work, combinations of injection height, initial vertical distribution, smoldering fraction, and offline and inline PLMRIM are tested to determine the more suitable settings for prediction of plume rise. Five case studies are set up for the evaluation of plume rise performance, and their respective initial plume rise profiles are shown in Table 2. The Nofire case represents the pollution condition when no BB emission is included, while the others allocate the BB emission from the FINN dataset. F0, F800 and F2000 represent the offline PLMRIM where the injection height is fixed at the generally accepted near-surface layer, 800 and 2000 m (J. Wang et al., 2013). This fixed-height method controls the plume top to be consistent; hence there is no hourly and daily variation in the plume top throughout the simulation period. FWrp uses the WRAP empirical equation to allocate the initial plume rise (WRAP, 2004). The plume top and bottom vary hourly with the buoyancy efficiency, with higher plume height during the hotter noontime, as illustrated in the initial plume profile in Fig. 2 (FWrp). However, the empirical ratio adopted for each burning grid is the same every day. IDef is the inline plume-in-grid system that comes with the CMAQ model (Gillani and Godowitch, 1999). Fire emission is fed into the model at each grid point, with plume top and bottom calculated through interaction of plume buoyancy efficiency and atmospheric stratification. The vertical distribution of CO plume on 12 March 2013 is shown in Fig. 2 (IDef), but the daily weather condition is expected to vary the vertical distribution. IWrp has updated IDef with the WRAP empirical specification on fire size. In this case, the plume can be distributed according to the diurnal buoyancy efficiency and near-surface smoldering fraction as specified by WRAP. With a more reasonable BB plume peak at noontime in Fig. 2 (IWrp), it is expected to improve the near-source concentration prediction of the model as seen from the initial plume profile. IWrp + EC is the same as IWrp but with the anthropogenic emission in PSEA replaced by the ECLIPSE dataset, as specified in Sect. 3.2.1. The initial emission profiles (within plume top and bottom) of all cases are distributed evenly according to the height of each vertical layer.

Table 3 Performance of modeled chemistry field with different settings of PLMRIM at a mountain site in the western North Pacific (LABS) and nPSEA (DAK). R: correlation coefficient; MFB: mean fractional bias; MFE: mean fractional error; MNB: mean normalized bias; MNE: mean normalized error. Bold values are model output that satisfied the standard of each index.

Figure 3 Comparison of PLMRIM (observation (black), nofire (blue), FWrp (green), IDef (orange), IWrp + EC (red)) of (a) hourly wind field and PM_2.5 at DAK and (b–d) hourly wind field and (b) PM₁₀, (c) CO and (d) O₃ at LABS in March 2013; gray shade highlights the high-pollution hour at LABS (CO > 300 ppb, PM₁₀>35 µg m⁻³). Wind field for observation (black) and simulation (red) is shown in vector form.

3 PLMRIM performance

3.1 Ground-based measurement stations

The model output is compared with the measurement data at a high-altitude background mountain station in the western North Pacific, the Lulin Atmospheric Background Station (LABS) (receptor; 2862 m a.m.s.l.; 23.47^∘ N, 120.87^∘ E), and Doi Ang Khang Meteorology Station (DAK) (source; 1536 m a.m.s.l.; 19.93^∘ N, 99.05^∘ E), marked in Fig. 1b and c. The DAK station is an upwind near-source BB location in nPSEA, located in Chiang Mai Province, Thailand, close to the border of Myanmar and Thailand. It is located away from the cities and mainly received air masses from burning regions in the upwind area (Hsiao et al., 2016; Pani et al., 2016), which made this site representative of the BB emissions from Myanmar, on the western side of Thailand (Khamkaew et al., 2016; S.-H. Wang et al., 2015). The hourly PM_2.5 data from the DAK station were collected during the 2013 7-SEAS spring campaign. Table 3 shows the performance of PLMRIM for daily PM₁₀, daily PM_2.5, hourly O₃ and hourly CO at LABS and DAK according to the model benchmark (correlation coefficient, R; mean fractional bias, MFB; mean fractional error, MFE; mean normalized bias, MNB; mean normalized error, MNE) suggested by the Taiwan EPA (see Sect. C in the Supplement). MFB results show that the pollutants are generally overestimated at these mountain stations. Unlike the case in the maritime continent that worked best with the F800 method (J. Wang et al., 2013), neither of the fixed-height methods (F0, F800, F2000) apply well for the nPSEA region. Only slight improvement is observed for the offline module (FWrp), with injection height varying according to the fire size. The inline modules (IDef, IWrp) have obvious improvement at both LABS and DAK. For the ground stations in Taiwan and Thailand (black markers in Fig. 1b and c), all models have underestimated the pollutant concentrations, while the IWrp has performed better than the default inline mechanism with higher correlation attained. The daily PM₁₀ at the northern Thailand PCD source stations for IWrp achieved R=0.84, an improvement from R=0.77 of FWrp, while daily PM_2.5 at the Taiwan EPA ground stations for IWrp achieved R=0.46, an improvement from R=0.26 of FWrp (see Table C1 for a detailed comparison). Adjustment of anthropogenic emission with ECLIPSE data (IWrp + EC) shows clear improvement of CO, especially at the stations in Taiwan but not in Thailand. The comparably insignificant emission amount of anthropogenic emission compared to the BB emission at the near-source BB sites in Thailand is attributed to the minor pollutant changes during the BB period.

Among all, the inline modules (IDef, IWrp, IWrp + EC) give the lowest bias and closest correlation with the measured ground station. This highlights the importance of atmospheric-stability-based PLMRIM to capture the plume rise variation at the source site. The boundary layer evolution throughout the day is very much distinctive for mountain valleys compared to the flat surface where burning usually happens. As highlighted previously (Chuang et al., 2016a; Dong and Fu, 2015b), the geographical lifting mechanism in nPSEA is the main factor enabling the BB emission to be carried into the subtropical westerlies and hence captured by LABS. Due to the similar performance among the offline and inline settings, the best-performing setup of the offline module (FWrp) and inline module (IWrp + EC) is selected to simplify the subsequent discussion.

Figure 3 shows the time series plots for the hourly wind field and PM_2.5 at the DAK source site and hourly wind field, PM₁₀, CO and O₃ at LABS. The high-pollution episode (marked by gray shade) fits well with the great contrast between the model fire and nofire scenarios, thus confirming that BB plumes are the main pollution source in the high-pollution episodes. From the time series plot, the hourly PM_2.5 at DAK (Fig. 3a) and hourly PM₁₀ (Fig. 3b) at LABS are well captured by the inline module compared to the offline counterparts. In Fig. 3b, the wind direction shifted to strong southwesterlies in the second half of March. It is followed by a rise in pollution level at LABS. The offline module (FWrp) has significantly overpredicted PM₁₀ at some peaks, even up 200 µg m⁻³. Fair agreement is obtained for CO (Fig. 3c) and O₃ (Fig. 3d), with slight overestimation when concurrent high PM₁₀ is modeled. Short-term peak values of 4–5 h are observed in all models for PM₁₀, CO and O₃. The systematic errors for these pollutants at the peak points are believed to be the uncertainties involving the FINN BB emission (Pimonsree et al., 2018). It is found that the performance of O₃ is relatively unaffected by the PLMRIM choice.

3.2 Aerosol vertical distribution

As illustrated in the shaded region in Fig. 3, the major period that affected LABS is during 18–28 March 2013. The transport time is known to be around 2–3 d (Chuang et al., 2015), but a longer time of 4–5 d is taken to account for the BB emission generation, lifting and dispersion at the source site. Hence, the vertical profile of the extinction coefficient of the model output at the DAK station is compared with the ground MPLNET L1.5a lidar data during 13–28 March in Fig. 4. In Fig. 4a, the MPLNET extinction coefficient has shown that the aerosol layer presents around 3–4 km height and is mainly confined below the boundary layer. Despite the lower concentration on the surface than the boundary layer, the concentrated amount of aerosol detected by MPLNET (14–15, 19–22 and 23 March) still agrees well with the DAK time series data in Fig. 3a.

Figure 4 (a) Vertical extinction coefficient profiles between 13 and 28 March 2013 at the DAK station from MPLNET with boundary layer height (white), (b) MPLNET 3-hourly average extinction coefficient, (c) IWrp + EC 3-hourly averaged model output, (d) FWrp 3-hourly averaged model output.

The 3-hourly average profile of the extinction coefficient from MPLNET data as well as IWrp + EC and FWrp model output during 13–28 March 2013 at the DAK station is illustrated in Fig. 4b–d. In Fig. 4b, the MPLNET extinction coefficient is low at the surface and peaks between 2.5–3.2 km. The model output has a lower elevation over the DAK station and has modeled a higher extinction coefficient, which is likely to be the accumulation effect due to lower-wind conditions. The model output in Fig. 4c and d shows that the maximum layers above the presumed cap (3.2 km) occurred most prominently during the evening to midnight and more often in offline than inline modules. The model shows that the offline module gives a time-invariant large value over the entire layers, while the inline module gives a greater approximation of the diurnal variation with the MPLNET result throughout the day. Therefore, during the daytime, the offline module has produced a higher plume height than the CMAQ inline module (Guevara et al., 2014). However, the 3-hourly averaged temporal variation in extinction coefficient of inline output fits well with that observed by the MPLNET system.

The extinction coefficient from MPLNET and model output data are only available for qualitative comparison due to their generically different derivations. The lidar system determines the extinction coefficient through the backscatter feedback from the release of the laser beam at 527 nm at every minute, while the CMAQ model used the mass reconstruction method to sum up the extinction coefficient of each model aerosol species in each layer (Mebust et al., 2003). The empirical assumption for each species and the lower vertical model resolution is attributed to the uncertainties in the modeled extinction coefficient that is typically higher than the value retrieved by MPLNET.

Figure 5 CALIOP vertical feature type and aerosol subtype during a continuous episode starting from (a–c) 19 March (06:02 LST), (d–f) 19 March (13:42 LST) and (g–i) 20 March (02:07 LST). The corresponding position of the satellite swath is marked in points of red and gray marked in (a, d, f) and altitude below 0 km in (b, c, e, f, h, i). Feature type: 0 – invalid, 1 – clear air, 2 – cloud, 3 – aerosol, 4 – strato, 5 – surface, 6 – subsurface, 7 – no signal; subtype of feature: ND – no data, 1 – marine, 2 – dust, 3 – polluted continental, 4 – clean continental, 5 – polluted dust, 6 – smoke.

Figure 6 Comparison of model PM₁₀ (µg m⁻³) cross-sectional profile corresponding to CALIPSO period and swath in Fig. 4. The range of the left panel is 0–120 µg m⁻³, and the range of the right panel is 0–300 µg m⁻³.

Figure 5 shows the CALIOP VFM at the midpoint of the BB pollutant transport route to the receptor during one of the episodes on 19–20 March 2013. On the morning of 19 March, the sensor (swath: Fig. 5a) captured the smog layer at a height of 4 km a.m.s.l. (above mean sea level) over the mountainous region (Fig. 5b and c). The aerosols detected are mainly made up of smoke and mixed polluted continental aerosols, which is the main burning emission source. It is known that the burning aerosols from the western part of nPSEA are orographically lifted by west-to-southwesterlies to a higher altitude depending on the terrain height (Cheng et al., 2013; Z. Wang et al., 2015). For the swath in Fig. 5d–f, the aerosol layers are detected at high levels up to 4 km during the midday. It is most certain to be transported over from nPSEA since the aerosol layer is detected over the sea, where burning does not occur. Secondly, the plume thickness is around 4 km despite the flat land surface, which is much higher than the source site, which usually ranges between 0–3 km. The aerosol layers are believed to be lifted to a higher level and also mixed to the surface over the land mask in southeastern China, which is later confirmed in the model result in Sect. 4. This region locates one of the largest cities and main industrial bases in Asia, the Pearl River Delta (PRD), which produces a large amount of anthropogenic emission. The potential vertical mixing is very likely to pick up the pollutants from the industrial base into the aerosol plume. Recently, it is proven through brute-force methods that the pollution from the PRD cluster arrived at the higher altitude in Taiwan during the winter season (Chuang et al., 2020). About 12 h later, when the swath (Fig. 5g–i) moved closer to Taiwan, the plumes moved towards north of 16^∘ N but still maintained a similar altitude that could be detected by the LABS station at ∼2.4 km a.m.s.l. in the model (Fig. 1). The plume is also found to continue gaining moisture content along the path.

A detailed comparison of vertical distribution for all sensitivity tests is given in Sect. D in the Supplement, but here we continue to discuss the FWrp and IWrp + EC cases. In general, the offline FWrp produces a much higher concentration of high-PM₁₀ aerosol layers compared to the inline IWrp + EC. Figure 6 shows the model PM₁₀ result for FWrp (range: 0–300 µg m⁻³) and IWrp + EC (range: 0–120 µg m⁻³) for the corresponding period of CALIPSO swath in Fig. 5. Comparison of Fig. 6a–d shows that the FWrp produces higher plumes, and IWrp + EC produces lower plumes since the former produces the initial plume profile on 19 March that is consistently high and less dependent on the atmospheric stability induced by mountain flow (Fig. S3). Farther from the source site (Fig. 6e and f), both runs predict a much lower aerosol layer around 2 km compared to the 4 km height captured by the CALIOP sensor. The underrepresentation of both systems along the transport path above the sea might be due to the moisture detrainment and entrainment process that is not accounted for in the current model (Paugam et al., 2016; Sofiev et al., 2012).

With a concentration difference of more than a factor of 2 between FWrp (up to 300 µg m⁻³) and IWrp + EC (up to 120 µg m⁻³), a more accurate value is captured at LABS by the IWrp + EC, as shown in Table 3. Regardless of the PLMRIM used, the top height of the plume is confined by an overhead upper-layer wind system. The system has created a strong shear and suppressed the lifting pertaining to the burning convective heat. This explains the invariant of plume height when different settings are used.

Figure 7 Comparison of daily total column AOD on 20 March (10:30 LST) of model output (a) IWrp + EC, (b) FWrp and (c) nofire with (d) MODIS data from Fig. 5. Vector profiles given in (a–c) are the surface wind profile.

The cross-sectional profile of PM₁₀ in Fig. 6 shows that the amount of emission produced by the offline method is substantially larger than the amount produced by the inline method. However, the vertical PM₁₀ value could not be verified due to the lack of measurement of vertical distribution of PM₁₀. The amount of PM₁₀ has directly contributed to the columnar AOD value, and the latter could serve as a good benchmark for the accuracy of model aerosol concentration. Hence, the total columnar AOD data provided by the $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ MODIS Terra Level 3 AOD product (MOD08_D3; Platnick et al., 2015) during the same period (20 March, 10:30 LST) are used for the verification of the aerosol concentration through the columnar AOD value. Figure 7 shows the second model domain (d02), which covers the transport route between the source (d04) and the receptor (d03) domains. The comparison between Fig. 7a and c is able to show the difference between fire and nofire cases, which is solely contributed by the biomass burning plumes from nPSEA. The figure also shows that the total column AOD produced by the inline module gives an approximation closer to that of MODIS. FWrp greatly overestimates the aerosol produced by the BB emissions, while the inline module gives a closer agreement in northern Thailand and southern Vietnam.

3.3 Reliability of inline PLMRIM

The variation in model performance might be caused by the compatibility of the emission inventory with the PLMRIM performance. The FINN dataset provides higher-resolution data for each fire (1 km²) compared to the other emission dataset (GFEDv4s: 0.25^∘; GFASv1.2: 0.1^∘). As the finest study domain at the burning source is downscaled to 5 km, the FINN dataset would have the nearest representation of the emission grid distribution. BB emission in nPSEA is mainly caused by small fires and prevailing dry conditions over the period (Giglio et al., 2013; Reid et al., 2013); hence the representation of the small fires (usually accounted from 500 m burnt area) in the emission inventory is relatively crucial. This might have been one of the reasons that it fits better in the inline calculation with the plume-in-grid concept. When the offline method is adopted (FWrp), the FINN emission dataset in the nPSEA region tends to overpredict by a factor of 4 (Fig. 3a). Previous literature has had to make an adjustment to the fire inventory to bring down the FINN emission amount of PM_2.5 and PM₁₀ that was overestimated by a factor of up to 2 to 3 at the source region (Pimonsree et al., 2018), and FLAMBE overestimates CO and PM₁₀ up to a factor of 3 at the LABS site (Chuang et al., 2015; Fu et al., 2012). From this study, it is seen that the prescribed heights in the offline method have overestimated the plume rise height under the dry weather condition, where the atmospheric stratification damps the pyro-convection through entrainment, while the inline module (IWrp + EC) considers the variability in atmospheric condition over the mountain region better.

Figure 8 Spatial distribution of near-surface PM₁₀ concentration on 19 March 17:00 LST over burning regions of nPSEA for the fourth domain (d04).

The inaccuracy of the offline module is likely to be caused by the role of the complex terrain in uplifting the smoke plume and the nature of the fuel loadings. The connecting slopes (0.2–1.8 km, as seen in Fig. 1c) cause the complication to boundary layer physics that governs the dynamics to transport the plumes formed in the valley pockets. Due to the unique topographic structure in nPSEA, the lifting and breaking away of burning emission plumes from burning area occur during the evening to night period. Therefore, mountain meteorology played an important role in the distribution of higher-level plumes. Moreover, the ability of PLMRIM to capture the boundary layer physics becomes essential in the mountainous region. Through the inline module with the WRAP initial plume profile (IWrp + EC), the natural buoyancy of fire together with the convective interaction of the atmosphere can correctly distribute the BB emission. The spatial distribution of PM₁₀ over burning regions in nPSEA is shown, with comparison made for the scenarios nofire (Fig. 8a), offline (Fig. 8b) and inline (Fig. 8c). Comparison of the figures shows that each sub-grid-scale fire hotspot more realistically represents the actual high concentration of emission emitted at the source (Fig. 8c) compared to the grid-following averaged-out effect in the offline method (Fig. 8b). Nevertheless, the current setting does not include the two-way aerosol–radiation and aerosol–radiation–cloud feedback. This will be further studied in future work looking at its importance in the cloud-laden SEA region (Tsay et al., 2016), as seen in the missing data due to the cloud cover in Fig. 6d.

4 Transport of biomass burning aerosol to Taiwan

The discussion below is performed using the model output of IWrp + EC and focuses on the high-pollution episodes observed at LABS during 13–28 March 2013, as seen in the gray shaded area of Fig. 3. In the source region of nPSEA, the complex land terrain has played a substantial role in the BB plume lifting. Figure 9 shows the evolution of the PM₁₀ concentration on 13 March 2013 at DAK but over nPSEA through the cross-sectional profile (Fig. 1c). During the day when the fires are active, BB emission is released from the surface (Fig. 9a and b). Along with the rising of planetary boundary layer height (PBLH), the BB aerosol mixes into the entire boundary layer. The residue layer starts to form during the transitional period between the day and night around 17:00 LST (Fig. 9c), when the ground surface cools down. When the atmosphere becomes stable into the night, the aerosol layer remains as the residue layer and does not move down with the boundary layer (Fig. 9d). The plume starts to be advected by the shear of the upper-layer flow at night on the downwind lee side of the hills. The descent of the boundary layer also confines the aerosol and causes a high concentration near the surface. The detachment of the aerosol layer therefore explains the two-layer plume feature from evening into the night in Fig. 4b and c. The dispersion of emission from the pockets is subjected to at least three systems: (i) strong westerlies from Myanmar flowing over the top of valley pockets that confined the emission (terrain structure shown in gray in Fig. 9), (ii) diurnal mountain valley breeze that might trap or disperse the emission, and (iii) local heating caused by the solar cycle that affects the plume rise and disperses the emission. Therefore, the amount of burning emission lifted is greatly coherent with the populated hills along the transport path.

Figure 9 The modeled vertical cross-section profile (BB' in Fig. 1c) up to 5 km over nPSEA on 13 March: PM₁₀ concentration (contour; IWrp + EC; µg m⁻³), horizontal wind profile (vector; m s⁻¹) in x direction and vertical wind profile in y direction (vector; cm s⁻¹), boundary layer height (dotted lines; in meters), terrain (shaded).

Comparing the model output data of the inline (IWrp + EC) and nofire cases, Fig. 3 shows that BB from nPSEA contributes 68±18 % to PM₁₀, 66±18 % to PM_2.5, 41±13 % to O₃ and 58±13 % to CO during the intense BB period (18–27 March) at LABS, while BB contributes 43±31 % to PM₁₀, 41±32 % to PM_2.5, 23±19 % to O₃ and 39±23 % to CO at LABS for the entire month of March 2013. The transport pathway of BB from nPSEA to LABS coincides with the anthropogenic emissions from nPSEA as well as southeastern China; BB aerosols from such emission regions are also captured in the model. Therefore, the actual amount might indicate a slightly lower contribution by BB aerosol than the derived contribution. There are several mechanisms identified in March 2013 to bring BB smoke to Taiwan.

Figure 10 Vertical cross-sectional AA' (Fig. 1b) profile for PM₁₀ (contour), wind in x–z direction (vector), PBLH (dotted lines) and terrain height (gray shade) on 24 March at 20:00 LST (a) with fire and (b) with no fire.

4.1 Westerlies to carry BB emission to LABS

In this case, the BB aerosol lifted is further carried by strong westerlies in the upper layer, around heights between 2–4 km towards LABS. This usually occurred during the night, when the atmospheric boundary layer is low and stable, as shown in Fig. 10. This is the commonly known mechanism that carries the BB plumes to higher ground in Taiwan. This condition occurred on 19–20, 24–25 and 27–28 March 2013. This is the commonly known scenario that is well studied due to the availability of measurements collected at LABS (Lee et al., 2011; Ou-Yang et al., 2014).

4.2 Mixing of BB emission with local pollution on surface

The land surface is heated up, and the boundary layer during the day grows as high as 1.5–2 km a.m.s.l. in western Taiwan, around 1 km a.m.s.l. on the windward side of the central mountain range and up to ∼4 km a.m.s.l. at LABS. When the BB plumes overpass are as low as the PBLH, then the BB aerosol is brought into the boundary layer and mixed to the ground, as shown in Fig. 11. The interaction of BB plumes with local pollutants depends on the loading of local pollutants that are present. The latter are subjected to the local weather system and the occasional Asian continental cold surge that might clean the accumulated pollutants. Such cases usually occur during the morning to noontime, when the land surface heats up, and PBLH develops. This condition occurred on 18, 19, 20, 21 and 28 March 2013. This is the main mechanism where BB aerosol affects western Taiwan. The detection of BB intrusion into surface sites in southwestern Taiwan is not a rare occurrence (Huang et al., 2013; Tsai et al., 2012). It was pointed out that cold surges might be responsible for the downdraft of the BB smoke plumes to the surface (e.g., Lin et al., 2017).

Figure 11 Similar to Fig. 10 but on 19 March 2013 at 08:00, 10:00, 12:00 and 14:00 LST (a) with fire and (b) with no fire.

4.3 Mixing of BB emission with local pollution above surface

Along with the sea–land heat difference, the sea breeze and mountain breeze are formed and enhance the uphill movement of local pollution in western Taiwan. In such a case, the local pollution is brought up to a high elevation to interact with the BB smoke plumes, as shown in Fig. 12. The local pollutants also brought uphill detaches from the planetary boundary layer when the surface cools down quickly. This residue layer of pollutants is then mixed into the BB layers and carried towards the east. Such cases usually occur during midday, when the local pollution plumes have moved up to the hill. This condition occurred on 17, 23 and 25 March. A larger number of fine nanoparticles from local sources is measured at LABS, especially during the morning, even outside of the spring burning season (Chen et al., 2013). Therefore, it is possible that mixing does occur when the local pollutants are transported up the hill through the valley breeze.

Figure 12 Similar to Fig. 10 but on 22 March 2013 at 11:00, 13:00, 16:00 and 18:00 LST (a) with fire (b) with no fire.

Among the three mechanisms, the BB aerosols have a more direct influence on the surface site in western Taiwan under the second mechanism. Such conditions occurred due to an Asian continental cold surge that moved the high-pressure system southeastwards. Under favorable upwind weather conditions, the dust can be lifted and transported downwind and concurrently be present with the BB aerosols. Such a situation is shown by the co-existence of two major pollution events (dust and BB) that reduced the surface O₃, NO_x and SO₂ aerosols over western Taiwan in 2006 (Dong et al., 2018). However, all three of these mechanisms are prone to altering the radiative forcing over western Taiwan. The future incorporations of the aerosol radiative forcing effect through one-way and two-way meteorology–chemistry processes of moisture detrainment and entrainment are necessary to understand the role of BB aerosol in the weather extremes in downwind regions. The cloud–aerosol interaction is particularly crucial to the study of the impact of BB aerosols in cloud-laden regions between nPSEA and Taiwan (Hsu et al., 2003; Tsay et al., 2016). The allocation of smoldering fraction in SEA will need to be improved to account for the tendency of small fires to smolder (Akingunola et al., 2018; Zhou et al., 2018).

5 Conclusions

In this study, several factors involved in the modeling of BB smoke plumes are tested in the WRF–CMAQ model, namely the injection height, initial vertical distribution profile of BB emission, inline PLMRIM and amount of anthropogenic emission. The conventional method used for the study region adopted the fixed-height allocation, which produces an excessive amount of emission over the entire transport route. The initial vertical allocation profile according to the WRAP empirical coefficient (IWrp) improves the surface concentration of the BB emission by the inclusion of the smoldering fraction compared to the default inline PLMRIM (IDef), while replacing the MIX (IWrp) emission in SEA countries with ECLIPSE (IWrp + EC) also improves the pollution concentration simulation at the downwind LABS, especially CO, which is the most important tracer of anthropogenic emission.

The model comparison shows that regardless of the injection height, the main deficiency of the fixed-height offline algorithm originates from its invariant vertical-layer allocation of BB concentration throughout the day. In the complex terrain over the nPSEA region, which is continuous and varies between 0.2 and 1.8 km, mountain meteorology played an important role in the distribution of higher-level plumes. The two-layer structure of the BB plumes observed in the MPLNET extinction coefficient profile at night is well captured by the inline PLMRIM (IWrp + EC), while the offline method (FWrp) gives a time-invariant large value over the entire layers. This highlights that the inline PLMRIM (IWrp + EC) is able to incorporate the diurnal boundary layer physics of the mountain to accurately represent the vertical distribution of the BB concentration in the source and downwind region. It is then clear that the amount of emission produced by the inline reasonably captures the columnar AOD distribution over the transport route between nPSEA and downwind Taiwan when compared to the MODIS columnar product. It is discovered that the inline module with the initial distribution profile of WRAP (IWrp + EC) is able to perform well at both the source and receptor sites compared to the offline module.

The model output shows that the BB plumes near nPSEA are emitted during the day within the BLH. Due to strong mountain valley wind, the smoke plume layers tend to detach from the BLH as residue layers when the surface cools down in the evening to night period. This is the layer of plumes that entered the free troposphere at approximately 1–3 km height and were further transported over to the western North Pacific and Taiwan. The plume layers clearly affect the Taiwan region via three conditions: (a) overpassing western Taiwan and entering the mountain area (LABS), (b) mixing down to western Taiwan, (c) transport of local pollutants upwards and mixing with BB plumes at LABS. The second condition involves the prevailing high-pressure system that is able to impact most of the population in Taiwan and would be an interesting case to explore in subsequent work.

However, care should be taken to select the BB emission inventory input when switching from the offline module to the inline module. The sub-grid-scale allocation of the BB emission requires fitting and testing of the BB emission inventory to make sure it reproduces the individual fires with distinct and realistic peaks. The work highlights the importance of atmospheric-stability-based PLMRIM and the accurate application of emission inventories to capture the plume rise variation at the source site with complex terrain. The correct representation at the nPSEA source site substantially affects the downwind BB concentration at mountain (LABS) and surface sites in Taiwan. It is also observed that the improved setting is able to represent the source site's vertical profile well; however, the height of the plume is reduced following the transport and evolution of the plume approaching Taiwan. This might be caused by the missing algorithm of the indirect and direct effect between aerosols and the high-cloud-cover region along the transport path. It leads to future exploration and incorporation of the effect of cloud–aerosol interaction over the cloud-laden region.