Nepal Ambient Monitoring and Source Testing Experiment (NAMaSTE): emissions of particulate matter and sulfur dioxide from vehicles and brick kilns and their impacts on air quality in the Kathmandu Valley, Nepal

Air pollution is one of the most pressing environmental issues in the Kathmandu Valley, where the capital city of Nepal is located. We estimated emissions from two of the major source types in the valley (vehicles and brick kilns) and analyzed the corresponding impacts on regional air quality. First, we estimated the on-road vehicle emissions in the valley using the International Vehicle Emissions (IVE) model with local emissions factors and the latest available data for vehicle registration. We also identified the locations of the brick kilns in the Kathmandu Valley and developed an emissions inventory for these kilns using emissions factors measured during the Nepal Ambient Monitoring and Source Testing Experiment (NAMaSTE) field campaign in April 2015. Our results indicate that the commonly used global emissions inventory, the Hemispheric Transport of Air Pollution (HTAP_v2.2), underestimates particulate matter emissions from vehicles in the Kathmandu Valley by a factor greater than 100. HTAP_v2.2 does not include the brick sector and we found that our sulfur dioxide (SO₂) emissions estimates from brick kilns are comparable to 70 % of the total SO₂ emissions considered in HTAP_v2.2. Next, we simulated air quality using the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) for April 2015 based on three different emissions scenarios: HTAP only, HTAP with updated vehicle emissions, and HTAP with both updated vehicle and brick kilns emissions. Comparisons between simulated results and observations indicate that the model underestimates observed surface elemental carbon (EC) and SO₂ concentrations under all emissions scenarios. However, our updated estimates of vehicle emissions significantly reduced model bias for EC, while updated emissions from brick kilns improved model performance in simulating SO₂. These results highlight the importance of improving local emissions estimates for air quality modeling. We further find that model overestimation of surface wind leads to underestimated air pollutant concentrations in the Kathmandu Valley. Future work should focus on improving local emissions estimates for other major and underrepresented sources (e.g., crop residue burning and garbage burning) with a high spatial resolution, as well as the model's boundary-layer representation, to capture strong spatial gradients of air pollutant concentrations.

1 Introduction

Air pollution is one of the most pressing environmental issues in South Asia. According to the 2016 Environmental Performance Index, air quality in Nepal ranked fourth worst in the world (Angel and Alisa, 2016) and Kathmandu, its capital and largest metropolis, is one of the most polluted cities in Asia. Kathmandu lies in a bowl-shaped valley with a floor elevation of ∼1300 m surrounded by mountains of 2000 to 2800 m. It is inhabited by approximately 2.5 million people with a steady population growth of 4 % yr⁻¹ (Muzzini and Aparicio, 2013). The primary sources of air pollution within the valley are uncontrolled emissions from vehicles, brick kilns, and biomass and garbage combustion coupled with dust originating from both local fugitive emissions and long-distance transport (Gronskei et al., 1996; Kim et al., 2015; Shakya et al., 2010; Stone et al., 2010, 2012). The unique topography coupled with high emissions of pollutants contribute to low air quality in the valley.

The rapid growth of the vehicle fleet is of particular concern. According to the Department of Transport Management (DoTM) in Nepal, the total number of registered vehicles in the Bagmati Zone (most of which operate in the Kathmandu Valley) increased from 292 697 to 922 831 (12 % yr⁻¹), between 2005 and 2015 (DoTM, 2017). Approximately 80 % of the total registered vehicles are motorcycles, while cars and pickup trucks account for another 13 %. These vehicles emit air pollutants, including carbon monoxide (CO), nitrogen oxides (NO_x = NO + NO₂), non-methane volatile organic compounds (NMVOCs), particulate matter (PM), elemental carbon (EC), and organic carbon (OC). Buses were estimated to emit more than 90 % of OC, EC, and NO_x among all vehicles, while motorcycles were estimated to emit large amounts of CO (50 %) and NMVOCs (66 %) by Shrestha et al. (2013). In 2010, the annual emissions of EC and OC from diesel-powered vehicles in the Kathmandu Valley were estimated at 2117 and 570 t yr⁻¹, respectively (Shrestha et al., 2013). Relative to estimates for the Kathmandu Valley in 2010 based on the global emissions inventory Hemispheric Transport of Air Pollution (HTAP_v2.2; Janssens-Maenhout et al., 2015), the above estimates for vehicle emissions of EC and OC are 80 and 20 times higher, respectively, than those estimated for all sectors combined in HTAP. Considering that Shrestha et al. (2013) did not include emissions from personal cars or trucks in their estimates, vehicle emissions in the Kathmandu Valley appear to be significantly underestimated in HTAP_v2.2.

More than 100 brick kilns of different types throughout the Kathmandu Valley produce over 600 million bricks per year (Quest Forum Pvt. Ltd, 2017; Gronskei et al., 1996; Weyant et al., 2014). The majority of these kilns (97 %) use the fixed-chimney Bull’s trench kiln (FCBTK) technology and its variations such as zigzag kilns. In the Kathmandu Valley, brick kilns typically operate for 6 months a year, generally from December to May. The main fuels burned in brick kilns are biomass and high-sulfur coal (Joshi and Dudani, 2008), both of which emit SO₂ and PM but in differing amounts. Burning biomass emits relatively greater amounts of PM, whereas coal combustion emits relatively greater amounts of SO₂ (Stockwell et al., 2016; Jayarathne et al., 2018). Weyant et al. (2014) estimate the total emissions from the brick industry in South Asia to be 120 Tg yr⁻¹ carbon dioxide, 2.5 Tg yr⁻¹ CO, 0.19 Tg yr⁻¹ PM_2.5 (PM with an aerodynamic diameter less than 2.5 µm), and 0.12 Tg yr⁻¹ EC. Pariyar et al. (2013) report that brick kilns contribute more than 60 % of total SO₂ and PM emissions in the Kathmandu Valley. Using a source apportionment method, Kim et al. (2015) found that brick kilns contribute 40 % of EC concentrations in the Kathmandu Valley in winter. Despite being one of the major sources of air pollution, brick kiln emissions are not included in the existing gridded inventory estimates. The developers of the Regional Emission inventory in ASia version 2 (REAS v2; Kurokawa et al., 2013) discussed that brick kilns are one of the major sources of EC, but their inventory did not include emissions from this sector. HTAP_v2.2 uses REAS v2 for Nepal and thus does not include emissions from brick kilns in their inventory either. Due to the lack of gridded emissions estimates, to date, no studies have explicitly simulated the impacts of brick kiln emissions on regional air quality.

