Spatiotemporal distribution and driving factors of atmospheric pollutants in the U-Chang-Shi urban agglomeration, Northwestern China

Air pollutant data

Air pollutant concentrations (PM_2.5, PM₁₀, O₃) were sourced from the CHAP dataset (2000–2022), developed through artificial intelligence-based fusion of MODIS multi-angle aerosol products with ground monitoring and emission inventory data (Wei et al., 2020; Wei et al., 2021a; Wei et al., 2021b; Wei et al., 2022). This 1-km resolution dataset, maintained by the University of Maryland research team and archived at the National Tibetan Plateau Science Data Center (NTPDC), addresses spatial gaps in conventional satellite retrievals.

Meteorological data

The 1-km resolution precipitation and temperature datasets (Peng et al., 2017a; Ding & Peng, 2020; Peng et al., 2019; Peng et al., 2017b) were constructed using Delta spatial downscaling methodology, integrating global climate products from the Climate Research Unit (CRU; 0.5° resolution) and WorldClim (high-resolution). These datasets underwent rigorous validation against observations from 496 meteorological stations across mainland China (including Hong Kong, Macao, and Taiwan) with demonstrated reliability. Spatial coverage excludes South China Sea islands but maintains continental continuity at 0.0083333° resolution (∼1 km).

Monthly potential evapotranspiration (PET) estimates were derived from China-specific 1-km temperature datasets (mean, min, max) archived at NTPDC (Peng et al., 2017a; Ding & Peng, 2020; Peng et al., 2019; Peng et al., 2021). Monthly PET was estimated using the Hargreaves formula. Daily wind speed data were obtained from NOAA’s National Centers for Environmental Information (NCEI). These station-based observations were then interpolated into a continuous national-scale raster dataset at a daily resolution using the Inverse Distance Weighting (IDW) algorithm, based on the geographic coordinates of the stations. Combining with administrative area of prefecture-level city crossing nationwide, after obtaining the statistic value of each station wind power density every day, then getting a day average velocity values of prefecture-level city in one day, month mean velocities got.

Normalized difference vegetation index data

The NDVI was derived from Landsat 5/7/8/9 imagery (2000–2022) via the Google Earth Engine platform. Annual maximum NDVI values were extracted at 30 m spatial resolution after preprocessing and smoothing, with outputs stored in GeoTIFF format (Didan, 2015).

Coefficient of determination

The coefficient of determination (R²) quantifies the proportion of variance in observational data explained by model simulations, and is defined as: (1)${\displaystyle {\mathrm{R}}^2=1-\frac{\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-{\hat{\mathrm{y}}}_{\mathrm{i}}\right)}^2}{\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}^2}}$ where (y_i) denotes the in-situ measurements from observational stations (PM_2.5, PM₁₀, O₃), $\left({\hat{\mathrm{y}}}_{\mathrm{i}}\right)$ represents the CHAP dataset, $\left(\overline{\mathrm{y}}\right)$ is the mean of observational data, and n is the sample size. The numerator, residual sum of squares (RSS), measures the discrepancy between the CHAP dataset and observations, while the denominator, total sum of squares (TSS), reflects the natural variability of observations.

Spearman correlation analysis

Spearman correlation analysis is a non-parametric statistical method for quantifying monotonic relationships between variables, which centers on calculating correlations from variable rank order (rather than raw values) and is applicable to non-linear or non-normally distributed data. The method effectively eliminates the effects of outliers and data distribution patterns by converting the data into rank series, and is able to robustly reveal trend consistency among variables (Spearman, 1904). Given its robustness to non-normal data and outliers, Spearman correlation analysis is widely applied in environmental science to uncover non-linear relationships between diverse geographical factors (Biswas, Chatterjee & Chakraborty, 2020; Zhang et al., 2015).

Within this research, the Spearman rank correlation coefficient (ρ) was employed to measure the magnitude and trend of monotonic relationships between PM_2.5, PM₁₀, O₃ concentrations and meteorological parameters, respectively. The statistical significance was set at a p-value threshold of less than 0.05. The analyses were carried out using SPSS software, and the visualization was completed through Origin. (2)${\displaystyle \mathrm{R}=\frac{\sum _{\mathrm{i}=1}^{\mathrm{N}}\left({\mathrm{X}}_{\mathrm{i}}-{\overline{\mathrm{X}}}_{\mathrm{i}}\right)\left({\mathrm{Y}}_{\mathrm{i}}-{\overline{\mathrm{Y}}}_{\mathrm{i}}\right)}{\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{X}}_{\mathrm{i}}-{\overline{\mathrm{X}}}_{\mathrm{i}}\right)}^2}\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{Y}}_{\mathrm{i}}-{\overline{\mathrm{Y}}}_{\mathrm{i}}\right)}^2}}}$ where: R is the correlation coefficient; Xi is the pollutant concentration in month i, µg/m³; $\overline{{\mathrm{X}}_{\mathrm{i}}}$ is the average of the pollutant concentration in all months; Yi is the value of meteorological elements in month i (when Y is precipitation, temperature, PET, wind speed, the unit is mm, °C, mm, m/s, respectively); $\overline{{\mathrm{Y}}_{\mathrm{i}}}$ is the average value of meteorological elements (the unit is the same as the above); N is the time series length (2000–2022).

