Long-Term Variability of Atmospheric Visual Range (1980–2020) over Diverse Topography of Pakistan

Abstract

A substantial drop in atmospheric visibility or visual range (VR) is documented in several parts of the world. We examined the long-term spatiotemporal variability in ground-based VR data for ten airports in Pakistan from 1980 to 2020. Average VR time series analysis was performed using nonparametric tests of Mann–Kendall (MK), Modified Mann–Kendall (MMK), Sen’s Slope (SS), and Sequential Mann–Kendall (SMK), followed by spatiotemporal mapping. Bad VR (<5 km) conditions prevailed over more than 60% of the study period. Noticeably, VR was better during 1991–2000 due to the high average annual rainfall (905 mm) during the 1990s. Plain areas such as Faisalabad experienced the greatest number of bad VR days with hourly (92.85%), daily (95.35%), monthly (90.97%), and seasonal (90.97%) measurements, followed by Sialkot and Multan. The VR of lowlands in central, southern, and eastern Punjab and northern, central, and eastern Sindh showed decreasing trends. In contrast, the VR in coastal areas of south and southwestern Sindh, in the western highlands, and parts of Khyber Pakhtunkhwa (KPK) province, have indicated positive trends. The results of this study will support policymakers in catering to the issue of declining VR in the region by supporting and developing strategies to limit the factor associated with VR.

1. Introduction

Visibility or visual range (VR) is the maximum distance at which an object can be identified in contrast to the horizon in the background [1,2]. The meteorological description of VR is employed in the Meteorological Aerodrome Report, also known as Meteorological Terminal Aviation Routine Weather Report (METAR), messages and aerodrome (airport) routine meteorological reports [3]. As an essential feature of atmospheric measurement, VR can serve as a visual air quality index [4,5]. The VR of a clean atmosphere may range from 145 km to 225 km. Moreover, in the absence of aerosols under Rayleigh scattering and gas absorption, VR can be up to 300 km [5]. The decline in atmospheric VR depends on the interaction between sunlight and atmospheric aerosols or particulate matter (PM), which can affect absorption or scattering depending on the aerosol’s size distribution, chemical composition, and intensity. The aerosols/PM can be natural, including dust particles, volcanic ash, and particles generated from wildfires. They can be anthropogenic (motor vehicles, energy plants, and biofuel burning) [6,7]. Nearly all PM in the atmosphere shows a little tendency of hygroscopicity, which enhances the PM’s size, volume, and mass [8]. Variations in the natural characteristics of particles, including their size and refractive index, immediately influence the capability of PM to scatter and absorb sunlight and consequently define VR [9]. Due to the rise in fossil fuel utilization in urban areas, PM2.5, sulfur dioxide (SO₂), and nitrogen oxides (NOx) severely impact urban air quality. PM2.5 can have harmful impacts on public health and also reduce atmospheric VR. A large number of aerosol sources depend on the degree of human activities, dust storms, the rainy season, and the burning of fossil fuels [10,11,12]. Furthermore, meteorological parameters like air temperature, wind direction, and speed can influence the spatiotemporal distribution of PM, affecting VR [8].

A substantial drop in VR in several parts of the world has been documented [4,5,13]. Several findings described a negative correlation between VR, temperature, relative humidity (RH), and wind speed [13,14]. Generally, various documents reveal that a decline in VR can occur during wet or dry atmospheric conditions. In dry atmospheric conditions, absorbing and scattering light, natural and anthropogenic resources of trace gases, and aerosols can reduce atmospheric VR [15]. Dense fog, haze, and poor air quality have been described over the East Asian region, Indian subcontinent, and Southeast Asian region. The rise of human-caused emissions due to development in the industrial sector, the high growth rate of urbanization, and the frequent use of fossil fuels for eating and cooking are the main causes of bad VR in these regions [16,17,18,19]. With the help of a statistical model, [20] used VR measurements from 674 meteorological stations in China and examined the long-term variations of PM2.5 pollution from 1957–2014. The author found that the PM2.5 levels changed 60–80 µg/m³ over the northern part of the North China Plain throughout the 1950s–1960s, rising to levels greater than 90 µg/m³. Declining trends of −4.6 µg/m³/year and −12.1 µg/m³/year have been noticed in Xi’ from 1979–2003 and 1979–1996, respectively [21].

Good VR (VR > 5 km) is important for the safety of all types of transport sectors comprising land, air, and sea [22,23]. Furthermore, VR can indicate meteorological and air quality variations and their influence on human health [15]. Poor VR can cause road accidents, flight delays, and visual impairment [10]. The health and safety concerns of extreme pollution events and poor VR have supported the development of VR monitoring networks [23]. The frequently employed VR technique introduced by [24] uses the association between VR and the extinction (EXT) of light due to aerosols. Additionally, empirical models built on the association between AOD and horizontal VR have been created in remote radiation sensing, like the Peterson, 6S, and MODTRAN models. Still, they have not been extensively utilized practically [18,25].

Presently, VR calculations are only accessible by ground-based equipment, like the Automated Surface Observing System (ASOS), and consequently, they require spatial coverage [19,26,27,28]. Several scales, measuring units, and instruments have been used to measure surface VR, but they have some limitations of complicated processing and high cost of operation. For instance, human observation is the most convenient way to estimate VR. Still, this method cannot work in extreme weather conditions (haze, rainfall, fog, etc.) that alter the contrast between the target and the horizon [11]. On the other hand, a transmissometer cannot measure VR at long distances (>50 km), and the nephelometer has some structural design limitations which do not allow the measurement of scattering beyond 170 degrees of scattering angle. Lastly, digital cameras are economical but cannot represent a region fully [29,30].

Similarly, the application of LiDAR for VR estimation has been progressively growing over recent years, and both instrumental advancements and data evaluation methodologies have been increasingly investigated to develop its operations [8,31,32]. For instance, in 1982, Gaumet et al. built a LiDAR system to determine atmospheric VR using the slope technique to reverse the trial data. They compared LiDAR VR calculations with those acquired by transmissometers alongside horizontal and slant pathways [33,34,35,36].

The health and security consequences of severe pollution incidents and decreased VR have encouraged the advancement of entire surface recording networks and innovative satellite sensors [23]. The significant development in satellite remote sensing and its role in meteorology cannot be denied. Still, it cannot be fully applied to VR assessment and estimation because VR is measured horizontally, while satellite remote sensing measures the attributes of the atmosphere vertically. However, these different viewing dimensions open the door to estimating and examining the visual range more effectively over the region. Therefore, we must utilize satellite remote sensing data to estimate surface VR over a large region.