The purpose of this study was to analyze the emissions and air quality impacts due to on-road vehicles and brick kilns in the Kathmandu Valley. Mues et al. (2018) reported that an emissions database is essential to improve the simulation of EC using regional chemical transport models in the Kathmandu region. We first estimated on-road traffic emissions using the latest number of registered vehicles and emissions factors generated for local conditions. We also created a point source emissions inventory for brick kilns using the newly measured emissions factors in the Kathmandu Valley. In the recent Nepal Ambient Monitoring and Source Testing Experiment (NAMaSTE) in April 2015, in situ emissions from several important under-characterized combustion emissions sources were measured, including brick kilns and motorcycles (Jayarathne et al., 2018; Stockwell et al., 2016). Emissions factors obtained from NAMaSTE were used to create the point source emissions inventory of brick kilns. We then modified the HTAP_v2.2 estimates with updated emissions from vehicles and brick kilns. Next, we conducted three simulations using Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) to explore the impacts of emissions from vehicles and brick kilns on the local air quality. These three simulations differed only in emissions scenarios for vehicles and brick kilns; all other model conditions were kept the same.

2 Method

2.1 WRF-Chem model description

We used the regional chemical transport model WRF-Chem version 3.5 in this study. The Regional Atmospheric Chemistry Mechanism (RACM) (Stockwell et al., 1997) is used for gas-phase reactions. Aerosol chemistry is represented by the Model Aerosol Dynamics for Europe with the Secondary Organic Aerosol Model (MADE/SORGAM) (Ackermann et al., 1998; Schell et al., 2001). MADE/SORGAM predicts the mass of several particulate-phase species, including sulfate, ammonium, nitrate, sea salt, dust, EC, OC, and secondary organic aerosols in the three aerosol modes (Aitken, accumulation, and coarse). This aerosol model has been widely used in previous studies (e.g., Gao et al., 2014; Kumar et al., 2012; Saikawa et al., 2011; Tuccella et al., 2012; Zhong et al., 2016). Photolysis rates are based on the Fast-J photolysis scheme (Wild et al., 2000). The Rapid Radiative Transfer Model (RRTM) (Mlawer et al., 1997) accounts for aerosol radiative feedbacks. Lin et al. (1983) and one and a half local Mellor–Yamada–Nakanishi–Niino Level 2.5 (Nakanishi and Niino, 2006) schemes are used to parameterize cloud microphysical and sub-grid processes in the planetary boundary layer (PBL), respectively. The horizontal winds, temperature, and moisture at all vertical levels are nudged to the large-scale meteorological fields from the National Centers for Environmental Prediction (NCEP) Global Forecast System final gridded analysis datasets.

The model domain covered large parts of the Himalayas, India, Nepal, and southwest China (Fig. 1). The domain included three levels of one-way nesting with horizontal grid spacing of 27, 9, and 3 km, each of which was centered on the Kathmandu Valley. The topography of the innermost model domain is complicated, with the Himalayan range sitting across west to east and separating the Indian subcontinent from the Tibetan Plateau. Even when we use 3 km spacing for the nested domain, the model is unable to resolve the very steep topographic features but this was the best we could do with this project, given the resolution of emissions available. There are 31 vertical levels from the surface to 50 mbar. The Model for Ozone and Related Chemical Tracers (MOZART) global chemical transport model (Emmons et al., 2010) was used to provide initial and lateral boundary conditions for chemical species in the outermost domain. We found that the inclusion of dust from MOZART led to overestimated aerosol optical depth (AOD) at Jomsom in Nepal and at Qomolangma (Mt. Everest) station for Atmospheric and Environmental Observation and Research, Chinese Academy of Sciences (QOMS_CAS) in Tibet, China. Therefore, in our simulations, dust concentrations were calculated online, using the Air Force Weather Agency (AFWA) emissions scheme (Marticorena and Bergametti, 1995) with zero initial conditions.

Figure 1 Nested model domains with the terrain heights in meters (color shaded) and the locations of four measurement sites: Jomsom, Nepal (blue asterisk); QOMS_CAS, China (pink dot); Bode, Nepal (red triangle); and the Tribhuvan International Airport, Nepal (black cross).

To analyze the impact of emissions from vehicles and brick kilns on air quality in the Kathmandu Valley, we performed a set of three nested model simulations, referred to as HTAP, HTAP_vehicle, and HTAP_vehicle_brick. The HTAP simulation used the original HTAP_v2.2 emissions inventory as inputs. The HTAP_vehicle simulation used HTAP_v2.2 with updated vehicle emissions (Sect. 2.2.1) as inputs. The HTAP_vehicle_brick simulation used the HTAP_v2.2 with updated vehicle emissions and the additional brick kiln emissions (Sect. 2.2.2) as inputs. The same meteorology and boundary inputs were used for all three nesting simulations as our focus is to understand the impact of emissions on the local air quality. We conducted each simulation for the 2-week period of 12–24 April 2015 during which observational data from the NAMaSTE field campaign were available for comparison (Sect. 2.3). The model was spun up for 5 d preceding the simulation period, which was sufficient to ventilate the regional domain.

