Source profiling of air pollution and its association with acute respiratory infections in the Himalayan-bound region of India

The studies related to air pollutants and their association with human health over the mountainous region are of utmost importance and are sparse especially over the Himalayan region of India. The linkages between various atmospheric variables and clinically validated data have been done using various datasets procured from satellite, model reanalysis, and surface observations during 2013–2017. Aerosol optical depth, air temperature, and wind speed are significantly related (p < 0.001) to the incidence of acute respiratory infections with its peak during winter. Model-derived particulate matter (PM2.5) shows high contributions of black carbon, organic carbon, and sulfate during winter. The wind roses show the passage of winds from the south–west and southern side of the region. Back trajectory density plot along with bivariate polar plot analyses have shown that most of the winds coming from the western side are taking a southward direction before reaching the study area and may be bringing pollutants from the Indo-Gangetic Plain and other surrounding regions. Our study shows that the accumulation of pollutants in the Himalayan valley is owing to the meteorological stability with significant local emissions from burning of biomass and biofuels along with long-range and mid-range transport during the winter season that significantly correlated with the incidence of acute respiratory infections in the region.


Introduction
Air pollutants are strongly associated with morbidity and mortality cases across the world, and accurate assessment of its spatial and temporal variation is important for epidemiological studies to estimate the health effects (Stafoggia et al. 2020). Air pollutants are causing major problems in modern civilization, having harmful toxicological bearings on the environment and human health (Ghorani et al. 2016). Aerosols suspended in the lower atmospheric layer (also termed as particulate matter (PM)) are a serious worldwide concern as they are linked with adverse human health conditions. There have been several epidemiological studies which revealed the association of airborne particulate matter with the incidence of acute and chronic respiratory diseases (Mohanraj and Azeez 2004). Over the last three decades, studies related to the exposure of airborne particulate matter have been documented in epidemiological studies related to human health impacts such as cardiovascular and respiratory diseases and are broadly distinguished into acute or chronic on the temporal framework (Brook et al. 2010;Olstrup et al. 2019). Numerous studies related to air pollutants and human health appraise the role of all types of air pollutants on respiratory blockades commencing from nose and throat irritation to bronchoconstriction and dyspnea (Kampa and Castanas 2008).
Air pollutants such as nitrogen dioxide (NO 2 ), sulfur dioxide (SO 2 ), carbon monoxide (CO), volatile organic compounds (VOCs), heavy metals, and respirable particulate matters are able to diffuse short or long distances and affect different organs of the human body, causing, among others, acute upper and lower respiratory to chronic respiratory tract Responsible Editor: Nicholas Apergis infections among children and adults (Kampa and Castanas 2008). Ambient air pollutants contribute about 8% to the Global Burden of Diseases (GBD) and are increasing continuously with time due to increases in fine particulate matter (particulate matter having an aerodynamic diameter less than or equal to 2.5 μm, i.e., PM2.5) (Cathryn 2017). Air pollutants such as NO 2 , SO 2 , CO, and other VOCs mainly generated from traffic and home burning appliances are very harmful to respiratory health (Faustini et al. 2014). There is a paucity of research related to the health effects of NO 2 exposures in smaller urban areas, and there is reason to believe that the effect of NO 2 varies from urban to non-urban areas (Stafoggia et al. 2020). Respiratory diseases are commonly encountered in medicine and pose a very high burden on the healthcare infrastructure (Victor et al. 2013). Acute respiratory infections (ARIs) are categorized on the basis of anatomic site as upper respiratory infections (rhinitis, sinusitis, rhino pharyngitis, laryngitis, and acute epiglottitis) and with greater severity of acute lower respiratory tract infections which include pneumonias, bronchitis, bronchopneumonia, and acute bronchiolitis and are one of the leading health problems and potent cause of morbidity and mortality among children in developing countries (Santos et al. 2017).
According to the reports of GBD, there are almost 2 million premature deaths in India due to ambient air quality and household air pollution. Thirteen out of 20 cities (across the world) having the highest annual level of PM 2.5 are in India, which lead to 21% increase in mortality (Gordon et al. 2018). Research work carried by India State-Level Disease Burden Initiative CRD Collaborators (Salvi et al. 2018) observed an increasing contribution of chronic respiratory diseases in India over the past quarter century with 10.9% of the total deaths and 6.4% of the total disability-adjusted life year (DALYs) in 2016 compared with 9.6% and 4.5% respectively in 1990, with air pollutants (53.7%) as the major contributing factor, followed by tobacco uses (25.4%) and occupational risks (16.5%). The recent study carried by the India State-Level Disease Burden Initiative Air Pollution Collaborators (Pandey et al. 2021) appraised the impact of air pollution on the health and economy in the states of India and revealed 1.67 million deaths in India are attributed to air pollution that contributes 17.8% of the total deaths in the country, along with the total economic loss of $36.8 billion which constitutes 1.36% of the country's GDP that could hinder the growth of $ 5 trillion economy by 2024.
Besides the variation in the air composition, meteorological factors also play a very significant role in the spreading of contagious airborne diseases and air pollutants. Temperature, humidity, atmospheric pressure, and wind speed are some of the important meteorological variables that determine the transport, dispersion, stability, and growth of air pollutants and airborne pathogens in the ambient atmosphere (Cui et al. 2015& Shrestha et al. 2016. Viruses such as human metapneumovirus (hMPV) and respiratory syncytial virus (RSV) have been found to cause epidemics in the winter season in cold and temperate regions (Darniot et al. 2018). Cold winter months are not only associated with an increased number of deaths but also have a substantial impact on morbidity (Rachel et al. 2014). In the UK, around 25,000 deaths occur annually during the winter months and mortality is considerably higher during the cold winter months (December-March) in the northern hemisphere (Rachel et al. 2014). There is a constant increase in hospitalization and mortality during the winter months, the bulk of which is related to cardiovascular and respiratory diseases (Caroline et al. 2001). Exposure to a severe cold environment has witnessed an increased rate of incidence of respiratory infections (Mourtzoukou and Falagas 2007). The effects of climatic change and air pollution and their impact on human health and the environment have been often discussed during the last two decades (Ariane et al. 2013).
The Himalayan region stretches over 2400 km from westnorthwest (Nanga Parbit 8125 m) to east (Namche Barwa 7755 meters) and is considered to be the roof of the world (Michal 2017). Due to human population burst in the last few decades, both the Himalayan ecosystem and human society are facing odd environmental changes (Prakash et al. 2018). Himalayan valley regions with unique geographical identity observe major seasonal changes in the atmospheric variables that are having adverse effects on respiratory health, mostly encountered due to low socio-economic status, overcrowding, and indoor and outdoor pollutants (Wani et al. 2020). Climate change, pollution, and continuous land use transformation in the Himalayan region have impacts on human health, where ecological changes and economic inequalities further influence the spread of various kinds of diseases (Rita 2019). Long-term trend analysis of atmospheric aerosols over the northwestern Himalayas has clearly witnessed the increasing trend of total suspended particulates (TSP) since January 1996 (Gajananda et al. 2005). Over the past two decades, air pollutants are on the rise and have worsened the regional air quality in the Hindu Kush Himalaya with invasion from Indo-Gangetic Plains (IGP). Biofuels used for cooking and heating such as wood and dung are the leading contributors of CO and PM in the Hindu Kush Himalaya (Saikawa et al. 2019). The unique geographical identity (physical, socio-economic, and cultural environment) of the Himalayan region is linked to different morbidity patterns associated with several respiratory tract infections and has seen an increasing trend in the incidence of ARIs including influenza, pneumonia, and bronchitis while a decreasing trend has been seen in tuberculosis cases from 2013 to 2017 (Wani et al. 2019). These geographical aspects of human health could be due to high variability in climatic conditions and atmospheric stability within the Himalayan region (Wani et al. 2020).
Air quality deteriorates significantly during the winter months in the valley region of Kashmir (Himalayas) with the highest PM 2.5 concentrations (~350 μg/m 3 ), which are about six times the Indian permissible limit (Zainab et al. 2018). Air pollution in the valley region of Kashmir is highest during the winter months due to low temperature and dry conditions along with elevated rise of biofuel emissions (burning of leaves and twigs). These biofuels are mostly a byproduct of agriculture and horticulture residuals and are found to be the potent causes of the rise in pollution level. Back trajectories show the role of westerly winds in the contribution of high PM2.5 levels in the Kashmir Himalayas originated mainly from Afghanistan and other surrounding areas (Zainab et al. 2018). Studies conducted at the Kullu-Manali region of lower Himalaya have shown that vehicular pollution is producing ultrafine particles which are getting readily absorbed by the lungs and are causing serious respiratory and neurological disorders (Sharma et al. 2011).
The mechanisms underlying the seasonal variations are not completely elucidated, but can be surely attributed to the changes in the outside and inside air temperatures, wind chill factors, exposure to sunlight, air pollution, pattern or food intake, and psychological conditions (Caroline et al. 2001). The relationship between the quality of air changes per hour and infection transmission in a closed environment is enigmatic and has a direct impact on the incidences of infectious diseases (Farhad 2013). A study conducted in seven cities situated in the crop-burning area of north India has shown the dependence of various pollutant concentrations on meteorology (Ravindra et al. 2019). The concentrations of pollutants such as carbon monoxide and volatile organic compounds have been found to be directly linked to atmospheric conditions, causing major respiratory damages to the population (Ravindra et al. 2019). Damage to respiratory and cardiac health due to criteria air pollutants have been found in the Delhi and also revealed the possible health damage to other places due to transport of these pollutants to other regions (Maji et al. 2018).
The Himalayan region is a very cold and temperate region which gets even colder in winter season with humid and stable atmospheric conditions (Shrestha et al. 2012). With the transport and buildup of harmful air pollutants, such stable atmospheric conditions contribute to the damage to the respiratory system, making it more susceptible and less defensive against further infections such as acute respiratory infections (ARIs) (Deepa et al. 2017). ARIs include serious respiratory health problems faced mainly by children and the elderly population and are caused mostly by viruses, bacteria, and fungi (Darniot et al. 2018).
In most of the studies done worldwide, ARIs are found to be linked with meteorological factors such as temperature, humidity, and wind speed (Chen et al. 2014& Mäkinen et al. 2009). In India, there have been wide research gaps correlating the incidence of respiratory diseases with the environmental factors and we have not found any study relating respiratory diseases with the meteorological variables and air pollution in the Himalayan region. Therefore, with very rapidly changing climatic conditions and transport of pollutants in the Himalayan region, the effect of air pollution and meteorological conditions needed to be studied with respect to respiratory health.

Study region
The Himalaya lies in the north of South Asia separating the great plains of India from the mighty Tibetan plateau interspersed by many longitudinal valleys. With northward slope movement, the Himalayan range acts as a great climatic barrier between India and Central Asia. Figure 1 shows the map and topographic structure of the study area. The present study pertains to the northwestern part of Kashmir valley (largest valley in the entire Himalayan range), which lies between the Pir Panjal range and Zanskar range of the Himalayas (average altitude ranges between 1038 and 5185 m above sea level), with majority of the population inhabiting the belt at 1038-1600 m above sea level (Wani et al. 2020), with approximately 6967 km 2 estimated to comprise three districts (Kupwara, Baramulla, and Bandipora) (Fig. 1). As per reports of JK Statistical Digest 2016, the study area with a total population of 22.2 lacs constitutes 32.2% of the total population of the Kashmir valley. About 36% of the total horticulture land (hectares) of the whole valley is contributed by this region. The study area has more rural characteristics and poor infrastructure along with prolonged backwardness. The physiographic attributes have contributed the study area popularly known as "gate of winds," thereby owing it great importance (Wani et al. 2020).

Meteorological data
Daily meteorological data (air temperature, atmospheric pressure, wind speed and direction, and rainfall) were provided by the Indian Meteorological Department (IMD Srinagar) for the period of 2013-2017. Planetary boundary layer height (PBLH) which determines the environmental capacity for the dispersion of atmospheric pollutants, was measured by parcel method and by temperature and humidity gradients, i.e., temperature, wind, turbulence quantities, and concentration of pollutants (Shi et al. 2020). Wind chill factor (WCF), an index to appraise the apparent temperature felt on the exposed skin by the combination of air temperature, wind speed, and relative humidity, was calculated by the formula given by National Weather Service, USA (Tanveer 2016). WCF = 35.74 + 0.6215T − 35.75 (V 0.16 ) + 0.4275T (V 0.16 ) where T = air temperature (°C) and V = wind speed (mph). Moreover, we have also used ERA5 reanalysis meteorological data (wind vector at 100 m AGL) to see the regional transport of air masses over the study region. The seasonal wind rose which shows the direction and speed of winds that blow in the region over a period of time was generated using MATLAB.

Atmospheric composition data
Data of monthly air pollution contribution to total PM 2.5 (black carbon, organic carbon, sulfate, dust, and sea salt) for the period of 2013-2017 were retrieved from Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) which is a global reanalysis product that assimilates satellite observations of atmospheric constituents and their interactions with other physical variables (https:// gmao.gsfc.nasa.gov/reanalysis/MERRA-2/data_access/) (Buchard et al. 2017). Aerosol optical depth (AOD) from a moderate-resolution spectroradiometer (MODIS; collection 6; https://ladsweb.modaps.eosdis.nasa.gov/) at 550 nm wavelength (combined data product of terra and aqua satellites) along with AOD obtained from a multi-imaging spectroradiometer (MISR) is used to understand the monthly and yearly aerosol loading in the atmosphere. All AOD products are monthly averaged over the given coordinates from 2013 to 2017. The details about MERRA-2, MODIS, and MISR products and its validity/comparison over the Indian region can be found elsewhere (Buchard et al. 2017;Navinya et al. 2020;Levy et al. 2013& Guleria et al. 2012).

Backward wind trajectory density and bivariate polar plots
Backward wind trajectory air mass analysis helps us in finding a geographically potential source region which could be adding the particulate matter at the receptor site. We have plotted the wind trajectory density on the map which is helping us in specifying the most frequent path taken by the wind through a specific region before reaching the study site. The HYSPLIT model (Draxler and Rolph 2010) backward trajectories were used for the wind density plots shown in Fig. VII (a-d). The HYSPLIT simulations were run for the months of January, April, July, and October, as a representative for winter, spring, summer, and autumn seasons, respectively. Selecting the end day of the respective month, we run two simulations for each day, i.e., at 12 a.m. and 12 p.m., and a total of 50 simulations were run for 25 days in each monthly file. Starting from the study site (34.25°N, 74.39°E), these wind trajectories were traced back to 72 h (3 days) at a height of 500 m. For each month, we downloaded the 5 years (2013-2017) of data and merged them together into a single file of the respective month for further analysis. These files contain the latitude and longitude of the backward trajectories computed on an hourly basis, i.e., each simulation contains 72 entries of latitude and longitude. The entire 5 years of latitude and longitude data (for each month) were used as a function of density on the map (Rolph et al. 2017). Using a color map on MATLAB statistical analysis tool, the number of wind trajectory passing through each point (density of wind trajectories) was specified.
We have also used the latitude and longitude data obtained from the HYSPLIT backward trajectory model for calculating the distance traveled by each wind before arriving at the study site. The end point of each simulation obtained at the 72 nd hour of back trajectory is taken as the source/starting point of each wind trajectory while the latitude and longitude of study site are taken as the end of the wind trajectory. Therefore, knowing the latitude and longitude of the starting and ending points of the wind trajectory, we have calculated the distance traveled by it using the Havesine formula (Prasetya et al. 2020). We obtained the distance for every trajectory and categorized it into three categories, i.e., local transport/LT (trajectories within 100 km distance), mid-range transport/MRT (trajectories within 100-500 km range), and long-range transport/ LRT (trajectories coming from a distance greater than 500 km) (Wangstrom et. al., 2008). Many studies (Rajput et al. 2020;Reizer et al. 2014) have found this Lenschow-type analysis (Lenschow et al. 2001) very useful to identify the pollution sources as LRT/MRT or local based on wind speed and direction.
Five years (2013-2017) daily data of wind speed, wind direction, and AOD is used to create bivariate polar plots for daily mean AOD over the region using ERA5 wind vector data at 100 m AGL. The data is segregated seasonally (winter, spring, summer, autumn) for all these 3 parameters. These bivariate polar plots show the dependence of AOD jointly on wind speed as well as on wind direction (Uria-Tellaetxe and Carslaw 2014). The AOD values are mapped over the wind speed values as a function of different color gradients (specified in the color bar). The wind speed values corresponding to the wind directions are scattered over the polar coordinates with circles of equal sizes. This method could help in clear identification of the direction of the pollutant source. We have used MATLAB software for this analysis. These plots mimicked the wind rose, but with the exhibited AOD values corresponding to wind speed and wind direction, it additionally tells us how much pollutants these winds are carrying with them, therefore clearly specifying the source of pollutants in the study site. Also combining these results with the LP, MRT, and LRT statistics, we can easily interpret how much of these pollutants carried by the winds are contributed from long-range, local, or regional transport.

Socio-economic determinants
To assess the magnitude of air pollutants emitted from different sources such as fuel used for cooking, formation of charcoal from horticulture and crop residues were retrieved from district census handbook 2011(www.censusindia.gov.in) and JK statistical digest 2016 (www.ecostatjk.nic.in), while to assess the transport pressure in the region, annual registration of vehicles was obtained from JK transport 2018 (www.jaktrans.nic.in).

Health data
ARIs including influenza-like illness (ILI) with ICD-10 code including moderate upper to severe acute lower respiratory tract illness with one of the following signs: fever > 38°C, cough, sore throat, abnormal breath sounds, tachypnea, dyspnea, chest pain, etc., which were hospitalized during the period of 2013-2017 in the study area, were taken into account for the study. Monthly data related to acute respiratory infections (ARIs) including influenza-like illness (ILI) were obtained from regional offices of the integrated disease surveillance program (IDSP) of Jammu & Kashmir under the Ministry of Health, government of India (www.dhskashmir. org), for the period of 2013-2017. The IDSP is a decentralized district-based system of weekly surveillance for communicable diseases from all existing health centers in the study area which constitutes 3 district hospitals, 3 sub-district hospitals, 110 primary health centers, and 19 community health centers in the region. The incidence rate defined as the number of new cases in the total population in a specified time (Marlies et al. 2010) was used to measure the frequency of morbidity pattern of ARIs that can be useful in predicting the risk factor of the population to develop a given disease.

Statistical analysis
Pearson's correlation coefficient is used to see the linear relationship between various environmental variables such as air temperature (Avg_T), atmospheric pressure (AP), relative humidity (RH), rainfall (RF), wind speed (Wind_S), wind direction (Wind_D), vapor pressure (VP), wind chill factor (WCF), and aerosol optical depth (AOD) with the incidence rate of ARIs. In this study, we carried regression analysis using both enter and stepwise multiple regression models to see the relationship of each independent variable with the dependent variable and to find the most significant predictors affecting the incidence of ARIs by taking into consideration time factor analysis, i.e., months by assigning dummy variables. A multicollinearity test was carried out to check the degree of relationship among the independent variables using the variance inflation factor (VIF) method. It is important to mention that the value of VIF = 1 indicates no correlation and VIF = 1-5 specifies moderate correlation, and if the value of VIF is > 10 and tolerance < 0.20, this indicates high collinearity between the independent variables, whereas tolerance is the inverse of VIF ranges (0-1) (Shrestha 2020). To check whether the statistical model is a good fit for the data, analysis of variance (ANOVA) F test was performed to test whether variables (independent and dependent) are statistically significant w i t h e a c h o t h e r . R o b u s t s t a n d a r d e r r o r u s i n g heteroscedasticity constant using the HC3 method was estimated to test the precise nature of the sample size taken (Long and Ervin 2013). In order to test the normality distribution of regression analysis data, i.e., impact of different environmental parameters on the incidence of ARIs, the probability plot (PP plot) technique was used. It compares the cumulative distribution function (CDF) of observed values to CDF of expected values. Normality of data can be asserted if the observed values follow the diagonal line (expected values) (Das and Imon 2016). All the statistical analysis was carried in Statistical Package for the Social Sciences (SPSS) 25.v.

Association of aerosol load and meteorological variables with incidence rate of ARIs
Numerous studies pointed out the association of acute and chronic health effects, i.e., respiratory and heart disorders with fine particulate matter (PM 2.5 ), and it has been carved that satellite-retrieved aerosol optical depth (AOD) has a spatial variation with PM 2.5 in some regions and AOD data can be used to estimate the spatial coverage of PM 2.5 and can be utilized in epidemiological studies (Zhiyong 2009). This study region surrounded by the Himalayas has been witnessing erratic climatic and unprecedented weather conditions with changing seasons that are directly affecting the health of the people as seen in the study. We have compared the monthly time series data of incidence rate of acute respiratory infections (IR_ARIs) with different environmental parameters (such as air temperature, wind speed, relative humidity, wind chill, rainfall, and AOD). Meteorological conditions with high air temperature and high humidity tend to cease the growth of influenza viruses and other infectious agents, and therefore, these meteorological variables can be a very useful tool for the spread and surveillance of such viruses (Yonglin et al. 2017). The monthly variations in meteorological parameters and IR_ARIs from 2013 to 2017 are shown in Supplementary information (Table S1). The average air temperatures from 2013 to 2017 were recorded (± 4.18 C) in the winter months, with low wind speed (2.2 mph), low relative humidity (59.4%), and low wind chill (below 2°C). With the beginning of the spring season, the region experiences a sharp increase of air temperature (13°C) (from the average recorded temperature in winter months from 2013 to 2017), wind speed, relative humidity, and wind chill, resulting in the low value of AOD which declines from 0.50 in winter to 0.28 and 0.26 in the spring and summer months, respectively, and again rises to 0.31 in the autumn months. Autumn is characterized by a temperature decline of nearly 8°C from an average recorded air temperature of~22°C in the summer months with atmospheric stability and dry weather conditions (least rainfall) attributing accumulation of air pollutants in the atmosphere. As shown in Fig. 2, the low PBLH which determines the environmental capacity for the dispersion of atmospheric pollutants is indicating the deteriorating air quality during the winter months. A strong association was seen in the variation of different meteorological variables and the incidence of ARIs as shown in Table S1, Fig. 2, and Fig. 3. During the winter months (January, February, November, and December), an extreme dip in the air temperature, relative humidity, low wind speed, and wind chill factor and high value of AOD may have resulted in high cases of ARIs among the population as per the medical records (Fig. 2). During the spring and summer months, cases of ARIs were recorded low as compared to the winter months from 2013 to 2017 and again show a sharp rise in the autumn months in the study region. Figure 3 shows the month-wise variation in the meteorological factors from 2013 to 2017 in the study region.
Apart from the above-discussed environmental variables, socio-economic factors such as education, occupation, income, housing environment, family type, room density, fuel used for cooking (Fig. S1), and smoking behaviors may influence the incidence of acute respiratory infections. In view of these factors, and non-availability of data related to socioeconomic factors and incidence of respiratory diseases, some members of our group had carried a field survey under the medical supervision through semi-structured questionnaire (Wani et al. 2020). A total of 1300 households from 45 sample sites were randomly selected by taking medical blocks as the unit of the study. Wani et al. (2020) have shown that a poor socio-economic structure induces higher vulnerability to respiratory problems in this region due to a weak immune system resulting from exposure to severe climatic conditions, improper nutrition, unhygienic environment, insufficient incomes, discrete occupational patterns, crowded population, low traditional standards of cooking, etc. An illiterate population with least comprehension and understanding showed a high susceptibility of developing respiratory infections as compared to more educated groups.

Statistical modeling
The monthly variations in meteorological parameters and IR_ARIs from 2013 to 2017 are shown in Supplementary information (Table S1). We carried linear regression using both enter and stepwise multiple regression methods. The Pearson's correlation coefficients (Table 1) indicate the relationship of independent variables (environmental parameters) with IR-ARIs. Air temperature (r = − 0.59), wind speed (r = − 0.55), relative humidity (r = − 0.59), vapor pressure (r = − 0.52), and wind chill (r = − 0.53) show negative correlation, Fig. 3 Monthly variations of meteorological parameters (temperature, wind speed, wind chill, rainfall, relative humidity) during 2013-2017 while AOD shows positive relationship with the incidence of ARIs (r = 0.59). No significant relationship is seen between atmospheric pressure (AP), rainfall (RF), wind direction (WD), and months taken as time factor. Full details of Table 1, i.e., significance level of association and total data sample taken, are provided in Table S2 (Supplementary information). A model summary generated by the linear regression test (using the enter method) taking into account the different environmental variables and IR_ARIs is shown in Table S3. R value (0.80) shows the correlation between dependent and independent variables, R square (0.64) shows the total variation for the dependent variable to be explained by the independent variable, and adjusted R square with (0.56) shows the variation of the sample results from the population in multiple regression. All the values generated, i.e., R, R square, and adjusted R square, here indicate good characteristics of the model (Akossou and Palm 2013). The strength of independent variables on the dependent variable using the enter method is shown in Table S4. As shown in Table S4, only AOD has a significant and positive relationship with the IR_ARIs, i.e., the unstandardized coefficient for AOD is 3.83 which indicates that for a one-unit increase in AOD, there is an increase of 3.83 in IR_ARIs per thousand of the population, while other independent variables have not shown any such significance. A multicollinearity diagnostic test (Table S5) shows high collinearity between variables, i.e., Avg_T, WCF, and VP with variance inflation factor (VIF) > 10. This high collinearity between these variables did not affect our results much as these related independent variables were not significantly contributing to the IR_ARIs that have been gauged through stepwise multiple regression analysis.
Stepwise multiple regressions are also carried out to find the most significant variables (independent variables) affecting IR_ARIs. Table 2 shows only 3 independent variables out of the total selected independent variables, i.e., RH, AOD, and wind_S are contributing more significantly to IR_ARIs with a statistically significant p value ≤ 0.05. A model summary of stepwise multiple regression signifies good quality of the model as shown in Table 2. ANOVA test (Table 3) shows all three contributing independent variables, i.e., RH, AOD, and WS, are statistically significant to the incidence of ARIs with p value < .0005. F values > 1 in all significant variables indicate good efficiency of the model. Table 4 shows the strength of independent variables on the dependent variable used (stepwise method). As seen in the table, all the variables have a significant relationship with the IR_ARIs with p value < 0.05. AOD shows a positive relationship, while wind_s and RH show a negative relationship, with the IR_ARIs. Unstandardized coefficient for AOD in stepwise multiple regression is 4.41 that indicates a one-unit increase in AOD; there will be an increase of 4.41 in the IR_ARIs per thousand of the total population. A multicollinearity diagnostic test carried in stepwise multiple regression shows no significant relationship between these significant independent variables as shown in Table S6. All values of VIF are < 10 which signifies no relationship between independent variables (Shrestha 2020).
Robust standard error carried out to check the level of unbiased error in the datasets. Least variation in the values of dataset indicates the precise nature of datasets taken for the study as shown in Table 5. A normal probability chart ( Figure S2; Supplementary information) shows the normal distribution of data, as the observed values (independent and dependent variables) follow the diagonal line of the expected values and indicate the normal distribution of data used in different statistical models. Overall, our statistical analysis shows that incidences of ARIs are significantly correlated with AOD (positive association), relative humidity (negative association), and wind speed (negative association).

Source profiling of local air pollution
Monthly variation in percentage contribution to PM 2.5 (black carbon, organic carbon, sulfate, dust, and sea salt) is shown in Fig. 4. High concentrations of black carbon and organic carbon are seen in the winter months from 2013 to 2017. A high pollution level during these months is mostly attributed to large-scale burning of biomass and biofuels utilized in cooking and charcoal formation. Generally, air quality deteriorates in Kashmir valley with the highest PM 2.5 concentrations (~350 μg/m 3 ) during winter (Zainab et al. 2018). This region shows the highest concentration of black carbon among all the high-altitude observation sites in India and Nepal. A high concentration of black carbon is found in the winter and autumn season in Kashmir Himalayas mainly because of large-scale burning of biomass such as twigs and shedding leaves of plants for the formation of charcoal (Mudasir et al. 2017).
Figure 5a-c shows local sources of air pollutants mainly burning of home appliances (fuel used for cooking), burning of horticulture residuals for the formation of charcoal, and transport sector in the region. The data obtained from the census of India show 78% of the total households in the study area are engaged with the traditional source of energy (firewood, crop residues, and cow dung) for cooking that can be attributed to the dominant horticulture and agriculture sector along with prolonged socio-cultural bindings and economic stagnation prevailing in the area. The dominant horticulture and agriculture sector can be considered to be the major contributing factor for ambient air pollution mostly BC and OC in the region as large scale of burning (leaves and twigs) for the formation of charcoal during the autumn and winter seasons. The study region with valley structure is surrounded by mountains from all sides and traps all the air pollutants mostly emitted from the burning of agricultural/horticulture residuals (with quantum jump seen in existing high-density orchards) and vehicular emissions in the region, along with traditional sources used for cooking. The prolonged backwardness in the region results into a higher dependency of households on fossil fuels and biofuels (firewood, crop residues, dry cow dung, coal, and dry leaves). This is enormously contributing in the emission of BC and OC during winter (Wani et al. 2020). As per the records of JK Statistical Digest 2016, the horticulture sector in the region is contributing 35.8% of the total horticulture land of the entire valley which is attributing to the emission of BC and OC in the lower layer of atmosphere by burning horticulture residues obtained from continuous and unavoidable process (pruning), along with shedding of leaves which is utilized in the formation of charcoal to keep orchard lands clean and protect them from pest and winter calamity. Burning of huge fire woods and dry leaves produces ample smoke mainly in the form BC and OC with least dispersion due to lower PBLH in the winter months. Due to the financial schemes provided by banking sectors, there has been an excessive rise in the vehicular registration in the region that can  be seen as a growing concern for the emission of air pollutants in the atmosphere and making inhabitants of the region susceptible to inhale harmful air pollutants. The winds mainly blow from the southwest, south, and southeastern side throughout the year at 100 m AGL. However, wind rose maps using IMD data at surface (Fig. S1) shows a slightly different prevailing wind pattern, i.e., dominant westerly, southwesterly, and southeasterly in all seasons. These small differences can be understood as altitude variation and location of the IMD measurement site. Zainab et al. (2018) have shown the role of westerly and southwesterly winds in the contribution of high PM 2.5 levels in the Kashmir valley, mainly originated from Afghanistan and other surrounding areas. We have also plotted the number of trajectories passing through a given area through backward trajectory density maps as shown in Fig. 7(a-d). In the month of January (winter proxy), most of the trajectories reaching the study region are coming from the west and south directions of the study area and are becoming densest in this area. The color bar represents the number of trajectories passing through a particular area.

Backward trajectory and bivariate polar plot analysis
We can see here that the number of trajectories is highest southward of the study area, similar to wind rose results. This is a very interesting result as a significant portion of IGP is lying in the southward area where these trajectories are densest. Since IGP is very agriculturally active and highly commercialized and comprises major urban cities, therefore these winds may be bringing a lot of pollutants with them in  Seasonal comparison of LT, MRT, and LRT shows that LRT dominates during winter (Jan~80%), spring (Apr5 5%), and autumn (Oct~58%), whereas MRT and LT together dominate during summer season (Jul~83%). Maximum LT is found during autumn season~14%. However, these results need not to be considered as pollutant load from different transport pathways. It may be possible that dominant LRT may bring more clean air masses, whereas weak MRT and LT may bring more polluted air masses to the study region. In order to delineate these ambiguities, we have analyzed bivariate polar plots, which show the dependence of AOD jointly on wind speed as well as on wind direction. Figure 8 shows the bivariate polar plots of daily mean AOD as a function of wind speed and direction (Carslaw and Beevers 2013) Figure 8(a) shows that higher AOD is associated with lower wind speed coming from the southeast direction, whereas lower AOD is associated with speedy winds coming from the northwest direction during winter. Combining these results with LT, MRT, and LRT indicates that LRT brings cleaner air masses as compared to LT and MRT, and also low wind speed conditions favor the accumulation of pollution during winter. It may lead to an increase in incidences of ARIs in the stagnant stable boundary layer during winter, which is also seen in our correlation analysis. In the autumn season (Fig. 8(c)), a higher AOD is associated with a higher wind speed coming from the northwest direction as well as relatively lower wind speed coming from the east and southeast. It indicates that the LRT and MRT can bring pollution from the northwest part of the IGP region where autumn is characterized by open burning of agriculture crop residues in fields (Kaskaoutis et al. 2014;Sembhi et al. 2020). In summer, high AOD is mainly associated with southerly and southeasterly winds with all wind speed ranges. AOD distribution during summer and autumn is less as compared to that during winter season. The increased PBLH and marine air masses coming via southwest monsoon may be the reason of low pollution load during the summer season. Overall, our backward trajectory density plot analysis with bivariate polar plots indicates that the local pollution from traditional burning activities along with transported pollution coming from the IGP region during the stagnant and stable winter season can increase the incidences of ARIs over the region.

Limitations
The study highlights the role of different meteorological parameters in the incidence of respiratory infections identifying potent causes and sources of different air pollutants in the study region. The study also pointed out the possible contribution of IGP in bringing the air pollutants to the valley which is quite significant. The main limitation of the study worth mentioning here is the non-estimation of the proportion of particulate matter in the atmosphere in different seasons due to non-availability of such data in the region that would have surely boosted this work more. There is also a need to find the source of aerosol particles in the region in order to quantify the contribution of IGP in bringing pollutants to the valley. Also, the correlation of meteorological variables with ARIs is found to be pathogen dependent in some studies. In future studies, we can also include the pathogenspecific casualties on respiratory health and correlate it with different meteorological variables.

Conclusions
The main findings of the study are that seasonal fluctuations in the environmental variables (air temperature, wind, relative humidity, and aerosol optical depth) are strongly associated with the seasonal incidence of ARIs in the Himalayan region.
We found high contributions of black carbon, organic carbon, and sulfate to total PM 2.5 during the autumn and winter months. The source of the pollutants can be attributed to both local and distant sources. Local factors include burning of biomass and biofuels for cooking and other purposes, formation of charcoal from horticulture and crop residues, and a quantum rise in the transport sector in a region with no major industrial setup. A significant amount of air pollutants from distant sources such as the west neighboring countries and the IGP region are brought mainly by long-range transport by the westerly and southwesterly winds. The wind back trajectory density analysis confirmed that winds before entering the valley deflect southward and enter the vast IGP where air mass is laden with highly concentrated pollutants and depart towards the longitudinal depression valley, where they get trapped. Finding the cause and mechanism behind such deflection will be very useful for not just understanding the climate and atmospheric circulation but also for making more effective health policies. With a low PBLH during winter, air pollutants mostly remain near the surface with least dispersion. Inhabitants of the study region get exposed to these harmful air pollutants, as seen through medical records (high cases of ARIs were seen in the winter and autumn months). A strong association was seen between different meteorological variables and incidence of acute respiratory infections in the study, most dominantly in the seasonal variation of air temperature, wind speed, and AOD. The air coming from polluted IGP may also bring pollutants and contribute to respiratory health problems in winter.
The study recommends that the pollution control board of the state must take the urgent initiative to install a particulate matterdetecting analyzer in order to continuously monitor the pollution data that will be benefited to carry more in-depth research and possible conclusion. The government also needs to raise the concern and must promote eco-friendly methods for disposal of horticulture produce (charcoal making, etc.) especially when the growth of horticulture in the entire valley is going to see a quantum jump due to the subsidized introduction of high-density apples in the valley; at the same time, the growth of vehicles is increasing exponentially in the state which adds to the emissions due to the combustion of fossil fuels.
Author contribution MAW and AKM designed the work. MAW, MA, and AKM collected the required data. Data analysis and interpretation were done by AKM, MAV, and SS. AKM, MAW, and IAM drafted the manuscript.
Funding There was no funding available for this work.

Declarations
Ethics approval and consent to participate Not applicable

Consent for publication Not applicable
Competing interests The authors declare no competing interests.