Pakistan’s tourist attractions are neglected because of bad VR and high PM levels. It is believed that urban areas of Pakistan are more polluted than rural areas because of the high level of anthropogenic activities. Despite anthropogenic pollution, dust also decreases VR in some regions of Pakistan. For instance, dust storms approaching the Thar Desert in the pre-monsoon season increase aerosol loading over the area [37,38,39]. High aerosol loadings can decrease VR and result in road traffic accidents (RTA). For example, in 2019, about 36,000 people died in road traffic accidents in Pakistan [40]. The motorway in Punjab is closed for many hours and sometimes even a few days due to winter smog.

In Pakistan, there is insufficient previous research work to address VR trends. The rapid expansion of big cities in the country results in increased events of bad VR. This work is the first attempt to investigate the variability of VR in Pakistan. Considering the diversification in topography, climate and economic development, population growth, urban development, and current air quality conditions in Pakistan, we collected and analyzed the VR trends that were not addressed in the past. We believe that our novel contribution will be beneficial in assessing the climatic dynamics of Pakistan and especially the air pollution status of these stations separately.

In this paper, the introduction section explains the basic concept of VR, its role in meteorology, the phenomenon of air pollution in association with bad VR, general causes and effects of bad VR, the literature review, different methods of VR data evaluations, and the situation of air quality in Pakistan. The data source, site description, trend calculations with Mann–Kendall and Modified Mann–Kendall tests, Sen’s Slope for trend magnitudes, Sequential Mann–Kendall for possible point change, and spatiotemporal mapping are described in the materials and methods section. The results section discusses the outcome of these applied methods with hourly, weekly, daily, monthly, and annual subcategories supported with tables, figures, and maps. The last section is comprehensively written to conclude the outcomes of applied methodologies.

2. Materials and Methods

2.1. Data Source

VR measurements and meteorological parameters (temperature, relative humidity, wind speed, and wind direction) were collected from 10 international airports in Pakistan (Figure 1) installed with an Automated Surface Observing System (ASOS). The ASOS systems are located at airports and serve as the nation’s primary surface weather observing network. The ASOS supports weather forecast activities, aviation operations, and climate research and development. These observations are typically called ‘ASOS’ or sometimes ‘AWOS’ or METAR. The ASOS stations’ primary function is to take minute-by-minute observations and generate basic weather reports. For this article, the temporal resolution of VR data was 30 s, which was later converted into hourly, daily, weekly, monthly, seasonal, annual, and decadal observations within the time frame from 1980–2020. However, many areas of the country have no VR data points, such as Baluchistan. Still, only one station, GWD in Baluchistan, has ground-level VR measurements acquired with the ASOS network.

Table 1. Location, elevation, and data range of selected stations for VR recorded from international airports in Pakistan.

Sr. No.	Station	Latitude (°)	Longitude (°)	Elevation (m)	Date Range
1	Faisalabad	31.37	72.99	185	2012–2020
2	Gwadar	25.23	62.33	29	2005–2020
3	Islamabad	33.62	73.10	508	1980–2020
4	Karachi	24.90	67.13	22	1980–2020
5	Lahore	31.52	74.40	217	1980–2020
6	Multan	30.20	71.42	123	1980–2020
7	Nwabshah	26.25	68.37	38	1980–2020
8	Peshawar	34.02	71.58	360	2011–2020
9	Sukkur	27.72	68.79	59.65	2005–2020
10	Sialkot	32.50	74.53	247	1980–2020

describes the location, elevation, and data range of selected stations for VR study in Pakistan.

2.2. Sites Description

The diverse topography of selected airport weather stations with variations in elevations is a profound blend of landscapes varying from alluvial plains to deserts, broad valleys surrounded by hills, and rivers, plateaus, and coastal areas near the Arabian Sea in the south (Figure 1). The elevation of these selected stations varies from 22 m to more than 500 m above sea level (masl). Among them, the lowest elevated station is Karachi (KHI), with hilly, coastal, and plain topography located in Sindh Province, followed by Gwadar (GWD), the coastal region, with an elevation of 29 masl. The most elevated station among all is Islamabad (ISB), a plateau in nature with an elevation of 508 masl, followed by Peshawar (PSH) at 360 masl, a broad valley surrounded by a plain area and hills. The remaining stations are lowlands with an elevation range of 38–27 masl.

Pakistan’s climate is a continental type of climate, characterized by extreme variations in temperature, both seasonally and daily, because it is located on a great landmass north of the Tropic of Cancer (between latitudes 25° and 36°N) [41]. Five stations (KHI, GWD, NWS, SKR, and MLN) have the hot desert type of climate with an average temperature range from 6.11–44.44 °C, followed by two stations (FSD and LHR) having the semi-arid type of climate with a temperature range of 6.11–40.56 °C. Furthermore, two stations (SKT and ISB) have the humid subtropical type of climate with a temperature range from 2.78–43.89 °C. Finally, there is only one station (PSH) with hot semi-arid climatic conditions where the temperature ranges from 3.33–41.11 °C. Few stations are in Pakistan’s economic and industrial hubs and heavily populated cities, such as SKT and FSD. Some other densely populated stations like LHR and ISB provide service industries and have business centers.

2.3. Applied Methodology

2.3.1. Long-Term VR Variability Analysis (Hourly, Weekly, Daily, Monthly, Seasonal, and Annual)

After collecting VR measurements, the raw data were converted to average values for variability analysis [42,43,44]. The dataset was sub-grouped into average hourly, weekly, daily, monthly, seasonal, and annual trends to gauge variations of VR measurements over Pakistan from 1980 to 2020. Two subdivisions of VR were chosen, ≤5 km and >5–10 km, representing bad and good VR [15]. Percentages of bad and good VR days from raw VR measurements and their summary statistics were computed to understand the trend and extreme values [5,15,39,45].

2.3.2. VR Trend

A Mann–Kendall (MK) trend estimator was used for hourly, weekly, daily, monthly, seasonal, annual, and decadal time series data to understand temporal trends of VR measurements over ten airport weather stations and Pakistan throughout the study period (1980–2020). The Mann–Kendall test was first presented by [46] and then developed by [47], and its application has been recommended by the World Meteorological Organization [48]. MK is suitable for time series that do not follow a specific statistical distribution and can be employed for datasets with missing and irregular sampling intervals [48,49]. Another advantage of this method is that it is not significantly affected by extreme values as observed in some time series. The null hypothesis of this test indicates a random trend in the time series, and the acceptance of the alternative hypothesis indicates a trend in the data series. First, it calculates the difference between observations and applies the sign function to extract the parameter S based on the following equation. For a time series, x_i (i = 1, n), the MK test statistic S is given by Equation (1):