2.2 Emissions

2.2.1 Emissions scenarios in 2015

We used three emissions scenarios (Table 1) to investigate the impact of emissions on local air quality in the Kathmandu Valley. The first emissions scenario is the same as the original HTAP_v2.2 (Janssens-Maenhout et al., 2015). HTAP is a gridded global emissions inventory combined with the regional inventories and gap-filled with the Emissions Database for Global Atmospheric Research (EDGAR v4.3) (Janssens-Maenhout et al., 2013). In Asia, HTAP_v2.2 uses the MIX inventory, a regional emissions inventory in Asia, which is also developed based on the “mosaic” approach including multiple existing national inventories (Li et al., 2017). The second emissions scenario utilizes the original HTAP_v2.2 with updated vehicle emissions (Sect. 2.2.2). The third scenario is built on the second scenario and adding emissions from brick kilns (Sect. 2.2.3). We used the latest available HTAP_v2.2 for 2010 as the baseline inventory, as this is the closest year to 2015 that we have the data for. The vehicle and brick kiln emissions were developed for the year 2015.

Table 1 Description for each emissions scenario in 2015.

2.2.2 Emissions from vehicles

Emissions from the road transport sector in the Kathmandu Valley were estimated using the International Vehicle Emissions (IVE) model version 2.0 (Davis et al., 2005). The IVE model is specifically designed to calculate emissions from motor vehicles in developing countries for local conditions and has been used extensively in several countries worldwide, including Nepal (Shrestha et al., 2013), India (Barth et al., 2007), China (Guo et al., 2007; Wang et al., 2008), and Iran (Shahbazi et al., 2016). Emissions were estimated by the product of the base emissions factors, the correction factors, and the distance traveled or the total starts for each of the vehicle categories. We classified vehicles into six overall categories (motorcycles, buses, cars, trucks, taxis, and three-wheelers), numerous technology-based subcategories, and two fuel types (gasoline and diesel). For fuel quality input values, we used unleaded gasoline with a sulfur content of 300 ppm and diesel with a sulfur content of 500 ppm. The base emissions factors are derived from emissions tests conducted mainly in the USA, along with data collected in developing countries. The correction factors consider local conditions (meteorology, altitude, inspection/maintenance program, etc.), fuel quality, and power and driving characteristics (vehicle specific power pattern, road grade, air conditioning usage, and start pattern). Local meteorology data, such as ambient temperature and relative humidity, were obtained from Weather Underground for April 2015. Driving characteristics for motorcycles, buses, taxis, and three-wheelers were adopted from surveys in the Kathmandu Valley, documented by Shrestha et al. (2013). Since we lack survey data for trucks and cars in Kathmandu, we used the data from Pune, India, for these two types of vehicles (Barth et al., 2007). Pune was the only representative city within South Asia, where the International Sustainable Systems Research Center (ISSRC) conducted a detailed study of vehicle activity. The travel distances were obtained from the vehicle registration number and vehicle kilometers traveled (VKTs). Daily VKT and the number of starts per day were adopted from Shrestha et al. (2013), which is specific for vehicles in the Kathmandu Valley. The number of vehicles in 2015 was taken from the government report of the DoTM (2017). Data used for the number of vehicles and VKTs for each vehicle category in 2015 are summarized in Table 2. Table S1 in the Supplement lists the detailed technology of each category, fraction of vehicles with different technology in fleet, and corresponding European vehicle emissions standards (Euro Standards).

Table 2 Number of vehicles, daily mileage, and starts for vehicles in the Bagmati Zone, Nepal, in 2015.

^a We obtained the number of vehicles in 2015 from the Department of Transport, Nepal. ^b VKT of truck/mini-truck are from Malla (2014).

We used the IVE model to estimate emissions of CO, NO_x, NMVOC, SO₂, PM, and some greenhouse gases. All emitted PM was assumed to be PM_2.5 because the ratio of PM_2.5 to PM₁₀ is 0.92 for diesel vehicles and 0.88 for gasoline vehicles in the EPA 2014 MOVES model. Studies such as Gillies et al. (2001) and Handler et al. (2008) have also found that 74 % and 67 % of PM₁₀ is PM_2.5 in on-road studies. Although we understand that assuming all emitted PM₁₀ to be PM_2.5 is potentially an overestimation, we believe that this is acceptable, given the lack of observational data in Nepal or in South Asia.

Because the IVE model does not directly estimate emissions of EC or OC, we used conversion factors derived from the study of Kim Oanh et al. (2010) to estimate these emissions. Kim Oanh et al. (2010) specifically focused on the emissions of diesel vehicles in developing countries and tested a large number of vehicles. For vans, we used an EC ∕ PM mass ratio of 0.46 and OC ∕ PM of 0.2, while for trucks and buses, we used EC ∕ PM of 0.48 and OC ∕ PM of 0.13. We collected a group of these conversion factors from different studies in Table S3. Our EC ∕ PM mass ratio is close to the median value of all the studies listed below. The total vehicle emissions from IVE were distributed spatially and temporally based on HTAP_v2.2. While we acknowledge that using conversion factors from one study ignores the potential uncertainty due to driving pattern, weather conditions, fuel quality, and vehicle characteristics, we also feel that our estimate provides a good middle ground, given the existing study results.

2.2.3 Brick kiln emissions

We identified the kiln types in the Kathmandu Valley by comparing the specific images of brick kiln types with Google Earth images dated April 2015. The monthly emissions of compound i (E_i,j, g month⁻¹) for a certain type of brick j were calculated as the product of the amount of fuel consumed (BK_j, kg-fuel) and the corresponding emissions factor (EF_i,j, g kg-fuel⁻¹):