Grey relational analysis

Grey relational analysis (GRA) quantifies inter-factor associations by evaluating geometric similarity between data sequences (i.e., the proximity of curve variation trends), where higher relational grades indicate more consistent dynamic changes between paired sequences (Deng, 1982). In this study, grey relational analysis was employed to assess dynamic correlations between pollutant concentrations and key meteorological factors. The analytical procedure involved three sequential steps: first, original data underwent normalization to eliminate dimensional discrepancies. Subsequently, grey relational grades were computed using MATLAB (The MathWorks, Natick, MA, USA) with a distinguishing coefficient set at 0.5. The resultant ranking of relational grades facilitated identification of dominant meteorological factors, which were then cross-validated against correlation analysis results. This methodological implementation strictly adhered to standardized protocols in grey system theory while maintaining compatibility with conventional statistical approaches.

(3)${\displaystyle {\zeta }_{\mathrm{i}}\left(\mathrm{k}\right)=\frac{\underset{\mathrm{i}}{min}\underset{\mathrm{k}}{min}\left|\mathrm{y}\left(\mathrm{k}\right)-{\mathrm{x}}_{\mathrm{i}}\left(\mathrm{k}\right)\right|+\rho \underset{\mathrm{i}}{max}\underset{\mathrm{k}}{max}\left|\mathrm{y}\left(\mathrm{k}\right)-{\mathrm{x}}_{\mathrm{i}}\left(\mathrm{k}\right)\right|}{\left|\mathrm{y}\left(\mathrm{k}\right)-{\mathrm{x}}_{\mathrm{i}}\left(\mathrm{k}\right)\right|+\rho \underset{\mathrm{i}}{max}\underset{\mathrm{k}}{max}\left|\mathrm{y}\left(\mathrm{k}\right)-{\mathrm{x}}_{\mathrm{i}}\left(\mathrm{k}\right)\right|}}$ (4)${\displaystyle {\mathrm{r}}_{\mathrm{i}}=\frac{1}{\mathrm{n}}\sum _{\mathrm{k}=1}^{\mathrm{n}}{\zeta }_{\mathrm{i}}\left(\mathrm{k}\right)}$ where ζi(k) is the grey relational coefficient; r_i is the degree of grey relational grade. y(k) is the normalized parameter value; xi(k) is the normalized comparison value, i is the number of comparison arrays (i = 1, 2...n), and k is the number of indicators for each comparison object (k = 1,2, …, m); ρ is the discrimination coefficient, which usually takes the value of 0.5; based on the results of previous research (Liu, Zhang & Li, 2018), the interval of the correlation value is divided into the degree of strength of the correlation, with [1, 0.8] as strong, (0.8, 0.6) as stronger, (0.6, 0.4) as moderate, (0.4, 0.2) as weaker, and [0.2, 0) as weak.

The Theil–Sen estimator

The slope n of pairs of data points was estimated using the Theil–Sen (TS) estimator (Sen, 1968; Theil, 1992), which is given by Eq. (5). (5)${\displaystyle \mathrm{TS}=\text{Median}\left(\frac{{\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}}{{\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{j}}}\right)}$ where xi and xj are data values at times ti and tj(i > j),  respectively. The slope that was calculated by the TS estimator is a robust estimate of the magnitude of a trend, and the use of the TS slope provides the special benefit of the rejection of inter-annual variability (Neeti & Eastman, 2011). This study employed the MK trend test, which assesses the similarity level between two sets of rankings assigned to the same set of objects. As indicated in references (Kendall, 1938; Hirsch, Slack & Smith, 1982; Valz & McLeod, 1990; Lanzante, 1996), this test relies on the number of inverted object pairs needed to convert one ranking order into the other. The data were first ranked according to the temporal sequence. Subsequently, each data point was sequentially taken as a reference point and compared with all subsequent data points in the time series, as described in reference (Douglas, Vogel & Kroll, 2000). A key challenge in applying time series analysis effectively lies in distinguishing between a genuine change in a remotely sensed variable (such as a vegetation index) and the presence of scene noise (such as smoke or cloud contamination), as noted in reference (Neeti et al., 2012).

The Mann–Kendall trend test

The Mann–Kendall (MK) trend test is designed to determine the existence of trends for each time scale independently. This test involves calculating the MK statistic (S) using the formulas in Eqs. (5) and (6), as cited in references (Dietz & Killeen, 1981; De Beurs & Henebry, 2004; De Jong et al., 2011; Sobrino & Julien, 2011).