$$\mathrm{S}=\sum _{\mathrm{i}=1}^{\mathrm{N}-1}\sum _{\mathrm{j}=\mathrm{j}+1}^{\mathrm{N}}\mathrm{sign}\left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}\right)$$

(1)

where N is the number of data points, x_i and x_j are the data values in time series i and j (j > i), respectively, and sign (x_j − x_i) is the sign function and is calculated applying Equation (2):

$$\mathrm{sign}\left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}\right)=\left\{\begin{array}{c}+1 \mathrm{if} {\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}>0\\ 0 \mathrm{if} {\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}=0\\ -1 \mathrm{if} {\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}<0\end{array}\right.$$

(2)

Since all the stations have more than ten years, the following relation is given in Equation (3), which is used to calculate the variance [50]:

$$\mathrm{Var}\left(\mathrm{S}\right)=\frac{\mathrm{N}\left(\mathrm{N}-1\right)\left(2\mathrm{N}+5\right)-{\sum }_{\mathrm{i}=1}^{\mathrm{m}}({\mathrm{t}}_{\mathrm{i}}-1)\left(2{\mathrm{t}}_{\mathrm{i}}+5\right) }{18}$$

(3)

In the above equation, m indicates the number of series in which at least one dataset is repeated, and t indicates the frequency of data with equal values. Finally, the Z-score is calculated, following [51,52,53], for the standardized MK test statistics (Z_mk) as given in Equation (4):

$${\mathrm{Z}}_{\mathrm{mk}}=\left\{\begin{array}{cc}\frac{\mathrm{S}-1}{\sqrt{\mathrm{Var}\left(\mathrm{S}\right)}}& \mathrm{if} \mathrm{S}>0\\ 0& \mathrm{if} \mathrm{S}=0\\ \frac{\mathrm{S}+1}{\sqrt{\mathrm{Var}\left(\mathrm{S}\right)}}& \mathrm{if} \mathrm{S}<0\end{array}\right.$$

(4)

In a two-tailed test, the null hypothesis is accepted to find the trend of data series if |Z|≤ Z_α/2 is the significance level and Z_α score is the normal standard distribution at the significant level α. Since the test was two-tailed, α/2 was used. This study used the Mann–Kendall test for a confidence level of 0.95. If the Z-score is positive, the data series trend will increase; if the Z-score is negative, the data series trend will decrease. The strength of monotonicity is assessed using Kendall’s Tau as given by Equation (5) [54]:

$$\mathsf{\tau }=\frac{\mathrm{S}}{\frac{\mathrm{n}\left(\mathrm{n}-1\right)}{2}}$$

(5)

The significance value (p-value) for each trend can be calculated by applying Equation (6):

$$\mathrm{p}=0.5-\mathsf{\phi } \left(\left|{\mathrm{Z}}_{\mathrm{mk}}\right|\right)$$

(6)

where φ () represents the cumulative distribution function of a standard normal variate.

2.3.3. Modified Mann–Kendall (MMK)

Autocorrelation in time series data leads to an increase in the variance of the MK test [55]. If autocorrelation exists in the data, the Modified MK tests such as Hamed and Ramachandra Rao (1998) or pre-whitening the data can be used [15,48]. In this study, the Hamed-and-Rao-modified MK test was used for the data that has autocorrelation effects. The variance Equation (3) is a modified version of Equations (7) and (8) using all significant lags [53]:

$$\mathrm{Var}{(\mathrm{S})}^{*}=\frac{\mathrm{n}}{{\mathrm{n}}_{\mathrm{s}}^{*}}\mathrm{Var}\left(\mathrm{S}\right)$$

(7)

$$\frac{\mathrm{n}}{{\mathrm{n}}_{\mathrm{s}}^{*}}=1+\frac{2}{\mathrm{n}(\mathrm{n}-1)(\mathrm{n}-2)}{\sum }_{\mathrm{k}-1}^{\mathrm{n}-1}(\mathrm{n}-\mathrm{k})(\mathrm{n}-\mathrm{k}-1){\mathrm{r}}_{\mathrm{k}}$$

(8)

where n is the actual number of observations, n_s* is the effective number of observations, n/n_s* is the variance correction factor, and r_k is the autocorrelation function of rank k of the observations at a 95% confidence interval.

2.3.4. Sen’s Slope Method (SS)

The magnitude of the trends is calculated by applying Theil and Sen’s median slope estimator [56]. This method provides an improved estimate of trends because it is not affected by the outliers as written in Equation (9) [48,53]:

$$\mathsf{\beta }=\mathrm{Median} \left(\frac{{\mathrm{Y}}_{\mathrm{i}}-{\mathrm{Y}}_{\mathrm{j}}}{\mathrm{i}-\mathrm{j}}\right) \mathrm{for} \mathrm{i}=1,\dots ,\mathrm{n}$$

(9)

where Y_i and Y_j are data values at times i and j (i > j), respectively.

2.3.5. Sequential Mann–Kendall (SMK) Test

Sequential Mann–Kendall (SMK) was also employed in the last section to examine the change point of VR measurements with time. The SMK test helps locate the year when a trend begins and changes. This test comprises two subsets, progressive series u(t) and backward or retrograde series u’(t), given in Equations (10)–(13). When a progressive series crosses ±1.96 (95% confidence limit) and does not return, it shows a considerable rising or declining trend. The rapid shift can be statistically significant when u(t) and u’(t) crisscross each other away from straight lines at ±1.96 (95% confidence limit). The crossing of u(t) with u’(t) in between the straight lines is regarded as not significant [48]. The sequential value of statistic U(t), which is the forward sequence, was calculated as:

$$\mathrm{U}\left(\mathrm{t}\right)=\frac{{\mathrm{t}}_{\mathrm{j}}-\mathrm{E}\left(\mathrm{t}\right)}{\sqrt{\mathrm{Var}\left({\mathrm{t}}_{\mathrm{j}}\right)}}$$

(10)

$${ \mathrm{t}}_{\mathrm{j}}=\sum _1^{\mathrm{j}}{\mathrm{n}}_{\mathrm{j}} \left(\mathrm{j}=2,3,\dots ,\mathrm{n}\right)$$

(11)

$$\mathrm{E}\left(\mathrm{t}\right)=\frac{\mathrm{n}\left(\mathrm{n}-1\right)}{4}$$

(12)

$$\mathrm{Var}({\mathrm{t}}_{\mathrm{j}})=\frac{\mathrm{j}\left(\mathrm{j}-1\right)\left(2\mathrm{j}+5\right)}{72}$$

(13)

2.3.6. Spatiotemporal Trends of VR variability

Lastly, with the help of the Inverse Distance Weighted interpolation method (IDW) [57], spatial distribution mapping of Mann–Kendall and Sen’s Slope trend estimators was generated to observe spatiotemporal trends over Pakistan [48]. IDW interpolation is a common method of interpolation in spatial analysis. This method uses a linear-weighted combination set of sample points to determine cell values where a greater weight will be assigned to the points closest to the target location. The unknown value, Z(S_o), at point S_o is given in Equation (14):

$$\mathrm{Z}({\mathrm{S}}_{\mathrm{o}})={\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{i}}\mathrm{Z}({\mathrm{S}}_{\mathrm{i}})$$

(14)

where n is the monitoring station, Z(S_i) is the value at the sampled locations S_i, and W_i represents the weight of S_i as shown in Equation (15):

$${\mathrm{W}}_{\mathrm{i}}=\frac{\frac{1}{{\mathrm{d}}_{\mathrm{i}}^{\mathrm{k}}}}{\left({\sum }_{\mathrm{i}=1}^{\mathrm{n}}\frac{1}{{\mathrm{d}}_{\mathrm{i}}^{\mathrm{k}}}\right)}\mathrm{i}=1,2,\cdots ,\mathrm{n}$$

(15)

where d_i is the horizontal distance between the interpolation points and the points observed, and k is the power of the distance.

3. Results

From 1980 to 2020, the total number of data samples was 80,187. The daily VR ranges from 0.05 to 10 km, with an average VR of 4.93 km. There are frequent daily drops in VR below 5 km. The decline in good VR can be due to changing metrological factors associated with wet or dry atmospheric conditions [15]. The overall percentage of hourly bad VR (57.11%) is greater than good VR (42.89%) with the 1st quartile of 3.96 km. The hourly VR ranges from 0.02 to 10 km, with an average VR of 5 km (Figure 2a). Out of ten stations, eight (80%) stations have an hourly average corresponding to <5 km, and only two (20%) stations (ISB and KHI) have an hourly average VR > 5 km. FSD experienced most of the highest bad VR (92.85%), followed by SKT (85.65%). During the night and early morning (1–4 am), FSD has the lowest VR < 2 km, which improves gradually till 8–11 am and then decreases after 11 am. Out of twenty-four, ten hours have experienced <3 km VR, contributing 41.67%, and four hours have <2 km. A maximum VR of 4.18 km is observed at FSD; hence, LHR never experienced good VR throughout the study period. On the other hand, KHI has experienced good VR (81.19%) most of the time. Seventeen hours with VR > 6 km is observed, contributing 70.83%, followed by five hours > 5 km. The highest VR of 7.06 km and 6.96 km is observed at 11 am and 10 am, respectively. Lastly, a minimum VR of 5.01 km is observed at 3 am. Other stations have also observed bad VR frequently, including GWD: 64.36%, ISB: 51.43%, LHR: 63.94%, MLN: 80.06%, PSH: 69.54%, NWS: 65.82%, and SKR: 78.10%.

The evaluation of the weekly time series also shows variations in VR at some stations throughout the week (Monday-Sunday). Overall, all Pakistan (PAK) VR shows 61.44% average bad VR days (38.56% average good VR days), with an average VR < 5 km. PAK VR ranges from 0.05–10 km, with a minimum–maximum average of 4.89–4.98 km. Six stations (FSD, MLN, NWS, PSH, SKR, SKT) have bad average VR (<5 km), while only two stations (ISB and KHI) have 100% good average VR all through the week. Out of ten stations, seven stations (FSD, GWD, MLN, NWS, PSH, SKR, and SKT) have bad VR days during the week, whereas the remaining three stations (ISB, KHI, and LHR) have good VR days (>50% VR days). FSD has experienced 95.33% bad days in seven days with an average VR < 2.50 km. The pronounced change in VR during working days and weekends for FSD is with minimum VR on Monday (1.39 km) followed by Friday (1.74 km). At FSD, maximum VR is observed on Thursday (2.70 km) and Sunday (2.66 km). In contrast, KHI and ISB have 72.52% and 54.01% good VR average VR > 5.5 km during the whole week. SKT, SKR, PSH, MLN, and NWS stations have an almost similar pattern of bad VR variations.

3.1. Daily Time Series Analysis

The daily time series analysis reveals that the average daily VR over Pakistan from 1980 to 2020 has significantly declined, having an average of 4.92 km. The overall percentage of bad VR (61.45%) is greater than good VR (38.55%) throughout the four decades. The minimum average VR is observed at SKR (0.05 km), followed by SKT (0.09 km). The calculations of daily VR indicate a significant decrease at seven (70%) stations (FSD, GWD, MLN, NWS, PSH, SKR, and SKT) that have an average bad VR < 5 km, except for ISB, KHI, and LHR, which have average good VR under the study period. Similarly, eight stations (80%) have more than 50% average bad days throughout the years, except ISB and KHI, which have 46% and 27.48% since 1980. FSD has the lowest daily VR (2.01 km), followed by MLN (3.98), whereas KHI has the highest (6.02 km), followed by ISB (5.85 km).

Daily minimum moving average VR shows that PAK has experienced bad VR most of the time (88.43%) compared to good VR (11.57%). It is declining, with the average rate of change at −0.01 m/year, and good VR is also declining at −0.03 m/year (Figure 2b). After 2010, the 88.07% minimum bad VR moving average was observed and increased to 89.79% in 2016, with an average increasing rate of 1.72% during the last five years. Furthermore, during the last three years, bad VR was 90.65%, and it increased to 100% after 2019, showing a 9.35% increase in bad VR during the last three years. Daily moving average results show that most of the time, bad VR (56.57%) also prevailed throughout the study period compared to good VR (43.43%), which is declining, with the average rate of change at −0.02 m/year (Figure 2c). Moreover, after 2000, 69.79% bad VR days prevailed, which increased to 78.43% in 2015, and 84.86% during the last two years (2019–2020). Bad VR has increased up to 15.07% during the last 20 years. It is also worth noting that only in the last year (2020), 80.59% bad VR days are observed with daily time series moving average.

3.2. Monthly Time Series Analysis

Monthly average VR time series data indicate a significant variation from January-December (Jan-Dec) during the study period in Pakistan (Figure 2d). The minimum VR is observed in Jan (0.05 km), followed by Dec (0.09 km), leading to the lowest average VR of 3.76 km in Dec (3.76 km). The highest proportion of bad VR days is observed in Dec (81.10%), followed by Jan (77.94%). Overall, seven months (March-September) have experienced a good average VR > 5 km, whereas the remaining five months (October- February) are under bad average VR conditions. Furthermore, the percentage of bad VR days (>50%) prevailed in almost all months except April (Apr), which is 47.80%. Minimum VR is observed at SKR in Jan (0.05 km) followed by Dec (0.09 km) at SKT, while the maximum is 10 km observed at 60% of stations, while others (40%) have 1.9 km to 9.99 km. Out of ten, seven stations (FSD, GWD, MLN, NWS, PSH, SKR, and SKT) have experienced <5 km VR during monthly time series assessments compared to ISB, KHI, LHR, which have average VR > 5 km during most of the months. Likewise, 70% stations have >50% bad days except for ISB, KHI, and LHR; among these seven stations, three stations (MLN, SKR, and SKT) have > 80% average bad days. At FSD, monthly shifts of VR measurements are observed from Jan-Dec. Overall, FSD has experienced 90.97% average bad VR days during all months. Minimum VR (0.65 km) is observed in Jan and November (Nov). Almost all the months have experienced bad VR except September (Sep), which has experienced good VR (5.01 km). One can also observe that the minimum VR during Jan increases until Apr, then goes down in May, and again improves until July. It again declined in August (Aug) and moved upward in Sep. After Sep, VR continuously declined until November, and a little improvement was observed in December (Dec). In contrast, KHI has experienced 72.51% average good VR days, with maximum good VR in Sep (6.41 km) followed by May (6.17 km), while minimum good VR occurred in Dec (5.68 km) and in Jan (5.74 km). Lastly, ISB has 54.14% average good VR days.

3.3. Seasonal Time Series Analysis

Seasonal time series variations (autumn, spring, summer, and winter) show 61.53% average bad VR days during all seasons over PAK. Minimum VR is observed in winter (0.05 km), followed by autumn (0.14 km). Spring and summer have experienced >5 km VR (Figure 2e). The highest number of bad days percentage is experienced during winter (74.78%), followed by the autumn season (64.93%). At SKR, minimum VR is observed during the winter (0.05 km) season, followed by SKT in winter (0.09 km). Eight (80%) stations have experienced bad VR during winter and autumn, followed by 60% bad VR in summer, whereas only 50% of stations have observed average bad VR in spring. Similarly, 80–90% of stations have bad VR days (>50%) during winter and autumn, respectively, while 70% of stations have bad VR during spring and summer. Overall, FSD has experienced 90.97% average bad days in all seasons and 100% during summer and winter, followed by SKT and MLN, with more than 80% VR bad days in all seasons. Among all stations, FSD has experienced pronounced variations and minimum VR during all seasons. The lowest bad VR is observed in winter (1.29 km), followed by summer (2.67 km), while the highest is observed in autumn (2.95 km). Comparatively, KHI has experienced 72.48% average good VR days during all seasons, with the highest (74.69%) in summer, followed by autumn (74.51%), and lowest (67.35%) in winter. However, the highest good VR is observed at ISB in spring (6.76 km) and summer (6.48 km), while bad VR is in winter (4.4 km).

3.4. Annual and Decadal Time Series Analysis

Annual time series results reveal bad VR for all of Pakistan. Overall, out of 41 years, 80.49% of average bad days are experienced with <5 km VR compared to good VR (36.59%). A minimum annual average (4.22 km) was observed in 2016 for Pakistan (Figure 2f). FSD has a minimum annual average VR (1.38 km), followed by SKT (1.50 km) in 2015. During all years, FSD, PSH, and MLN have experienced 100% average bad VR days, followed by SKR with 93.75%. Some stations like KHI and GWD have good VR days (75.61% and 61.54%) compared to bad VR days (24.39% and 38.48%). Almost all stations have shown continuous bad VR during the last decade except KHI and PSH, which have good VR after 1992 and 2018, respectively. Lastly, some stations have observed continuous bad VR after the 1990s, including MLN after 1991 and SKT after 1999. However, two stations (GWD and ISB) have experienced a blend of good and bad VR during the study period. FSD is the station that has experienced bad VR with variations throughout the study period (2012–2020). The highest VR was observed in 2014 (4.17 km), followed by 2012 (3.66 km), while the lowest was in 2015 (1.38 km) and 2016 (1.49 km). In contrast, KHI has experienced good VR observations since 1992, and no other station has experienced prolonged good VR annual observations since that time. The highest good VR was experienced in 1999 (7.99 km) and 2001 (7.63 km), while the lowest bad VR was observed in 1991 (4.04 km), followed by 4.08 km in 1985.

The examination of the decadal time series shows a significant decline in good VR throughout the decade from 1980–1990 < 5 km. The overall percentage of bad VR was 70.19% compared to good VR, which was 29.82%. The decadal VR from 1980 to 1989 ranged from 0.19 to 10 km, with an average of 5.05 km. From 1991 to 2000, an improvement in good VR can be observed, which was 89.99% compared to bad VR (40.10%). VR from 1991 to 1999 ranged from 0.19 to 10 km, with an average of 5.80 km. From 2001 to 2010, bad VR again prevailed with a percentage of 63.80, whereas good VR was 36.19%. During these ten years, VR ranged from 0.19 to 10 km, averaging 4.69 km. The situation of bad VR was worse from 2011 to 2020. During this decade, bad VR was 68.81% compared to good VR, which was 31.18%. VR in these years ranged from 0.05 to 10 km, with an average of 4.45 km.

3.5. VR Monotonic Trend Evaluations

We performed the Mann–Kendall (MK) test for trend analysis and Modified Mann–Kendall (MMK) test for serially correlated data using the Hamed and Ramachandra Rao (1998) variance correction approach. Overall results of both tests show decreasing trends of VR all over Pakistan during the last four decades (Figure 3). Sen’s Slope (SS) statistics also indicate decreasing trend magnitude (km/year) at most stations throughout the study period. However, some stations have significant decreasing trends, some have significant increasing trends, and some have no monotonic trends. The detailed MK, MMK, and SS results are listed in Tables S1–S4.

MK and MMK statistics on hourly average VR data and SS magnitude show that out of ten stations, the trends of three stations (FSD, GWD, and SKT) are not significant with the monotonic trend. Seven stations (63.63%) have significantly average decreasing monotonic trends (km/year). At all stations, the SS magnitude of significant trends is almost the same (−0.001 km/year). Nine stations (81.81% stations) show a significant decreasing magnitude, one has no trend, and only one (FSD) has a positive magnitude value (0.021 km/year), which is also statistically insignificant at a 95% confidence level. A significant positive trend was not observed at any station, including Pakistan, with hourly average VR data.

The average weekly VR data results reveal significant and insignificant trends. Three (30%) stations (FSD, GWD, ISB) have insignificant increasing monotonic trends (no trend), and only one station (SKT) has an insignificant decreasing trend. In contrast, only 20% (only 2) stations (KHI and PSH) have a significant increasing trend. On the other hand, 50% stations (LHR, MLN, NWS, SKR) and PAK have significant decreasing trends. Average SS results show that only two stations (SKR and LHR) have the highest significant negative magnitude (−0.003, −0.002 km/year), and only one station (PSH) has a significant positive magnitude (0.001 km/year). The remaining stations either show static magnitude (0.000 km/year) with average weekly VR measurements or their magnitudes are insignificant. Five stations (KHI, LHR, MLN, PSH, and SKR) have significant increasing or decreasing trends during seven days of the week, while some have a blend of significant increasing or decreasing trends. For instance, ISB has significantly increasing trends on Tuesday, Wednesday, and Thursday. NWS has significant decreasing trends on Monday, Tuesday, and Thursday, while SKT has significant decreasing trends only on Tuesday. FSD and GWD did not experience any significant trend throughout the week except on Friday; only FSD had an increasing trend. The magnitudes of these significant trends are also significant. Neutral SS (0.000 km/year) is observed at KHI and MLN stations throughout the week. ISB and NWS also have neutral SS magnitude for significant trends. FSD and GWD do not have significant magnitudes in all seven days except SKT, which only has a significant neutral magnitude on Tuesday.

As per average daily VR data, four stations (LHR, MLN, NWS, and SKR) have significant decreasing trends, which shows that during the last 41 years, 40% of stations have decreased VR trends, while four stations (FSD, ISB, KHI, and PSH) have significant positive trends. Only two stations (GWD and SKT) have insignificant positive trends. The average trend magnitudes based on the average daily VR ranges between −0.001 km/year and 0.012 km/year. The only significant decreasing magnitude is observed at LHR, while the only significant positive increasing magnitude is observed at FSD. However, five (50%) stations have neutral SS values, and two stations have insignificant SS magnitude (GWD and SKT), contributing only 20%.

Monthly average VR shows decreasing trends (Jan–Dec) for Pakistan and five stations (GWD, LHR, MLN, NWS, and SKR). However, except for LHR, SKR, and PAK trends, MK, MMK, and SS results are insignificant at eight stations (FSD, GWD, ISB, KHI, MLN, NWS, PSH, and SKT), which means that these stations have no monotonic trend. The highest significant decreasing trend is observed at LHR (−0.361), followed by SKR (−0.353). MMK results with GWD data show an insignificant positive monotonic trend with monthly measurements compared to MK, which shows an insignificant decreasing trend. Overall, 50% of stations and PAK have insignificant positive upward VR trends. The magnitude at eight stations has positive and insignificant negative values per km/year, whereas only two stations (LHR and SKR), including PAK, have negative significant declining SS.

The VR measurements with the seasonal average trends reveal significant negative trends for most stations, except SKT, which has insignificant decreasing trends. The highest significant declining trend is observed at LHR (−0.113 km/year), followed by SKR (−0.093 km/year). Only three stations (FSD, KHI, and PSH) have positive VR trends; in PSH, trends are insignificant. The SS magnitudes of only two stations (GWD and PSH) have positive insignificant magnitude, while other stations (KHI, MLN, NWS, SKR), including PAK, have significant neutral magnitudes (0.000 km/year).

MK and MMK analysis with annual average VR reveals that only three stations (FSD, GWD, and NWS) have insignificant positive upward trends. In contrast, six stations (ISB, LHR, KHI, MLN, SKR, and SKT) have significant downward trends, except for PSH, which has significant declining trends. Significant negative magnitude is observed at PSH (−0.002 km/year) and PAK (−0.001 km/year), while all other stations have insignificant SS values. One can also observe that MMK shows different results as compared to MK. The reason for this difference is that the MKK algorithm is developed for serially correlated data which uses a variance correction approach. In MKK, data are initially detrended, and the effective sample size is calculated using significant serial correlation coefficients [55].

We noticed significant decadal trends pointing to possible teleconnections between VR and climatic drivers such as IPO and PDO, as documented by [58,59,60]. However, further investigations on teleconnections with VR and climatic variables are beyond the scope of this paper. Analysis with 1980–1990 data shows significantly increasing trends with neutral magnitude values (0.000 km/year) for all PAK, while the remaining three decades have decreased significantly, meaning that VR trends have declined since 1991. The highest significant decreasing trends (−0.03 km/year) and magnitude (−0.002 km/year) are observed during the last decade.

Table 2. Trend statistics for decadal timescale visibility over Pakistan. Bold values in the table indicate positive trends; italic values are decreasing trends; neutral values have no trend; PAK = Pakistan, MK = Mann–Kendall, MMK = Modified Mann–Kendall, SS = Sen’s Slope.

Station	Test	1980–1990	1991–2000	2000–2009	2010–2020
PAK	MK	0.04	−0.02	−0.02	−0.03
	SS	0.000	−0.000	−0.001	−0.002
	MMK	0.04	−0.02	−0.02	−0.03
	SS	0.000	−0.000	−0.001	−0.002

lists the MK, MMK, and SS values for all four decades.

3.6. Spatiotemporal VR and Trend Mapping

We applied the IDW interpolation method to get raster data from point station data to observe the spatiotemporal distribution of VR measurements for daily and seasonal time series (Figure 4). The same was done to compute trends for average hourly, weekly, daily, monthly, seasonal, and annual VR for all of Pakistan (Figure 5).

The spatial distribution of daily VR mapping unveils that almost all parts of Pakistan have experienced bad VR (<5 km) except ISB and some southern and southwestern parts of Sindh and southeastern parts of Baluchistan (Figure 4a). Central parts of Punjab, including FSD, have experienced the lowest VR (2.08–2.98 km), followed by the maximum part of Punjab, including MLN and SKT stations. Furthermore, both Khyber-Pakhtunkhwa (KPK) and Baluchistan, excluding the southeastern parts of both provinces, have observed daily average bad VR during the study period.

In autumn (SON), bad VR is observed over central, eastern, and southern Punjab along with central, western, and southwestern Baluchistan and KPK. Good VR is observed in southwestern and southern parts of Sindh, including KHI station and ISB (Figure 4b). The lowest average autumn VR is observed at FSD and its adjacent areas, while ISB and KHI have experienced good VR this season.

In spring (MAM), bad VR was observed over the maximum parts of Punjab, Sindh, and Baluchistan. These regions have experienced bad VR < 5 km, whereas the southern and western Sindh areas, including KHI station, southern, southeastern, and southwestern parts of Baluchistan, and some central and northern parts of KPK, have good VR (Figure 4c). Moreover, the areas of good VR can also be noted in ISB, Gilgit Baltistan (GB), and Azad Jammu and Kashmir (AJK). The lowest VR is observed in FSD, followed by MLN and SKR regions.

In summer (JJA), bad VR is observed over most parts of Punjab, KPK and Baluchistan, in the northern and eastern parts of GB, and central, northern, northwestern, and eastern Sindh. The lowest VR is observed in FSD, followed by MLN (Figure 4d). The area of good VR is observed in central and southwestern parts of GB, in AJK, in LHR, and in ISB, central, southern, and western parts of Sindh, including KHI and southeastern parts of Baluchistan.

In winter (DJF), bad VR is observed almost all over Pakistan except southwestern parts of Baluchistan, including GWD, and southwestern parts of Sindh, including KHI. It is worth noting that GWD and KHI are also affected by desert dust brought from the Arabian Peninsula [38]. The lowest VR is observed in FSD, followed by MLN (Figure 4e).

Overall, hourly MK VR trend mapping shows that most of Pakistan, including most parts of Sindh, KPK, Baluchistan, and northern and southern Punjab, have decreasing MK values. Hence, VR decreased during the last 41 years from −0.06 to −0.35 km/year. The Sistan basin in Baluchistan is a major world source of desert dust [59] which is commonly transported over all parts of Baluchistan [60]. Only a few areas in Punjab and Sindh have shown increasing trends (Figure 5a). SS values also show the decreasing magnitude of VR measurements in many parts of KPK, Baluchistan, Sindh, and Punjab (−0.001 to −0.322 km/year) except the central and northern Punjab, which is showing a positive increasing magnitude of VR (Figure 5b).

Many parts of the country have experienced decreasing trends with average weekly VR data, including parts of Punjab and Sindh (Figure 5c). Similar decreasing and increasing magnitudes are also observed at the same stations. However, some regions of the country have neutral SS values, meaning no increase or decrease in magnitude (Figure 5d).

Daily trend mapping also shows that central, southern, and eastern Punjab and upper parts of Sindh have experienced decreasing VR trends. Other parts of the country have upward VR trends (Figure 5e). Almost all of the country has neutral SS values except some central parts of Punjab that have positive magnitude (0.012 km/year) (Figure 5f).

Monthly trend mapping shows a blend of decreasing and increasing VR trends. However, these values are significant or non-significant. Maximum parts of Punjab, Sindh, and Baluchistan have decreasing magnitude compared to KPK, northern areas, and southern parts of Sindh, which have increasing VR trends (Figure 5g). Monthly SS values also have a mixture of positive and negative magnitudes, which may be significant or non-significant. Almost half of the country (Punjab, Sindh, upper parts of Baluchistan, and lower parts of KPK) has decreasing magnitude. In comparison, the other half (central and lower parts of Baluchistan, northern areas, along with central and upper parts of KPK) has an increasing magnitude (Figure 5h).

Seasonal trend mapping shows decreasing trends in most of Punjab, Sindh, and Baluchistan, while KPK, southern Sindh, and northern areas have positive magnitude (Figure 5i). SS magnitude in seasonal mapping shows that maximum parts of the country either have negative magnitudes or neutral values, except central and southern parts of Baluchistan, which have positive SS values (0.085 km/year) (Figure 5j).

Lastly, annual trend mapping also has a blend of decreasing and increasing MK values which may be significant or insignificant. Many parts of Punjab and KPK have negative VR trends compared to parts of Baluchistan and southern Sindh, which have positive trends (Figure 5k). SS magnitude in annual mapping shows that the maximum parts of the country have negative magnitudes except for central and southern parts of Baluchistan, which have positive SS values. No neutral magnitude is experienced in any part of the country with annual trend mapping (Figure 5l).

3.7. Possible Change Point Detection Analysis

Figure 6 illustrates the monthly SMK for the PAK VR data. The monthly VR time series results show that all months have at least one change point. These monthly change points may be significant or non-significant. Moreover, the results also show that VR measurements have significantly decreased over many months in the last few years. For instance, there was a significant VR decline in Jan after the year of 2016. Mar (2013, 2014), May (2005), and Dec (2007) also have significant change points. Likewise, in February (Feb), a significant VR increase was observed after the 1990s, while a significant decrease was observed after 2000 and later in 2011. During the 1990s, most months had increased levels of VR, which declined after 2000. Mar and Apr show a significant decline during the last 5–6 years after 2014 and 2015, respectively. May and June have similar upward and downward VR measurements. A small increase in VR was also observed in July (2005). Aug, Sep, and Oct also showed significant upward VR from 1995 till 2000; after 2010, these values declined significantly. Nov has a significant slight change before and after 1985. One can also observe that some months have various insignificant change points along with insignificant increasing and decreasing VR measurements. Jan did have one insignificant change point (abrupt change) in 2006, and no significant increase was observed during Jan. The last two months (Nov and Dec) do not show any significant decrease before 2006. All other months also have insignificant change points. The detailed seasonal and annual SMK results are illustrated in Figures S1 and S2.

4. Discussion

There is a pronounced change in VR for hourly and weekly time scales at FSD due to the operational activities of factories with industrial emissions, SO₂, NO_x, and fine particles [61]. These emissions are aggravated due to emissions related to the daily commute of people between home and work, where 61.31% of the total vehicles in Punjab are in FSD [62]. Due to the arid type of climatic conditions, natural and anthropogenic activities such as high temperatures, less rainfall, suspended dust particles, transboundary air pollution, air masses from desert areas (Thar, Thal, Cholistan, Kharan, and Nara deserts), fossil fuel burning, vehicular and industrial emissions at MLN, NWS, and SKR stations that degrade their regional air quality level ultimately degrade VR [12,63,64,65,66,67,68,69,70,71,72]. On the other hand, KHI and ISB have good VR most of the time. KHI is a coastal city and a busy port that periodically encounters the sea and land breeze, specifically from the Arabian Sea and the Bay of Bengal. Coastal wind sweep contributes to decreasing levels of lower-tropospheric pollution. ISB is a plateau in nature with an elevation of 508 masl and a humid subtropical climate with an average VR of 5.75 km (Figure 3a). Central and northern parts of the country have humid characteristics due to highlands and forest cover, and humid and subhumid regions have longer wet spells compared to arid regions, which results in low air pollution and good VR compared to lowland regions [73,74].

The minimum VR rate of decline is slower than the maximum VR, as depicted in Figure 2b. The gap between maximum and minimum VR decreased rapidly (Figure 2b) in 41 years, which may approach zero in a couple of decades, resulting in absolutely diurnal change (Figure 2a).

KHI has experienced good VR because wind speed in KHI ranges from 8–14 mph, whereas in FSD, it is 4–6 mph [75]. This results in the diffusion of aerosols and dust particles in KHI more frequently than in FSD. Seasonal VR variations in FSD can be linked to monthly variations in production phases in FSD and seasonal shifts to complete seasonal production in advance [76,77]. Although there are many cottage and automobile industries in KHI, the city is also known for its software-outsourcing center, with an annual growth rate of nearly 6.5 percent, resulting in a lower industrial emission rate than FSD [71].

All through the country, minimum VR is observed in winter (0.05 km), with the highest number of bad days (74.78%), followed by autumn (0.14 km), with 64.93%. Almost all stations have average bad VR during winter. In winter, highly populated areas in Punjab (FSD, MLN, SKT), Sindh (SKR), and nearby villages use biomass burning as a major heating source. Additionally, large-scale crop residue burning in the presence of winter-induced low mixing height in stable air masses limits the particulate diffusion. It formulates secondary particulates to reduce VR [78,79]. However, KHI and GWD experience good VR in winter, while KHI and ISB experience good VR in autumn. The higher concentration of aerosols in cold air results in frequent severe smog events in LHR, FSD, and some other areas of Punjab. The low boundary layer in cold temperatures raises near-surface pollutants’ concentration [80,81]. The reasons for good VR in KHI in different seasons are high seasonal wind speed in each season compared with FSD, rainfall frequencies, and moderate temperature [63,72].

Six stations (FSD, SKT, SKR, NWS, MLT, and GWD) experienced bad VR during spring. This can be linked to dust storms during this season. Agricultural pre- and post-harvest burning and anthropogenic activities in urban and rural areas may be responsible for the high pollution over the region [68,71]. The high pollution levels and bad VR over GWD are due to severe long-range dust storms in spring and summer. The source regions of these storms are in Iraq, and they enter from the southwest corridor of Pakistan. The other reasons may be low rainfall, high temperature, hygroscopic development of aerosols, biomass-burning activities, and crop harvesting events [82,83]. The reasons for poor air quality and bad VR at most stations in summer and autumn are high temperature, wind speeds, absence of rain, and dust storm frequencies. The airflow in summer blows moist air, strong sea breeze, and moderate humidity levels in KHI. In autumn, a decline in atmospheric pollution (CO, O₃, etc.) is primarily caused by stronger wind speeds, rainfall, moderate temperatures, and cloud cover [77,78,79,80].

Out of 41 years, good VR prevails for 15 years, and all good VR is observed before 2002, which is only 20% of the total study period. FSD, PSH, and MLN have experienced 100% average bad VR days, followed by SKR with 93.75%, while KHI and GWD have more good VR days (75.61% and 61.54%) compared to bad VR days (24.39% and 38.48%). In 1980, SKT, ISB, and LHR had good VR, while KHI, NWS, and MLT experienced bad VR, as shown in Figure 2f.

Significant decadal declining trends of VR were observed, which may have possible teleconnections with IPO and PDO, as documented previously [59,60]. The possible reasons for bad VR and poor air quality during the 1980s are the consumption of oil and petroleum products in the industrial sector, low-quality fuel in motor vehicles and the energy sector, and consumption of biomass for cooking and heating purposes in rural areas of Punjab and Sindh. Pakistan had its wettest decade in the 1990s [63,84,85,86,87]. A decade-by-decade analysis shows that, on average, annual rainfall in Pakistan was 905 mm during the 1990s, which is 22% higher than the next rainiest decade of the 1980s, when average rainfall was 793 mm. The country’s driest decades were the 1960s and 2010s, when annual average rainfall stood at 773 mm and 797 mm, respectively, during both these periods. The possible causes of bad VR from 2000–2010 were high urban population growth rate, energy crisis at the country level, utilization of fossil fuels in the thermal energy sector and for domestic purposes, anthropogenic activities on a massive scale, rapid increase in private vehicles, widespread industrial development, and poor management of air quality monitoring organizations [61,63,66,88].

In addition to the reasons above, the rapid increase in bad VR of FSD, MLN, SKT, and PSH can be due to an increase in poor air quality during 2010–2020 triggered by mega-construction projects in megacities like the Orange Line Metro Train System and Metrobus System in Lahore, The Pakistan Metrobus System between Islamabad and Rawalpindi, the Multan Metrobus System, the Green Line–Karachi Metrobus, and the Trans Peshawar bus rapid transit system in Peshawar [71,73,88].

5. Conclusions

In this research work, we analyzed bad VR in Pakistan. The variations of VR at various stations are due to natural, meteorological, and anthropogenic triggering factors that affect VR. Overall, bad VR prevails over the region compared to good VR during the study period. Bad VR prevailed in three decades (1980–1990, 2001–2010, and 2011–2020), except during 1991–2000. FSD ranked first, with the highest percentage of bad VR with hourly (92.85%), weekly (92.33%), daily (95.35%), monthly (91.67%), and seasonal (100%) measurements, followed by SKT and MLN. In contrast, KHI, LHR, and ISB have experienced good VR, though they are populated cities and have high pollution levels. Still, their geographical location, climatology, weather pattern, and human activities differ from other stations. MK, MMK, and SS revealed a decreasing trend of VR all over Pakistan. The VR in Punjab and Sindh is decreasing except in the south and southeastern Sindh and KPK. The bad VR in Pakistan implies extremely poor air quality conditions, which are dangerous for the inhabitants of Pakistan. This extensive study based on available data may be considered a baseline VR study. It will be helpful to the ministry of climate change Pakistan and policymakers to mitigate air pollution problems and their triggering factors.