$$\begin{array}{}\text{(1)}& E_{i,j}={\mathrm{BK}}_j\times {\mathrm{EF}}_{i,j}.\end{array}$$

The calculated emissions were mapped based on the location of each kiln.

The amount of fuel burned for each brick kiln BK_j is estimated using the following formula:

$$\begin{array}{}\text{(2)}& {\mathrm{BK}}_j=P_j\times W_{\mathrm{brick}}\times E_{\mathrm{brick}}/U_{\mathrm{fuel}},\end{array}$$

where P_j is the production of brick kiln j (number of bricks produced per month per kiln), W_brick is the average weight of a brick (kg brick⁻¹), E_brick is the specific energy consumption of a brick (MJ kg⁻¹), and U_fuel is the specific energy density of fuel (MJ kg-fuel⁻¹). Since the production of each brick kiln was not available, we estimated the monthly mean production using one-sixth of the annual average production, as brick kilns in the Kathmandu Valley usually operate 6 months a year. The annual mean production was obtained from a report submitted to the government of Nepal (SMS Environment and Engineering Pvt. Ltd, 2017). Most brick kilns in the Kathmandu Valley are fueled by high-sulfur coal (70 %), which is supplemented with sawdust (24 %), and wood and other fuels (6 %) (Joshi and Dudani, 2008). Considering that the dominant fuel for brick firing is coal and that the EFs used here correspond to emissions from coal-fueled brick kilns, we used the specific energy density of coal to estimate emissions. Values and references for each variable in Eq. (2) are presented in the Supplement (Table S2).

EF values measured from a zigzag kiln during the NAMaSTE field campaign are applied in Eq. (1) to calculate emissions estimates of various trace gases and PM. Table S1b lists the species, associated EF, and references that we estimated the emissions for. We used EFs of a zigzag kiln rather than a clamp kiln for all types of kilns in the valley because zigzag is the most common kiln type identified in the valley. We will provide more accurate emissions estimates in our emissions inventory when EFs of different types of kilns become available.

2.2.4 Other emissions

Anthropogenic emissions of other gaseous pollutants (CO, NO_x, NH₃, SO₂, and NMVOCs) and PM (EC, OC, PM_2.5, and PM₁₀) from major sectors are taken from HTAP_v2.2 for 2010. HTAP_v2.2 is the most recent global emissions inventory that includes emissions from various sectors such as energy, industry, agriculture, residential (including both heating and cooking), aircraft, and shipping and has the highest spatial resolution. However, as mentioned earlier, HTAP does not currently include brick kiln emissions in their estimates and it also excludes large-scale biomass burning and crop residue burning. For nonresidential “open” biomass burning emissions, we therefore used emissions from the Fire INventory from NCAR (FINN) inventory for the year 2015 (Wiedinmyer et al., 2011). For biogenic emissions, we used the Model of Emissions of Gases and Aerosols from Nature (MEGAN) version 2.1 (Guenther et al., 2012). Dust emissions are calculated online, using the AFWA emissions scheme, as described above.

2.3 Observations and statistical methods for comparisons

We compared our simulations with the surface observations of air temperature, relative humidity, and wind speed at two sites in the Kathmandu Valley. The meteorological data at the Bode site at a height of 23 m were collected during the NAMaSTE field campaign and the data at the Tribhuvan International Airport site at standard meteorological monitoring heights (2 m for air temperature and relative humidity, and 10 m for wind speed and direction) were provided by the Department of Hydrology and Meteorology of the Ministry of Population and Environment of the government of Nepal. The Bode site is located in the eastern part of the Kathmandu Valley at a latitude of 27.689^∘ N and longitude of 85.395^∘ E. The altitude is about 1337 m. The airport is approximately 4 km west of Bode (Fig. 2), at approximately the same altitude.

Figure 2 Spatial distribution of brick kilns in Nepal (a) and Kathmandu Valley (b). Red, orange, blue, and green dots denote the fixed-chimney Bull's trench kiln (FCBTK), Hoffmann kiln, vertical shaft brick kiln (VSBK), and zigzag kiln, respectively.

The daily ground-based AOD values at 550 nm were obtained from the AErosol RObotic NETwork (AERONET). We used Level 2.0 for the QOMS_CAS site in China and the Level 1.5 for the Jomsom site in Nepal. In addition, we also compared the space-based AOD values retrieved from the MODerate Resolution Imaging Spectrometer (MODIS) instrument aboard the Terra satellite with the simulated AOD from WRF-Chem. MODIS provides AOD retrievals at a resolution of 10 km ×10 km. In this study, we used Level 2 and Collection 6 aerosol optical thickness at 550 nm. Concentrations of EC and SO₂ at Bode, Kathmandu, were sampled at a height of 20 m during the NAMaSTE field campaign on 12–24 April 2015. We also compared simulated and observed surface SO₂ concentrations at several other sites in the valley (Kiros et al., 2016). The observed surface SO₂ is 8-week mean concentrations between 23 March and 18 May 2013 from Kiros et al. (2016). They were measured at 15 sites in the valley, including five urban sites (Bode, Indrachowk, Maharajgunj, Mangal Bazaar, Suryabinayak), four suburban sites (Bhaisepati, Budhanilkantha, Kirtipur, Lubhu), and six rural sites (Bhimdhunga, Nagarkot, Naikhandi, Nala Pass, Sankhu, Tinpiple) (Kiros et al., 2016).

The overall performance of WRF-Chem in simulating meteorological data and air pollutants against observations was evaluated using the correlation coefficient (r), the normalized mean bias (NMB), the mean fractional bias (MFB), the mean fractional error (MFE), and the root-mean-square error (RMSE). The evaluation is based on 2-week statistics using the daily mean values weighted for the day and night sampling times at each site.

3 Emissions comparison