(6)${\displaystyle \mathrm{S}=\sum _{\mathrm{k}=1}^{\mathrm{n}-1}\sum _{\mathrm{j}=\mathrm{k}+1}^{\mathrm{n}}\mathrm{sgn}\left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{k}}\right)}$ (7)${\displaystyle \mathrm{sgn}\left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{k}}\right)=\left\{\begin{array}{cc}-1,\hfill & {\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{k}}<0\hfill \\ 0,\hfill & {\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{k}}=0\hfill \\ 1,\hfill & {\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{k}}>0\hfill \end{array}\right.}$

Here, xj and xk denote the data points at time j and k(j > k) respectively, and n represents the total number of data points. However, to statistically measure the significance of the trend, it is essential to calculate the probability associated with the statistic S and the sample size n. When n ≥ 10, the statistic: S approximately follows a normal distribution, with its mean and variance defined by Eq. (7), as stated in references (De Beurs & Henebry, 2004; An, 2007). (8)${\displaystyle \mathrm{Var}\left(\mathrm{s}\right)=\frac{\mathrm{n}\left(\mathrm{n}-1\right)\left(2\mathrm{n}+5\right)}{18}}$ (9)${\displaystyle \mathrm{Z}=\left\{\begin{array}{cc}\frac{\mathrm{S}-1}{\sqrt{\mathrm{Var}\left(\mathrm{S}\right)}},\hfill & \mathrm{S}>0\hfill \\ 0,\hfill & \mathrm{S}=0\hfill \\ \frac{\mathrm{S}+1}{\sqrt{\mathrm{Var}\left(\mathrm{S}\right)}},\hfill & \mathrm{S}<0\hfill \end{array}\right.}$

In this context, S represents the test statistic, Z is the standardized test statistic, and n signifies the length of the time series. At a specified significance level ɛ, when |F| > F1-ɛ/2, the time series data show a distinct change trend. In this study, temporal trends are considered significant when F < −1.96 or F > 1.96. The Z statistic follows the standard normal distribution: a positive value implies an upward trend, while a negative value indicates a downward trend.

HYSPLIT trajectory simulation analysis

HYSPLIT is a commonly used model for analyzing pollutant distribution laws and transportation in the atmosphere, which is constructed by the NOAA of the United States (Draxler & Hess, 1998). The model adopts Lagrangian particles diffusing and Eulerian grids dispersing methods to realize three-dimensional air flow paths and quantify pollutants’ diffusion, dilution, wet dry deposition etc., based on trajectory calculation engine and diffusion prediction module, which can perform analysis from an hour’s scale up to seasons, taking meteorological data at different scales. In this experiment, global forecast meteorology data driven model with lateral space interval 0.25 × 0.25° and vertical layer stratification number 38, describing mid-latitude atmospheric circulation conditions in Eurasia under winter weather, was selected. There are release height AGL such as 500 m green color, 1,000 m blue color, 3,000 m red color in our experiments for two-way propagation processes between seven days’ back tracing and predictive steps further along (Rolph, Stein & Stunder, 2017).

The overall research design and the integrated methodological workflow employed in this study are summarized in Fig. 2. This schematic illustrates the sequential process from multi-source data acquisition and preprocessing to the application of statistical analyses and models, culminating in the final results and conclusions.

Dataset accuracy verification

This study takes the U-Chang-Shi urban agglomeration as the study area, and there are 12 state-controlled stations in this area, namely, seven monitoring stations in Urumqi, seven monitoring stations in Xinjiang Academy of Agricultural Sciences farms, monitoring stations, railway bureau, 31st middle school, 74th middle school, Midong District Environmental Protection Bureau and fee collection office, three monitoring stations in Changji Prefecture, two monitoring stations in Sunshine School and Aiqing Poetry Hall, and two monitoring stations in Shihezi City, and the Agricultural Water Building in Wujiaqu. Three monitoring stations in Changji Prefecture; two monitoring stations in Shihezi City, at the Sunshine School and the Aiqing Poetry Hall; and the Agricultural and Water Building in Wujiaqu, We first validated this dataset by linearly fitting the PM_2.5 dataset for December 2018, January 2019, and February 2019, as well as July, August, and September 2019, to the air quality values provided by the China Air Quality Online Monitoring and Analysis Platform (https://www.cnemc.cn/). Prior to the fitting, we sampled and extracted the values of the Tif images of the dataset at the coordinates of the 12 state-controlled stations in the study area using Arcgis software. From the fitting results of each station (Fig. 3), the R² values range from 0.67 to 0.99.

This R² range can indicate that the independent variable explains the dependent variable to a high degree and the fitting effect is significant. It can be proved that the CHAP dataset of the University of Maryland has a good applicability within the U-Chang-Shi region in the interior of the arid zone.

Monthly and seasonal spatial distribution characteristics of air pollutants

As shown in Figs. 4 and 5, observational data indicate a clear seasonal pattern in O₃ concentrations in the U-Chang-Shi region. The monthly mean concentrations increase significantly from late spring to summer (April–August), from approximately 92.15 µg/m³ in April to a yearly maximum of 144.86 µg/m³ in July and 145.27 µg/m³ in August. This represents an overall increase of nearly 58% from April to midsummer. In contrast, O₃ concentrations in autumn and winter seasons (September–December) remain much lower, with monthly means dropping from 89.47 µg/m³ in September to 57.29 µg/m³ in November and as low as 14.02 µg/m³ in December, representing more than an 84% decline compared to the summer peak.

Spatially, the high-concentration zones (above 130 µg/m³) in summer cover the majority of the region, particularly in industrial areas with intensive NOx and (Volatile Organic Compounds (VOC) emissions. Winter maps show extensive low-concentration zones (below 30 µg/m³) dominated by blue and light-yellow areas. The most polluted period appears in summer, followed by spring, while autumn and winter exhibit the lowest levels.

The spring transitional period (March–May) shows a steady increase (from 84.37 µg/m³ in March to 97.84 µg/m³ in May, a rise of ∼16%). Once autumn temperatures begin to drop (September–November), mean O₃ concentrations fall sharply (>35% within two months).

Overall, the intra-annual fluctuation in O₃ concentration within this region is substantial, ranging from 14 µg/m³ to 145 µg/m³.

These large variations indicate that multiple factors influence air quality: in addition to local industrial production, transportation emissions, and human activities, long-range transport of pollutants also contributes to the observed O₃ levels. This seasonal pattern is consistent with earlier findings (An, 2007), confirming that O₃ pollution reaches its maximum intensity in summer, moderate in spring, and minimal in autumn–winter periods. The spring transitional rise is probably due to gradually rising temperatures and increasing VOC emissions. The rapid autumn decline is driven by the reduction of photochemical reaction rates.

The monitoring results (Figs. 6 and 7) show that the winter season (December–February) is the most polluted period in the U-Chang-Shi region. In January, the mean concentration at Changji reached 171 µg/m³, representing the highest recorded value. High-value contamination zones (>150 µg/m³) appear in western Shihezi, central Urumqi, northwestern Changji, the northern part of Shawan, and the southern part of Wujiaqu during winter.

For particulate matter, spring is the cleanest season: may recorded a mean concentration of 31.52 µg/m³, which is about 40% lower than March (52.67 µg/m³). Summer concentrations are mostly under 30 µg/m³—over 80% lower than the January peak. Autumn values gradually rebound from 30.15 µg/m³ in September to 46.38 µg/m³ in November, an increase of nearly 50%. Across the year, PM_2.5 fluctuates between 14 and 171 µg/m³.

Such seasonal variation is directly linked to meteorological conditions, emission intensities, and the effect of topography in trapping pollutants during the stagnant winter period. The winter peaks are influenced by heating activities and unfavorable atmospheric dispersion, while the low concentrations in summer are aided by stronger mixing and precipitation scavenging.

Figures 8 and 9 show that PM₁₀ concentrations also peak during winter, with the highest value (225 µg/m³) measured in Urumqi in January, and the lowest value (43 µg/m³) in July. Spring mean concentrations decreased by 13.53% from March to May, while autumn increased by 17.75% from September to November. The seasonal sequence is winter > spring > autumn > summer. During spring, an unusual high-value PM₁₀ zone is visible at the border between southeastern Urumqi and Turpan, while most of the rest of the region has moderate levels.

The higher PM₁₀ levels in spring and autumn compared to PM_2.5 are due to the significant contribution of coarse particles, including dust events and soil particles. The distinctive southeastern Urumqi–Turpan high-value zone is likely linked to local topography, prevailing wind directions, and dust source proximity.

Annual concentration trends and statistical characteristics of air pollutants

Figure 10 presents a comprehensive analysis of the spatiotemporal evolution of PM_2.5 and PM₁₀ concentrations in the “U-Chang-Shi” region from 2000 to 2022. Figure 10A illustrates the interannual variability of PM_2.5 concentrations in the “U-Chang-Shi” region from 2000 to 2022. Over this 22-year period, annual mean PM_2.5 levels exhibited a fluctuating decline, decreasing from 50.18 µg/m³ in 2000 to 40.26 µg/m³ in 2022 (17.92% reduction). Concurrently, PM₁₀ concentrations followed a similar trajectory, declining from 102.33 µg/m³ to 77.67 µg/m³ (24.09% reduction; see Fig. 10B). Notably, both pollutants displayed synchronized temporal patterns, with peaks observed in 2012 followed by substantial reductions post-2015.

City-specific analyses revealed heterogeneous trends from 2000 to 2022. PM_2.5 and PM₁₀ concentrations decreased significantly in Urumqi (by 41.21% and 43.72%, respectively), Shawan (30.83% and 34.10%), and Changji (22.83% and 26.83%). In contrast, the declines in Shihezi were much smaller, with a PM₁₀ reduction of only 4.49%. Notably, Wujiaqu was an exception, experiencing a 10.09% increase in PM_2.5 with its PM₁₀ level remaining largely stable (a slight decrease of 1.57%).

The temporal dynamics of PM_2.5 from 2000 to 2022 exhibited three distinct phases (Fig. 10C): a period of irregular increase until 2012, a sharp decline from 2013 to 2015, and a steady decrease following a brief rebound in 2016. Wujiaqu City constituted an exception to this regional pattern, showing a continuous upward trend. Similarly, PM₁₀ concentrations followed a comparable trajectory, characterized by a primary peak between 2012 and 2014, a decrease around 2017, and a sustained decline thereafter (Fig. 10D).

From the perspective of Urumqi’s regional air pollution and environmental management, 2012 was the peak of regional pollutant concentration average, and the beginning of air pollution and environmental management focus. “Coal to gas” project began to implement, Urumqi City in six months to complete the transformation of 12,900 tons of steam tons, the particulate matter (PM_2.5, PM₁₀) can also be seen in the graph of a significant decline. In the process of pollution management, 2016, 2017 particulate matter (PM_2.5, PM₁₀) appeared in the concentration value of the phenomenon of a small rebound after the timely containment, and strong, effective management of pollutants, and after 2020 leveled off. This shows that the remediation of the atmospheric environment has achieved gradual results.

Figure 11 illustrates the contrasting temporal dynamics of O₃ concentrations across the study region. As summarized in the heatmap (Fig. 11A), O₃ levels remained relatively stable from 2000 to 2016, with annual variations ranging from −0.14% to +0.6%. During this period, Shihezi and Wujiaqu consistently exhibited lower baseline concentrations compared to other urban centers. A notable regime shift occurred after 2016, as clearly depicted in the temporal trend analysis (Fig. 11B), marked by accelerated growth rates of 4.2–4.6% annually during 2017–2022. The trend graph particularly highlights anomalous spikes in 2017, with Urumqi and Shihezi experiencing increases of 7.1% and 16.7% respectively—a pattern potentially linked to extreme thermal episodes or intensified industrial activity. Although a transient decline emerged in 2020, the post-2016 upward trajectory remained unaltered, with all cities exceeding 7% annual growth during peak phases. The persistent growth momentum observed after 2016, especially the rapid increase starting in 2017 shown in Fig. 11B, appears largely unabated. This trend may be attributed to the surge in industrial emissions in recent years. For instance, the commissioning of the chemical park in Shihezi City in 2017, as reflected in the notable 16.7% concentration increase for that year, provides a plausible explanation for the observed spatial–temporal pattern.

Spatial distribution patterns of multi-year average pollutant concentrations

Figure 12A reveals pronounced spatial heterogeneity in multi-year average PM_2.5 concentrations (2000–2022) across the “U-Chang-Shi” region, with values ranging from a minimum of 30.70 µg/m³ in the southern mountainous areas to a maximum of 75.90 µg/m³ in Urumqi’s metropolitan core. High-value zones (>70 µg/m³) form a contiguous east–west corridor stretching from Urumqi through Changji Jundong to Shihezi, corresponding to regions of dense population, intensive industrial activity, and heavy traffic emissions. Areas with medium concentrations (50–65 µg/m³) are primarily located in northwestern Changji, Shawan, and Wujiaqu, while low-value zones (<40 µg/m³) are distributed mainly in peripheral natural landscapes, particularly the high-altitude southern mountains where cleaner air is maintained due to strong atmospheric dispersion and limited anthropogenic sources.

This stark north–south gradient reflects the combined influence of topography—with mountain ranges acting as barriers to pollutant transport—and land-use patterns, where urban-industrial belts generate persistent PM_2.5 hotspots, contrasting with the low-emission natural terrain in the south.

Figure 12B illustrates the multi-year (2000–2022) average spatial distribution of PM₁₀ in the “U-Chang-Shi” region, with concentrations ranging from a minimum of 71.80 µg/m³ in the southern mountainous areas to a maximum of 142.82 µg/m³ in the core urban zone of Urumqi. High-concentration zones (>130 µg/m³) form a continuous pollution belt encompassing Urumqi, Shihezi, and Wujiaqu, and extending to the Urumqi–Turpan border. Surrounding areas show a gradual decrease in concentration toward the south of the Wuquan Canal, with markedly lower values (<85 µg/m³) south of the Shawan Mountain range.

A particularly notable anomaly is evident in the southeastern Urumqi mountain area, where PM₁₀ levels remain unexpectedly high (∼125–135 µg/m³), a feature absent in other particle size distribution maps. This anomaly is possibly related to complex topography—alternating high valleys and low hollows with vegetation (e.g., poplar forests) and sulfur-rich geological zones. Under such terrain conditions, airflows shift abruptly, especially near nightfall, when cold air descends the slopes and accumulates in valley floors, forming a stable stratification that traps particulates.

This persistent spatial pattern—urban-industrial hotspots contrasting against cleaner high-altitude southern regions—highlights the combined influence of anthropogenic emissions and terrain-driven atmospheric processes in shaping PM₁₀ distributions.

Figure 13 shows that, compared to PM_2.5 and PM₁₀, elevated O₃ concentrations occupy a more extensive and clearly defined spatial range across the “U-Chang-Shi” region, with values ranging from 62.70 µg/m³ to a maximum of 96.50 µg/m³. The highest O₃ levels are concentrated in the southeastern part of Urumqi—particularly in areas adjacent to Turpan City—and the southern part of Changji City, where average concentrations exceed 90 µg/m³.

Medium-value regions (approximately 75–85 µg/m³) are typically adjacent to these hotspots, extending along the southwestern mountainous zones near the northern foothills of the Tianshan Mountains. These areas present sharp concentration boundaries following the topographic transitions in the mountainous terrain.

In contrast, low-value areas (<70 µg/m³) are primarily located in the central “U-Chang-Shi” region, including the main urban core of Urumqi, where complex interactions between atmospheric dynamics and precursor emissions limit O₃ accumulation compared to surrounding elevated zones.

The spatial gradient—high values in the topographically elevated southern margins and low values in central depressions—reflects the influence of terrain-driven meteorological processes (e.g., enhanced photochemical activity in high-altitude, sun-exposed slopes) combined with local precursor distribution patterns.

MK trend of air pollutants

Overall, the TS Slope values (Fig. 14B) exhibit a high consistency with the MK-Z values of PM_2.5 (Fig. 14A) in terms of spatial distribution. The central and northern regions show a larger slope, and the MK-Z is significantly positive (|Z| ≥ 1.96, corresponding to p < 0.05), reflecting a significant increasing trend in PM_2.5 concentration, with the maximum annual increase reaching up to 4.04 µg m⁻³ yr⁻¹. In the southern and southwestern mountainous areas, the slope is negative and the MK-Z is significantly negative (p < 0.05), indicating a significant decreasing trend, with the maximum annual decline reaching −4.80 µg m⁻³ yr⁻¹. Some marginal and central areas show Z values close to 0 (—Z—< 1.96, corresponding to p ≥ 0.05), indicating that the change trends are statistically insignificant. Notably, the MK-Z values and TS Slope values in the Uga District and Shihezi City are both at a high level, and the significance test (p < 0.05) reflects that PM_2.5 concentrations in these areas have significantly increased in recent years.

The TS Slope values (Fig. 14D) and MK-Z values of PM₁₀ (Fig. 14C) exhibit broadly consistent spatial patterns. In the central–northern part of the region, both TS Slope values are positive and MK Z values are significantly positive (|Z| ≥ 1.96, p < 0.05), indicating a statistically significant increasing trend, with annual growth rates reaching up to 3.11 µ g m⁻³ yr⁻¹. In the southern and southwestern mountainous areas, TS Slope values are negative and MK Z values are significantly negative (p < 0.05), reflecting a significant decreasing trend, with the largest decline reaching −4.25 µg m⁻³ yr⁻¹. In some marginal and central areas, MK Z values are close to zero (|Z| < 1.96, p ≥ 0.05) and TS Slope values are near zero, suggesting statistically insignificant long term changes. Notably, Wujiaqu City and Shihezi City have both high TS Slope values and significantly positive MK Z values (p < 0.05), highlighting a pronounced and statistically significant upward trend in these areas.

Figure 15 presents the spatial heterogeneity in both the statistical significance and magnitude of O₃ concentration trends across the study region from 2000 to 2022. As shown in Fig. 15A, the per-pixel MK-Z values analysis reveals a distinct north–south divergence in trend significance. Significantly increasing O₃ trends (Z ≥ 1.96, p < 0.05) are observed in the northeastern parts of Changji and Urumqi, as well as in central Shawan, while significantly decreasing trends (Z ≤ −1.96, p < 0.05) dominate the southern and western regions. The central and transition zones show no significant trend (—Z—<1.96, p ≥ 0.05). Figure 15B complements this by quantifying the rate of change via the TS Slope, with high values (up to 3.35 µg m⁻³ yr⁻¹) in the central and northeastern areas and low to negative values (down to −0.30 µg m⁻³ yr⁻¹) in the south and southwest. Regions with significant Z-values but modest slopes reflect gradual yet persistent O₃ increases, whereas areas combining high significance and large slopes indicate rapid, substantial changes. This spatial pattern underscores the influence of regional-specific drivers on O₃ pollution dynamics.

In regions where MK-Z values indicate significant increases but TS Slope values remain relatively low, the trend represents a gradual yet significant increase. This may be attributed to steady growth in emission source intensity or inhibition of ozone production by topographic factors (e.g., mountainous terrain). Conversely, areas with significant MK-Z decreases but small TS Slope declines exhibit gradual yet significant decreases, potentially due to stable ecological purification capacity or limited effectiveness of mitigation measures. In the central urban agglomeration and adjacent areas, both MK-Z and TS Slope values approach 0, defining a zone of insignificant trends and slow changes. These regions are influenced by stable anthropogenic emissions and unaltered natural factors, leading to a relatively balanced ozone production-depletion dynamic.

Correlation analysis between pollutants and meteorological factors with GRD

Spearman correlation analysis (Fig. 16A) revealed strong and significant negative correlations between PM_2.5 and NDVI (ρ = −0.89), T, (ρ = −0.87), PET (ρ = −0.83), and WIND (ρ = −0.77), with P (ρ = −0.71) also moderately negatively correlated. These results suggest that stagnant meteorological conditions—characterized by low PET, low T, weak Winds, and sparse vegetation—are conducive to PM_2.5 retention, while higher precipitation enhances wet scavenging. GRA results (Fig. 16B) indicated a different emphasis, ranking P (0.6229) highest, followed by NDVI (0.56846), PET (0.55139), T (0.52758), and WIND (0.52737), underscoring precipitation’s dominant role in wet deposition and vegetation’s importance in pollutant capture. To assess robustness, we repeated Spearman analyses under |ρ| ≥ 0.3 and |ρ|≥ 0.7 thresholds and compared with GRA rankings. The order of dominant drivers (precipitation–vegetation–PET) remained consistent, confirming that our identification of key meteorological factors is stable across analytical approaches and parameter settings.

Spearman correlation analysis (Fig. 16C) showed that PM₁₀ had the strongest negative correlations with NDVI (ρ = −0.78) and T (ρ = −0.73), indicating that higher vegetation coverage and warmer conditions can reduce PM₁₀ levels through interception, adsorption, and enhanced convection. PET (ρ = −0.67) and WIND (ρ = −0.56) also showed negative relationships, suggesting that increased atmospheric motion supports pollutant dispersion, while P (ρ = −0.56) reflected moderate wet scavenging effects. GRA results (Fig. 16D) ranked P (0.63685) highest, followed by NDVI (0.6056), PET (0.60073), WIND (0.59184), and T (0.57164), again emphasizing precipitation’s primary role in PM₁₀ removal. Robustness checks using Spearman thresholds of —ρ— ≥ 0.3 and —ρ— ≥ 0.7 produced rankings consistent with GRA, with P, vegetation, and PET always emerging as the top three drivers. This alignment demonstrates that the meteorological influences on PM₁₀ are stable across different statistical criteria.

Spearman correlation and grey relational analysis (Figs. 17A–17B) together revealed the main meteorological drivers of O₃ variation. PET emerged as the most influential factor in both methods (ρ = 0.90; GRD = 0.727), indicating a dual role as a monotonic predictor and a synchronous trend driver, likely reflecting its control over surface energy balance and photochemical reaction rates. WIND ranked second in importance, with a high Spearman coefficient (ρ = 0.83) and substantial GRD (0.681), suggesting a near-linear role in pollutant transport and dispersion. T showed the largest discrepancy between methods—its high Spearman correlation (ρ = 0.88) contrasted with the lowest GRD (0.553)—implying that thermally driven photochemical processes may be subject to non-stationary behavior or temporal variability. P (ρ = 0.66; GRD = 0.695) and NDVI (ρ = 0.76; GRD = 0.615) both exhibited moderate linear correlations but relatively higher grey relational degrees, suggesting that their influence on O₃ likely operates through lagged or non-linear pathways. A robustness check—repeating Spearman analysis under correlation thresholds —ρ— ≥ 0.3 and —ρ— ≥ 0.7 and comparing with GRA rankings—confirmed PET and WIND as consistently dominant drivers, with P, NDVI, and T contributing secondary effects. These results highlight the complementary value of combining linear (Spearman) and non-linear (GRA) frameworks to capture both immediate correlations and underlying trend synchrony in detecting O₃ driving mechanisms.

Spearman correlation analysis quantifies the statistical correlation between meteorological elements and pollutants, reflecting the strength of linear correlation of variables and the immediate impact of meteorological elements on pollutants. Grey correlation analysis is based on the trend similarity between meteorological elements and pollutant sequences, analyzing the potential mechanism of elemental dynamics on pollutants. Due to the differences in the logic and methodology of statistical correlation and trend similarity, the results of the impact assessment of meteorological factors on pollutants are different in terms of the ranking of the degree of correlation and the characterization of their effects. Furthermore, related studies have also mentioned the differences between the methods (Zheng et al., 2023).

PM_2.5 and PM₁₀ (Fig. 18) exhibited a strong positive correlation (ρ = 0.89), reflecting their overlapping anthropogenic emission sources. Coal-fired power plants and motor vehicle exhausts release both fine (PM_2.5) and coarse (PM₁₀) particles, and incomplete fossil fuel combustion often produces a bimodal particle size distribution. Meteorological influences on the two pollutants are largely consistent: under stable atmospheric conditions, low WIND, and temperature inversion layers, the diffusion of fine particles and deposition of coarse particles are inhibited, leading to their synchronous accumulation. During rainy periods, both are subject to wet deposition through humidity-driven scavenging.

A strong negative correlation was observed between PM_2.5 and O₃ (ρ = −0.84), indicative of a photochemical competition for common precursors, such as nitrogen oxides (NO_x). PM_2.5 also exerts a radiative shielding effect by absorbing and scattering ultraviolet radiation, thereby weakening key photolysis reactions essential for O₃ formation.

For PM₁₀, mineral aerosols (e.g., dust) can catalyze O₃ decomposition (O₃ → O₂). This effect is especially evident during northern spring dust events, where increases in PM₁₀ concentrations correspond to decreased O₃ levels. As coarse particles (2.5–10 µm) are less directly involved in photochemical reactions, the negative PM₁₀–O₃ correlation likely reflects indirect associations, such as dust events accompanied by strong WIND, rather than direct atmospheric chemistry.

HYSPLIT forward and backward trajectory analysis

The HYSPLIT_4 developed by the NOAA Air Resources Laboratory was employed to simulate air mass backward and forward trajectories, in order to investigate potential transboundary pollution transport affecting the study area. The model was driven by the Global Data Assimilation System (GDAS) meteorological dataset with a horizontal resolution of 0.25° × 0.25°, a vertical layering of 38 levels, and a temporal resolution of 3 h.

The trajectory simulations were initialized from the central Urumqi site (43.97°N, 87.34°E) at three release heights above ground level (AGL): 500 m (boundary layer), 1,000 m (lower free troposphere), and 3,000 m (high-level transport layer), to account for vertical variability in transport pathways. Both 72-hour backward and forward trajectories were computed at 1-hour time steps to capture synoptic-scale transport events. The emission assumption was set as an instantaneous and homogeneously mixed release of a non-reactive tracer, with no consideration of dry/wet deposition or chemical transformation, as the primary aim was to track air parcel origins and dispersion pathways rather than perform chemical fate modeling.

Model uncertainties may arise from several factors, including (i) errors in meteorological input data, especially in complex terrain; (ii) sensitivity of the calculated trajectories to the choice of starting location, time, and release height; (iii) cumulative divergence of trajectories at extended simulation lengths beyond 72 h; and (iv) model assumptions that exclude atmospheric chemistry, deposition processes, and small-scale turbulence. To assess sensitivity, additional simulations were conducted with varied release heights (±500 m) and durations (48 h vs. 72 h). Results indicate that higher-altitude trajectories are generally more stable, whereas near-surface trajectories in winter are more variable due to valley winds and thermal inversions. These uncertainty sources were considered when interpreting the trajectory results in the discussion section.

As shown in Fig. 19A, regarding back-tracing trajectories, we find that low-level air parcels released at 500 m AGL originated from the northwestern region, starting from eastern Kazakhstan, then approaching the source area over the next few days, about 800 km. At the same time, this air mass descends gradually from 1,500 m to 500 m along the way. This can be attributed mainly to the airflow descending through the northern slope of Tianshan Province due to terrain. As the high-level trajectory moves eastward towards Central Asia, specifically from the northern part of the Caspian Sea within the Volga River Basin, there is almost no movement at high altitudes, with a very small upward speed reaching up to only 0.2 cm/s, indicating stable conditions for upper-level westerlies.

Further analysis shows that the large-scale circulation is affected by strong west winds, which contribute significantly to the high-level advection transfer caused by the Siberian High Pressure system (a high pressure center greater than 1,040 hPa), forcing low-level northerly cold air southwards. On one hand, mesoscale topography has an important impact because of the obstruction of mountains around Tien Shan Mountain Range or forcing the mountain uplift phenomenon. On the other hand, lake Balkhash’s heat contrast also affects the direction of low-level flow. Additionally, boundary layer diffusion such as nocturnal inversion processes inhibits the vertical diffusion of low-level air currents, causing them to stagnate near the transport layer.

As shown in Fig. 19B, the low-level trajectory shows southeast movement and forms a zigzag line about 600 km, mainly affected by the friction of the northern foothills of Tianshan Mountains and local circulation; the middle-layer acceleration to the east and reaching the western part of Mongolian Plateau is due to high-level momentum downward transmission; intermediate-level accelerated eastward to the western region of the Mongolian Plateau, indicating that high level also contributed to the momentum transportation. In general, the high-level trajectory is mainly controlled by strong westerlies, and its track is relatively straight, transporting the longest distance and finally arriving at North China. Trajectory scattering is caused by the influence of westerly jet stream, the effect of pressure gradient on Mongolia, and the influence of the topography of Tianshan Mountain; the curved low-level route may be related to the temperature, barometer in Dzungarian Basin, and mountain leeward wave dynamic conditions.

These simulations can provide us with more information to better understand transboundary pollutants transported from Central Asia during winter, especially for the pollution transport in central Eurasia region. The result has demonstrated the prominent westerly transport channel at an altitude of 3,000 m, which could provide some basic evidence to study long-range high-level sand-dust aerosol transport in further works.