3.1 Vehicle emissions

The numbers, mileage, and starts for different vehicle types in the Kathmandu Valley, based on the vehicle registration information, are summarized in Table 2. Using the IVE model, we estimated monthly vehicle emissions for CO, SO₂, NO_x, NMVOCs, EC, OC, and PM_2.5 (Table 3). Relative to the vehicle emissions estimates for 2010 by Shrestha et al. (2013), our estimates are about 2 to 4 times higher due to (1) the increases in numbers of vehicles between 2010 and 2015, and (2) the inclusion of diesel trucks that were not considered in the earlier study. The total number of vehicles in this study is about 70 % higher than that of Shrestha et al. (2013). In addition, the estimated running EFs for trucks are the highest among seven categories of vehicles for all air pollutants, with values 4 to 5 times higher than those for buses, which ranked the second highest (Tables S0 and S1). Although trucks account for only 2.3 % of the total numbers of vehicles, they are the major contributor to pollutant emissions due to their substantially higher EFs. Trucks account for more than 80 % of monthly total emissions for both PM_2.5 and NO_x, and 50 % or more for the other pollutants.

Table 3 Total emissions from vehicles, brick kilns in the Kathmandu Valley during April 2015 estimated by this study versus corresponding emissions from vehicles, and all sources considered in HTAP emissions inventory.

Comparison of these new emissions estimates with those from the ground transport sector in the HTAP_v2.2 emissions inventory reveals that CO, NO_x, NMVOCs, EC, OC, and PM_2.5 are significantly underestimated in HTAP (Table 3). For example, the emissions estimates of PM_2.5, OC, and EC calculated using the IVE model were factors of 186, 100, and 375, respectively, greater than those in HTAP. In contrast, SO₂ emissions estimates agreed well and only differed by 17 %. Our revised emissions of PM_2.5, EC, CO, and NO_x from vehicles drive the substantially greater total emissions of these species in the Bagmati Zone relative to those based on the HTAP (Fig. 3).

Figure 3 The monthly mean surface emissions of five pollutants, PM_2.5, EC, CO, NO_x, and SO₂, from all sources in April 2015 used in WRF-Chem for the three simulations. The star indicates the location of the Bode site.

3.2 Brick kiln emissions

We found 112 brick kilns in the Kathmandu Valley, consistent with a previously reported total of 110 (SMS Environment and Engineering Pvt. Ltd, 2017). Figure 2 shows the spatial distribution of these brick kilns. Approximately 40 % of brick kilns are located in the southern portion of the valley, 35 % in the eastern portion, and the rest are in the western portion. We identified four types of brick kilns, including FCBTK, Hoffmann kiln, vertical shaft brick kiln (VSBK), and zigzag kiln. Out of the 112 kilns, the dominant types were zigzag (63) and FCBTK (46), while there were only three Hoffmann kilns and VSBK kilns combined. Based on Eq. (2), the average fuel consumption of these kilns was about 9700 t month⁻¹ or 58 200 t yr⁻¹, which is close to 65 100 t yr⁻¹ estimated in a previous report (SMS Environment and Engineering Pvt. Ltd, 2017).

Table 3 summarizes the estimated monthly total emissions of major air pollutants from brick kilns in the Kathmandu Valley. Of these species, those with the greatest mass emitted were PM_2.5 (135 t month⁻¹), followed by SO₂ (123 t month⁻¹) and CO (98 t month⁻¹), while emissions for other pollutants were less than 20 t month⁻¹. Table 3 also compares emissions from brick kilns and those from all other sectors in the HTAP_v2.2 emissions inventory. Our brick kiln SO₂ and PM_2.5 emissions estimates are each equivalent to 68 % and 16 %, respectively, of the total emissions in the HTAP estimates. The increase in PM_2.5 and SO₂ emissions due to adding the brick kiln emissions can be seen clearly in Fig. 3. For EC, OC, CO, NO_x, and NMVOCs, the brick kiln sector contribution is less than 3 % of our updated total emissions.

4 Model results and evaluation

4.1 Meteorology

The different emissions scenarios did not impact simulated meteorological conditions. Therefore, we only present the statistical analysis of our HTAP_vehicle_brick simulation, using the HTAP inventory with updated emissions for both vehicles and brick kilns. Model-simulated 2 m temperature and relative humidity, as well as 10 m wind speed, are compared to observations at two sites in the Kathmandu Valley: Bode and the Tribhuvan International Airport. Figure 4 shows the comparisons of predicted daily-averaged quantities with observations and Table 4 presents the statistical indices of comparisons for each site. The temperature at Bode is simulated with a correlation of 0.8 and a small negative NMB of 4.7 %. At the airport, the correlation of 0.7 is close to that at Bode, but a larger bias is observed, with NMB of 15.8 %. The model systematically underestimates relative humidity with a correlation of 0.5–0.6 and NMB of −40.8 % to −34.0 %, due to an underestimation of both minima and maxima. In a previous WRF-Chem study in the Kathmandu Valley (Mues et al., 2018), an underestimation of relative humidity was also clearly observed near the ground. The temporal correlation coefficient of daily 10 m wind speed is 0.7 at the airport and 0.8 at Bode. Although the model reproduces the daily variability well, it overestimates the wind speed at both sites. The model performs better at Bode (NMB = 67 %) than at the airport site (NMB = 176 %) in simulating wind speed. The modeled mean wind speed at the airport site is about 1.68 m s⁻¹ higher than the observation, with a RMSE of 1.74 m s⁻¹. WRF-Chem usually has difficulty simulating wind speed over complex mountain terrains; a larger bias over mountain regions was also found in previous studies (Mar et al., 2016; Mues et al., 2018). Zhang et al. (2013) explained that the overestimation in wind speeds is likely caused by poor representation of surface drag exerted by unresolved topographical features in WRF-Chem. A closer look at the wind observational data reveals the presence of a local, thermally driven diurnal wind circulation that controls the airflow regime in the Kathmandu Valley during weak gradient synoptic-scale flow and arguably significantly impacts the air quality in the valley (Fig. 5). During the day, the winds at both Bode and the airport are predominantly from the SW quadrant along the axis of the nearby river, with hourly-averaged magnitudes of up to 5 m s⁻¹. In contrast, the katabatic winds during the night are much weaker, generally under 1 m s⁻¹ and with a prevailing easterly component. Note that winds at Bode are consistently stronger than those at the airport. One of the most likely reasons for that is the differing measurement height at both sites: 23 m above surface at Bode and 10 m at the airport.

Table 4 Statistical performance of model simulation for daily surface temperature, 10 m wind speed, and surface relative humidity at the Tribhuvan International Airport (Airport) and Bode.

Figure 4 Comparisons of observed (blue dots) and modeled (red lines) daily mean 2 m temperature, 10 m wind speed, and 2 m relative humidity at two sites (Airport and Bode) in the Kathmandu Valley.

The model generally reproduces the wind direction shift at both sites quite well (Figs. 6 and 7); however, the wind speed magnitude is substantially overpredicted. The overestimation is not that large during the day but it is severe during the nighttime hours, suggesting serious differences in the structure of simulated and observed nighttime boundary layers. We should note here that the model does not directly predict the wind field at 10 m height; instead it is extrapolated from the first model level using Monin–Obukhov similarity theory (Jimenez et al., 2012), and therefore is highly influenced by the PBL scheme used in the simulations.

Figure 5 Hourly-averaged observed U- and V-wind components during the period of 11–24 April at the (a) Bode site and (b) Tribhuvan International Airport. Daytime winds are shown with red closed triangles and the open blue circles denote winds during the night. Note that the winds at the Bode site and the airport are measured at a height of 17 and 10 m above the surface, respectively.

Figure 6 Hourly-averaged simulated U- and V-wind components during the period of 11–24 April: (a) at 10 m height above surface and (b) at the first model level (28 m height) at the model grid point closest to the Bode site. Daytime winds are shown with red closed triangles and the open blue circles denote winds during the night.

Figure 7 Same as Fig. 6 but for the model grid point closest to the Tribhuvan International Airport.

4.2 AOD

AOD is a column-integrated measurement and, thus, not directly correlated with concentrations of near-surface PM. However, in the context of model validation, AOD is useful as a general indicator of near-surface air quality. Figure 8 depicts the spatial pattern of the 2-week mean AOD observed by MODIS and the modeled AOD at the 550 nm wavelength. In general, the simulated AOD values are higher in the southern part of the domain than those in the northern part. MODIS AOD also shows a similar spatial distribution in the southern part but most data in the northern part are missing due to cloud coverage. It is clear from the figure that adding brick kiln emissions has little impact on AOD in the Kathmandu Valley, while modifying vehicle emissions leads to a significant increase in modeled AOD. The average difference in AOD between the simulations with and without the revised vehicle emissions is 18 % in the Kathmandu Valley. As discussed in Sect. 3.2, diesel engines emit a large amount of EC. These aerosols strongly absorb sunlight at all UV–vis wavelengths and, consequently, contribute to higher AOD values.

Figure 8 The 2-week average AOD: (a) retrieved from the MODIS Terra satellite, (b) simulated using WRF-Chem with original HTAP emissions, (c) simulated using WRF-Chem with HTAP with updated vehicle emissions, and (d) simulated using WRF-Chem with HTAP emissions plus updated vehicles and brick kiln emissions. The star indicates the location of the Bode site.

The time series of simulated versus observed daily mean AOD at the two AERONET sites within the model domain (Jomsom, Nepal and QOMS_CAS, China) and the corresponding performance statistics are presented in Fig. 9. The model tends to overestimate the lower observed values at the QOMS_CAS site, with an MFB of 56 %. At Jomsom, the model predicts the lower measured AOD values reasonably well during 15–24 April. However, the model misses the peak on 14 April, when the observed AOD is near 1.0, an indication of severe air pollution. Instead, our model predicts a somewhat lower peak on the preceding day. This high AOD value was driven by emissions from a wildfire located southeast of Jomsom. It can be seen clearly from the model simulation (Fig. S1) that the fire caused high surface PM and CO concentrations near the burning area on 12 April. Since the prevailing wind direction on 13 April in the simulation was from the southeast and with the overestimated wind speed, the smoke was transported to the northwest and increased the simulated AOD value at the Jomsom site earlier than observed.

Figure 9 Comparisons of observed (blue dots) and modeled (red lines) daily mean AOD at two AERONET sites (QOMS_CAS and Jomsom).

4.3 EC

Figure 10 shows the spatial pattern of simulated EC, averaged during the simulation period. The average simulated EC concentration in the Kathmandu Valley during the 2-week simulation period was approximately 6.2 µg m⁻³, higher than concentrations in the surrounding regions. Vehicles (primarily diesel trucks and buses) contribute approximate 85 % of total EC emissions, whereas brick kilns account for only about 0.11 %. Consequently, the simulation that includes brick kiln emissions (HTAP_vehicle_brick) does not improve the model performance in predicting EC at Bode relative to that with updated vehicle emissions (HTAP_vehicle).

Figure 10 The 2-week average surface EC (a, b, c) and SO₂ (d, e, f) concentrations obtained from three simulations: HTAP, HTAP_vehicle, and HTAP_vehicle_brick. The star indicates the location of the Bode site.

Figure 11 depicts the time series of observed and simulated surface EC concentrations. The average EC concentration observed at Bode during the campaign was 5.6 µg m⁻³ during daytime, 10.82 µg m⁻³ during nighttime, and 8.32 µg m⁻³ for the 24 h average. Two factors contributed to the nighttime increase in surface concentrations: (1) the diminishing mixing-layer depth and (2) the air flow circulation shift in response to surface cooling. As the night progresses, the turbulent mixing in the developing nocturnal boundary layer is suppressed and air pollutants are confined in a shallow layer close to surface. The shift in the wind direction is also conducive to increased surface concentrations as the Bode site is located downwind with respect to the major cluster of brick kilns in the area (Figs. 2 and 5). In contrast, during daytime the mixing-layer depth increases in response to solar heating and promotes the vertical mixing throughout the depth of the whole boundary layer, leading to a notable decrease in surface concentrations. The Bode site during the day is also upwind from the brick kiln cluster, which leads to further reduced concentrations. This idealized scenario holds true for most of the simulation period with the exception of the days between 13 and 16 April, when the local air flow circulation is disrupted by a large-scale disturbance. During that period, the wind speeds are generally below 2 m s⁻¹ even during the day, suggesting that the peak in the surface concentrations (Fig. 11) could be related to suppressed boundary-layer mixing.

Figure 11 Comparisons of observed (blue dots) and modeled EC concentrations in the daytime and nighttime and daily mean for the three scenarios at Bode. Observed values are taken during the NAMaSTE campaign (Jayarathne et al., 2018)

Observed EC concentrations are strongly underestimated with poor correlation when the original HTAP emissions are used as model input (MFB = −125 %, r=0.19 for 24 h averaged EC, Table 5). This is similar to Mues et al. (2018), who reported an underestimation of EC concentrations by a factor of 5 when using HTAP emissions in their simulations. The simulation with updated vehicle emissions (HTAP_vehicle) shows reduced bias and better correlation (MFB = −73 %, r=0.61 for 24 h average), although the model still underestimates EC concentrations.

Table 5 Statistical measures calculated for three model simulations with different emissions inputs for EC. Obs (µg m⁻³) and model (µg m⁻³) are 2-week mean daily average values of observed and modeled EC, respectively. r is correlation coefficient between observation and model simulations; NMB (%) is the normalized mean bias between observations and model simulations; MFB (%) and MFE (%) are the mean fractional bias and mean fractional error; RMSE is the root-mean-square error between observations and the model (µg m⁻³).

The model captures the daytime low concentrations very well during 18–21 April, but it fails to predict the observed peak on 15–16 April. Such high EC concentration episodes can be caused by either stagnant meteorological conditions or enhanced emissions. We examined the precipitation and wind speed at Bode during the high and low episode periods (Table S2). On 16 (high EC) and 19 April (low EC), no rainfall was observed. The daytime wind speed on 16 April (2.4 m s⁻¹) was slightly lower than that on 19 April (2.7 m⁻¹) but the EC concentration on 16 April was 6 times higher. This suggests that a sporadic emissions source (e.g., garbage or biomass burning) that was not accounted for in our model could be responsible for the observed high concentration episode.

The nighttime surface concentrations are significantly underestimated by the model throughout the entire period. One possible cause for this underestimation is illustrated in Fig. S2, showing the diurnal evolution of the modeled mixing layer height. During the night, the simulated boundary layer remains well mixed up to a height of 500 m, which facilitates the vertical transport of pollutants, and consequently, leads to lower surface concentrations.

Figure 12 Comparisons of modeled and observed SO₂ concentrations at 14 sites in the Kathmandu Valley. The modeled SO₂ is the 2-week mean daily SO₂ concentrations averaged from 12 to 24 April 2015. The observed SO₂ is the 8-week mean SO₂ concentrations between 23 March and 18 May 2013 reported in the study of Kiros et al. (2016).

It is also possible that we might have underestimated the EC emissions from the brick kilns near Bode. In our emissions inventory we use an average kiln productivity to estimate the emissions; thus if the productivities of these kilns are higher than the average, or if their efficiency is lower, their emissions could have been underestimated in the inventory. Another possibility is that additional sources are still missing, for example, garbage burning, biomass burning, and/or diesel generators. Stockwell et al. (2016) found that garbage burning in the Kathmandu Valley may produce significantly more EC emissions than previously thought. Due to power outages, especially in the dry season, the use of generators was still prevalent in the valley in 2015. A study conducted by the World Bank (2014) found that nearly 200 000 small power generators, powered by diesel, were used for pervasive power shortages in Nepal. Nevertheless, these comparisons suggest that there is still a need for further improvement of constructing local emissions inventories in the Kathmandu Valley.

4.4  SO₂

Figure 10 presents the 2-week average SO₂ concentrations for the three simulations. The average SO₂ concentration in the HTAP_vehicle_brick simulation is approximately 3.6 µg m⁻³ in the Kathmandu Valley, higher than the surrounding areas in Nepal. Revised vehicle emissions have little impact on SO₂ concentration in the valley, but brick kiln emissions contribute 50 % of simulated SO₂ concentrations. The 2-week mean SO₂ concentration measured at the Bode site was 39.7 µg m⁻³. The model largely underestimates the observation, with the modeled mean SO₂ concentration of 5.3 µg m⁻³. SO₂ is mainly a primary pollutant, directly emitted from sources. SO₂ concentrations are highly related to its emissions. The large discrepancy between the model simulation and observation at this site is probably because our model resolution is not able to capture spatially concentrated high emissions of SO₂ near Bode. Kiros et al. (2016) measured SO₂ concentrations at 15 sites in 8 weeks from March to May 2013 in the Kathmandu Valley and observed high SO₂ concentrations close to brick kilns. In particular, the average SO₂ concentration measured at Bode was the highest (39.2 µg m⁻³) and 2–6 times higher than those at the other 14 sites (4.7–15.8 µg m⁻³) in the Kathmandu Valley. They explained that the elevated surface SO₂ concentration at Bode is mainly caused by nearby brick kilns, which are fueled by coal. We found that there are 12 brick kilns located within the 4 km distance from the Bode site. Although we have included emissions of these brick kilns in our model inputs, we may still have overestimated dilution and/or underestimated emissions from these brick kilns by using monthly average productivity and average emissions factors of zigzag brick kilns. We hope to improve our emissions inventory of brick kilns when more information on individual brick kilns becomes available.

Figure 12 compares modeled 2-week mean SO₂ concentrations with observed 8-week mean SO₂ measured between 23 March and 18 May 2013, reported in the study of Kiros et al. (2016). None of these sites exceeded the Nepal national air quality standard of 70 µg m⁻³ for the 24 h mean, but SO₂ concentrations at the Bode site were almost twice as high as the WHO standard of 20 µg m⁻³. Since our own NAMaSTE campaign only collected SO₂ at the Bode site, we also included the study of Kiros et al. (2016) to illustrate the magnitude difference in observational data at different locations within the Kathmandu Valley. The 2-week mean SO₂ concentration from NAMaSTE in 2015 was 39.7 µg m⁻³ at Bode, while the 8-week mean in 2013 by Kiros et al. (2016) was 39.2 µg m⁻³, showing similarities, giving us confidence that comparing the magnitude difference among sites was possible, despite the difference in observed years. We used these 2013 measurements to represent the ambient SO₂ concentrations during our simulation period. Including brick kiln emissions improves model prediction of SO₂ concentrations at all sites. The simulated SO₂ concentrations with brick kiln emissions are closer to observations compared with those without them, although the model still underestimates SO₂. This underestimation is probably due to brick kiln SO₂ emissions. We applied an emissions factor of 12.7 g kg⁻¹ of fuel measured from zigzag kilns (Stockwell et al., 2016) to all types of brick kilns. These were the only available observational data in Nepal at the time of this study. A more recent study by Nepal et al. (2019) reported that the mean value of the SO₂ emissions factor from zigzag kilns is 24 ± 22 g kg⁻¹ of fuel, which is almost twice as high as that used in our study. If we doubled our SO₂ emissions for brick kilns, the modeled SO₂ concentrations would be much closer to the observations. Assuming the linear relationship in SO₂, the average difference between the observed and modeled SO₂ concentrations would drop from 4.4 to 2.8 µg m⁻³. We plan to revisit our brick kiln emissions inventory as more emissions factors become available. Our study highlights the importance of improving the emissions factor of SO₂ for brick kilns in Nepal. The difference between observation and model simulation ranges from 0.8 to 10 µg m⁻³ for all sites except Bode where the difference is 34.4 µg m⁻³. This result suggests that surface SO₂ concentrations in the Kathmandu Valley are highly variable and are influenced by nearby sources. Future simulations with a higher spatial model resolution and emissions inputs may help to resolve the strong spatial gradients in SO₂ concentrations.

5 Summary and future work

In this paper, we modified the HTAP emissions inventory for Kathmandu's road transport sector. We also developed a point source emissions inventory for brick kilns in the Kathmandu Valley, and examined the impacts of emissions from on-road vehicles and brick kilns on local air quality for April 2015. Emissions from vehicles were updated using the IVE model to reflect the most recent vehicle registration information and the local vehicle technology and driving conditions in the Kathmandu Valley. We found that PM emissions from the road transport sector in the HTAP_v2.2 inventory are largely underestimated. The IVE-estimated EC emissions are 375 times higher than those estimated in HTAP. Our brick kiln emissions estimates were created to account for one of the most important missing sources in the existing emissions inventories. We found that emissions from brick kilns contribute 68 % of the total SO₂ emissions in the Kathmandu Valley in the HTAP_v2.2 inventory.

Using the original HTAP emissions results in large underestimations of both surface EC and SO₂ in the Kathmandu Valley. Our revised vehicle emissions significantly reduced model bias and improved model–observation correlation for surface EC concentrations. We found that surface EC concentrations increased by 50 % on average due to our revised on-road vehicle emissions estimates. Conversely, brick kiln emissions contributed approximately 50 % of the modeled surface SO₂ concentrations in the Kathmandu Valley. Although model performance has been enhanced considerably, by using revised vehicle emissions and by adding newly created brick kiln emissions, the model still underestimates the observed EC by 73 % and SO₂ by 87 % at the Bode site during the simulation period. The large underestimation at Bode could be a result of the site's proximity to large point sources or assuming average EFs for these point sources, but additional sources not included in our inventory could also be important for improving the model performance. More information on the production rates of individual brick kilns and emissions factors for each major type of brick kiln could significantly improve the inventory and comparisons. It is important that the complex topography and meteorology with limited observational data limit the degree of model evaluation currently possible in the Kathmandu Valley. The concentrations of pollutants are highly dependent on the measurement locations and topography of their surroundings and more observational data at a finer scale within the valley are essential to better evaluate the local chemical transport models. Despite the uncertainties, the results here suggest that emissions from brick kilns are substantial and current estimates of emissions underestimate total emissions by omitting this source.

Our main objective was to improve the emissions inventories for on-road vehicles and brick kilns and assess the impacts of revised emissions on local air quality. Our results demonstrate that the existing emissions inventories need significant modification for the road transportation sector. Missing sources in the Kathmandu Valley such as brick kilns are also important in predicting local air quality. We suggest that more efforts are needed to improve local emissions information by updating emissions estimates from major sources and developing an emissions inventory including underrepresented sources, such as crop residue and garbage burning. A more comprehensive and accurate emissions inventory allows the local government to identify and define key emissions sources in the Kathmandu Valley. The improved emissions inventory is urgently needed to robustly evaluate the effectiveness of various future policies on emissions mitigation in this region.