Complex network analysis of fine particulate matter (PM_2.5): transport and clustering

In this paper, complex network theory has been applied to reveal the transport patterns and cooperative regions of fine particulate matter (PM_2.5) over China from 2015 to 2019. The results show that the degrees, weighted degrees, and edge lengths of PM_2.5 cities follow power law distributions. We find that the cities in the Beijing–Tianjin–Hebei–Henan–Shandong (BTHHS) region have a strong ability to export PM_2.5 pollution to other cities. By analyzing the transport routes, we show that a mass of links extends southward from the BTHHS to the Yangtze River Delta (YRD) regions with 1 or 2 d time lags. Hence, we conclude that earlier emission reduction in the BTHHS region and early warning measures in the YRD region will provide better air pollution mitigation in both regions. Moreover, significant links are concentrated in wintertime, suggesting the impact of the winter monsoon. In addition, all cities have been divided into nine clusters according to their spatial correlations. We suggest that the cities in the same clusters should be regarded as a whole to control the level of air pollution. This approach is able to characterize the transport and cluster for other air pollutants, such as ozone and NO_x.

1 Introduction

The Earth system behaves as a complex self-regulating system comprised of atmosphere, hydrosphere, cryosphere, lithosphere, and biosphere, with highly nonlinear interactions and feedbacks between the component parts (Steffen et al., 2015). With the increasing understanding of interactions between physical, chemical, biological, and human processes, a new “science of the Earth”, known as Earth system science (ESS), has been initiated (Steffen et al., 2020). Facilitated by its various tools and approaches, ESS has introduced some new concepts and theories, the most important of which is the concept of the Anthropocene (Malm and Hornborg, 2014). In the Anthropocene era, haze events have occurred frequently in China, and the problem of air pollution has received wide attention from the government, scholars, and the public in China (Huang et al., 2014; Sheehan et al., 2014).

Fine particulate matter (PM_2.5) is the primary cause of haze pollution (Ding et al., 2016; Cai et al., 2017). It has adverse influences on human health, atmospheric visibility, and global climate change (Liang et al., 2016; Liaoet al., 2017). PM_2.5 pollution is generated from both anthropogenic and natural sources, including primary aerosols and secondary aerosols that are produced in the atmosphere through the chemistry of precursor gases (Squizzato et al., 2012). In recent years, it has also been recognized that air pollution in a given area is influenced not only by the air pollutant emissions within the same region but also by the transport of air pollutants from other regions. Based on trajectory clustering methods, Li et al. (2019) concluded that regional PM_2.5 transmission has become the key factor driving severe haze in Beijing. By using the positive matrix factorization approach, Khuzestani et al. (2017) revealed that remote transmission accounted for approximately 77 % of the PM_2.5 concentration in the Ordos region. Furthermore, PM_2.5 transmissions are also examined using model simulations. For example, Wang et al. (2014) quantified the regional contribution of PM_2.5 in southern Hebei by using Mesoscale Modeling System Generation 5(MM5) and the Models-3 Community Multiscale Air Quality (CMAQ) modeling system. Zhang et al. (2017) investigated the effect of regional pollution transport based on the GEOS-Chem chemical transport model and its adjoint. These studies suggest that curbing air pollution has not been a local issue, and the regional coordination could be an effective approach to improve the air quality of the regional atmospheric environment. In 2012, The 12th Five-Year Plan on Air Pollution Prevention and Control in Key Regions approved proposed dividing China into three key regions to jointly prevent air pollution, i.e., the Beijing–Tianjin–Hebei (BTH) region, Yangtze River Delta (YRD), and the Pearl River Delta (PRD), and major urban agglomerations such as Lanzhou–Xining, Wuhan and surrounding areas, and Shaanxi and Guanzhong city (MEP, 2012). However, this kind of region division ignores the nonlinear transport characteristics of PM_2.5 concentrations; furthermore, considerable discrepancies exist in the above studies of PM_2.5 transmission in different cities and regions during different air pollution periods. For example, the transport from the BTH region to the YRD is significant during the hazing periods (Huang et al., 2020). High PM_2.5 to the southwest and south of Beijing is related to the PM_2.5 transmission in Baoding and Hengshui in Hebei Province and Dezhou, Liaocheng, Heze, Jining, and Zaozhuang in Shandong Province, respectively (Li et al., 2015). Hence, PM_2.5 transport for the whole of China over a long-time period is not been fully understood; furthermore, the traditional approaches adopted in the above studies did not fully consider the nonlinear transport processes between cities.

Methods are required that help to unveil the transport processes at the national scale. In addition, it is important to quantify their spatial and temporary interactions between cities. During the last 2 decades, complex network theory has been applied to reveal the statistical and dynamic topological features in complex systems (Fountalis et al., 2014; Feldhoff et al., 2015). In complex networks, geographical locations are considered to be nodes. Links represent communications between time series of nodes, and their strength is measured by the cross-correlation between records (Castrejon-Pita and Read, 2010). The network-theory-based approach has been used to study teleconnection patterns (Zhou et al., 2015; Boers et al., 2019; Ying et al., 2019), El Niño events (Yamasaki et al., 2008; Ludescher et al., 2013, 2014), the North Atlantic Oscillation (Guez et al., 2012), Atlantic Multidecadal Oscillation (Wyatt et al., 2012), and Rossby waves (Wang et al., 2013; Ying et al., 2020). This approach is also useful in the studies of atmosphere environment systems, especially enabling us to investigate the nonlinear spatiotemporal dynamics between air pollution agents. Such nonlinear relationships are critical for assessing the intrinsic dynamics of atmospheric pollution systems, but traditional statistical or model simulation methods are difficult to reveal. The network-theory-based approach has been used to uncover the correlation pattern of PM_2.5 concentrations (Zhang et al., 2018), analyze the PM_2.5 spillover routes in BTH region cities (Li et al., 2019), discriminate between urban and rural tropospheric ozone (Rafael et al., 2019), and quantify the interaction between upper-air conditions and surface PM_2.5 concentrations (Zhang et al., 2019). It is obvious that complex network methods are valuable tools for depicting and quantifying air pollution transmission and cluster among cities. In addition, for traditional model simulation, numerous parameters are needed in the simulation process. In contrast, complex network theory is performed based on time series of field observations, so the estimation process is faster and more economic.

In the present study, we attempt to explore the transport and cluster of PM_2.5 based on complex networks, and in the next section we introduce the data and methods. The patterns of PM_2.5 concentrations and their transport features and demarcation regions are presented in Sect. 3. Finally, the summary and discussion are detailed in Sect. 4.

2 Data and methods

2.1 Data

The daily PM_2.5 concentration data for 336 cities over China from 1 January 2015 to 31 December 2019 are used in this study. These raw data were acquired from the China National Environmental Monitoring Centre (CNEMC). We then pre-processed these data according to the needs of the ambient air quality standard for the validity of air pollutant concentration data. Specifically, the missing values in the PM_2.5 data are excluded. Following this, the error values, e.g., negative values and those larger than 900 mg m⁻³ on a given day for a given year, are removed, and for these years we deleted the data corresponding to those days. As a result, we obtained data for 360 valid days per year (data on 9 January, 1 April, 6 July, 5 September, and 29 November are removed), and thus the total length is 5×360 (1800 d).

The anomaly records of PM_2.5 are adopted, where the anomalies are obtained by subtracting the daily averages and dividing them by the corresponding standard deviations to remove the seasonal cycle.

2.2 Methods

The network construction includes three steps. First, we calculate the weight of the edges between nodes. Second, we apply a shuffled procedure to identify a certain threshold. Third, we calculate network typological metrics to determine the interaction strength between two nodes. Below, we detail each step.

2.2.1 Step 1: the calculation of the weight links between nodes

The anomalous PM_2.5 time series of each node i is represented as δS_i(t), where i is the node index. Similar to earlier studies (Gozolchiani et al., 2011; Ying et al., 2020), we define X_i,j(τ) as the time-delayed cross-correlation function for PM_2.5 nodes (i and j), i.e., δS_i(t) and δS_j(t). For τ>0,

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& X_{i,j}\left(\mathit{\tau }\right)=\\ & {\displaystyle \frac{〈\mathit{\delta }{S}_i(t-\mathit{\tau })\mathit{\delta }{S}_j\left(t\right)〉-〈\mathit{\delta }{S}_i(t-\mathit{\tau })〉〈\mathit{\delta }{S}_j\left(t\right)〉}{\sqrt{〈{\left(\mathit{\delta }{S}_j(t-\mathit{\tau })-〈\mathit{\delta }{S}_j(t-\mathit{\tau })〉\right)}^{\mathrm{2}}〉}\cdot \sqrt{〈{\left(\mathit{\delta }{S}_j\left(t\right)-〈\mathit{\delta }{S}_j\left(t\right)〉\right)}^{\mathrm{2}}〉}}}\end{array},\end{array}$$

where τ denotes the time lag, which is in a range between 0 and +30 d, and $X_{i,j}\left(\mathit{\tau }\right)=X_{j,i}(-\mathit{\tau })$. The bracket is the average over the time period of our study. We quantify the strength of the correlations as follows (Gozolchiani et al., 2011; Guez et al., 2014):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{W}_{i,j}^{\mathrm{pos}}={\displaystyle \frac{\max \left(X_{i,j}\right)-\mathrm{mean}\left(X_{i,j}\right)}{\mathrm{SD}\left(X_{i,j}\right)}},\text{(3)}& {\displaystyle }& {\displaystyle }{W}_{i,j}^{\mathrm{neg}}={\displaystyle \frac{\min \left(X_{i,j}\right)-\mathrm{mean}\left(X_{i,j}\right)}{\mathrm{SD}\left(X_{i,j}\right)}}.\end{array}$$

In this approach, max( ), mean( ), min( ), and SD( ) denote the maximum, minimum, mean, and standard deviation of the cross-correlation function X_i,j(τ), respectively. The deviations in the link identification caused by persistence or autocorrelation in the records are reduced through dividing the SD(X_i,j). We defined the maximum and minimum of X_i,j as $P_{i,j}^{\mathrm{pos}}$and $P_{i,j}^{\mathrm{neg}}$, respectively; ${\mathit{\tau }}_{i,j}^{\mathrm{pos}}$ and ${\mathit{\tau }}_{i,j}^{\mathrm{neg}}$ represent the maximum and minimum values of X_i,j(τ), respectively; and the sign of ${\mathit{\tau }}_{i,j}^{\mathrm{pos}}$ (or ${\mathit{\tau }}_{i,j}^{\mathrm{neg}}$) represents the direction of each positive (or negative) link. When ${\mathit{\tau }}_{i,j}^{\mathrm{pos}}>\mathrm{0}$, the link is regarded as coming from node i and pointing to node j. When ${\mathit{\tau }}_{i,j}^{\mathrm{pos}}<\mathrm{0}$, the link is regarded as pointing away from node j toward node i. The direction is undefined when ${\mathit{\tau }}_{i,j}^{\mathrm{pos}}=\mathrm{0}$. The definitions are similar for the negatively weighted links.

The adjacency matrix is defined as follows:

$$\begin{array}{}\text{(4)}& {\mathbf{A}}_{i,j}^{\mathrm{pos}}=\left(\mathrm{1}-{\mathit{\delta }}_{i,j}\right)H\left(W_{i,j}^{\mathrm{pos}}-Q\right)\end{array}$$

where δ_i,j is the Kronecker delta introduced to avoid self-loops in the network and H(x) is the Heaviside step function ($H(x>\mathrm{0})=\mathrm{1}$ and $H(x<\mathrm{0})=\mathrm{0}$). Q denotes a certain threshold value. The definitions are similar for the negatively weighted links. We constructed networks by pruning the links for which the statistical significance was below a certain threshold (Guez et al., 2014). The threshold is determined according to the shuffle method, which is explained in detail in the next section.

2.2.2 Step 2: the identification of the critical threshold

In the shuffled case, the order of years is permutated and the order of days within each year is maintained for each pair of nodes (Ying et al., 2020). For each link, we selected one of two nodes randomly and then shuffled this time series by persisting the order of days in each year and changing the permutation of years several times. We then calculated the cross-correlation function and weight links for the shuffled datasets. The shuffling procedure represents the properties of statistical quantities and the autocorrelations of the original records, which may introduce unrealistic links. We only considered the link weights in the original network that are significantly higher than values in the shuffled case as a real link; otherwise, they are classed as spurious links. According to the principles mentioned above, Fig. 1 depicts the research process and integration of analytical tools.

Figure 1 The flow chart of the method with complex network analysis. Please note that the above figure contains disputed territories.

2.2.3 Step 3: the determination of network typological metrics

The degree is the most common application for measuring complex networks. A link that points toward a node is referred to as an “in-degree” link, and a link that points away from a node is considered as an “out-degree” link. The “in-weights” and “out-weights” for degrees of node i are denoted as In(w)_i and Out(w)_i, representing the total incoming (or outgoing) weighted links, respectively, as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{In}(w{)}_i=\sum _j{A}_{j,i}{W}_{j,i},\text{(6)}& {\displaystyle }& {\displaystyle }\mathrm{Out}(w{)}_i=\sum _j{A}_{i,j}{W}_{i,j}.\end{array}$$

The in- and out-weighted degrees represent a node's dependence on its surrounding nodes and the influence of the node on the surroundings nodes, respectively. Nodes with higher values in the network indicate a larger amount of connection with other nodes, whereas zero values indicate that the node is isolated.

Figure 2 Positive link weights as a function of geographical distances D_i,j for (a) $W_{i,j}^{\mathrm{pos}}$ and (b) $P_{i,j}^{\mathrm{pos}}$ for real (blue) and shuffled (red) data. Panels (c, d) are the same as panels (a, b) but for negative links.

The Girvan–Newman algorithm is used to explore regional division in the networks. In binary networks, the quality of community structure is typically measured by the modularity (Q) function (Newman, 2006). A high value of Q suggests a strong division of a network into clusters. Nodes in the same community may have the same properties. The Q in networks is defined as follows:

$$\begin{array}{}\text{(7)}& Q={\displaystyle \frac{\mathrm{1}}{\mathrm{2}M}}\sum _{i,j}\left[\left({\mathbf{A}}_{i,j}-{\displaystyle \frac{k_i{k}_j}{\mathrm{2}M}}\right)\mathit{\delta }\left({\mathit{\sigma }}_i,{\mathit{\sigma }}_j\right)\right],\end{array}$$

where k_i and k_j are the weight of node i and j, respectively, A_i,j is the adjacency matrix, δ is the membership function, and M is the number of edges.

Figure 3 PDF of positive (a) and negative (b) link weights for the original data and shuffled data. The blue lines represent the original data, and the dashed red lines denote the shuffled cases.

Figure 4 (a) PDF for degree values (red dots) and the power law fit curve (black line). (b) PDF for weighted degrees (red dots) and the power law fit curve (black line). (c) PDF for edge lengths (km) (red dots) and the power law fit curve (black line).

3 Results

3.1 Characteristics of the PM_2.5 network

The function of positive link weights $W_{i,j}^{\mathrm{pos}}$ and geographical distances D_i,j for the original and the shuffled networks are shown in Fig. 2a. $W_{i,j}^{\mathrm{pos}}$ values in the original network are greater than those in the shuffled network, indicating that the stronger positive links are the result of information transport of PM_2.5 concentrations and the similarity of weather patterns (Liu et al., 2022). For the relation between the largest cross-correlation $P_{i,j}^{\mathrm{pos}}$ versus D_i,j, we observe that the values in the shuffled case are smaller than those in the original case (Fig. 2b), which is in agreement with the pattern of $W_{i,j}^{\mathrm{pos}}$. In the negative case (Fig. 2c and d), there is no distinct difference between the original network and the shuffled network.

Figure 3 shows the probability density function (PDF) of links in the original network and the shuffled network. The PDF of positive link weights has a long tail in the original data, which is not presented in the link weights of the shuffled networks. The PDF of negative link weights is a signature of random behavior, which further indicates that the many significant positive links are not likely to occur by chance. As a result, we consider links that are separated from the shuffled links. Both $W_{i,j}^{\mathrm{pos}}$ and $P_{i,j}^{\mathrm{pos}}$ can be used as a measure of the strength of links. In our analysis, positive link weights of 4.2 are the threshold and accordingly gain the adjacency matrix of the network.

In the network, 284 cities are connected by PM_2.5 concentrations with 3930 links among cities. The clustering coefficient measures the probability that the adjacent nodes of a node are connected. If one city has a high clustering coefficient, there are close connections between its neighbors. In this paper, the clustering coefficient is 0.46. We also analyze the shuffled network with the same number of edges. A total of 337 cities are connected, and the value in the shuffled network is 0.07, suggesting that PM_2.5 cities are more connected to each other. The density of networks is 0.05 in the original network, while the value is 0.03 in the shuffled network. It reflects the degree of completeness of the network, and high values mean strong connections between cities. The average path length is 3.15 and 4.61 for the original and shuffled network, indicating that cities transport the PM_2.5 concentrations to other cities crossed almost three other cities. PM_2.5 cities have a higher clustering coefficient and lower average path length compared with the shuffled network, demonstrating that cities with higher PM_2.5 concentrations can quickly affect their surrounding cities.

Figure 5 Distribution of in-weighted degree (a) and out-weighted degree (b) in the network of each node for positive cases.

The degree of a node is one of the most important statistical properties in networks. The weighted degree characterizes the total strength of correlation of the node with surrounding cities. The PDF of degrees, weighted degrees, and edge lengths of the nodes are shown in Fig. 4. It is found that the degrees, weighted degrees, and edge lengths conform to power law distributions that are associated with some climate and weather phenomena such as the tropical circulations and cyclones (Pierrehumbert, 1986). The power law exponents are 1.3, 1.2, and 1.5, with R² values 0.71, 0.70, and 0.63, respectively. These links are heterogeneous, with few nodes possessing the majority of links in the network. Most of the PM_2.5 concentration links remain confined to a handful of cities.

Figure 6 Distribution of in-weighted degree (a) and out-weighted degree (b) in the network of each node for seasonal data.

Figure 7 (a) Map of PM_2.5 transport links among the monitoring cities in China. (b) The cluster regions of PM_2.5 concentrations. Different colors represent different communities.

To examine a node's dependence or influencing role on the other nodes, we analyze the patterns of in- and out-weighted degrees. The direction of links is determined by the sign of the time delay, which quantifies the incoming or outgoing nodes. Links with zero time delay represent bidirectional links. The in-weighted degree of each node measures incoming links towards the target city, and high values indicate a stronger export effect from source cities to the target city. Out-weighted degrees denote the strength of outgoing links to other cities, and higher values suggest that more cities transfer their PM_2.5 concentrations to the target city. Figure 5 presents the spatial distribution of in- and out-weighted degrees for the whole-year datasets. Different colors represent the ability to transmit. The dark red indicates strong transmission capacity, while white indicates weak transmission capacity. Regions in the Beijing–Tianjin–Hebei–Henan–Shandong (BTHHS) region, YRD, and northwestern China show significant synchronicity with the rest of the provinces in terms of PM_2.5 mass concentrations. These regions correspond to regions with high mean PM_2.5 concentrations. Furthermore, we observe that the distribution of the in-weighted degree is similar to that of the out-weighted degree, which indicates that these cities are both recipients and senders in the networks. This suggests that their pollution is not only due to the local emissions but also imported from other cities. Therefore, solving air pollution should not only rely on reducing emissions in a single city but also on developing inter-city cooperation. Compared with the out-weighted degrees, in-weighted degrees are stronger over the BTHHS region. These cities (sending cities) can also export PM_2.5 concentrations to other cities (recipient cities). In addition, the values of in- and out-weighted degrees display remarkable differences in different seasons, as shown in Fig. 6. The weighted degrees in summer and autumn are small (Fig. 7b and c). In winter and spring, the values of in- and out-weighted degrees are significant (especially in wintertime), and their patterns are similar to that of the whole year.

3.2 Routes and clustering of the PM_2.5

Both in- and out-weighted degrees offer information in terms of nodes (cities). It is reported that urban air quality can be substantially influenced by atmospheric transport of PM_2.5 pollution from distant cities. An analysis of the edges can contribute to revealing the transport routes of PM_2.5 among cities. A recent study found that PM_2.5 concentrations over a distance of 1000 km were related to a typical cyclonic scale within the Rossby waves (Zhang et al., 2019). Here we discuss the transport path within 1000 km and only focus on positive time lags. This is since they are typical links that are related to different climate processes, and they thus enable detailed comparisons with the previous literature. The transport routes show that southward propagation is predominant in the sub-network (Fig. 7a; Zhang and Cao, 2015). We focus on two groups of connections that belong to different regions. The first group is made up of links that traveled from the Gobi Desert over southwestern parts of Mongolia and the Badain Jaran Desert to the BTH region. The second group is made up of links transported from the BTHHS region to the YRD region, and these links show a 1 or 2 d time lag. This is consistent with previous studies obtained from the WRF-Chem model (Huang et al., 2020). The outbreak of YRD pollution usually peaks with a time lag of 1–2 d after that in the BTHHS region.

In addition, we also analyzed the transport routes in different seasons (Fig. 8). The transport routes are significant in autumn and winter, especially in wintertime. It means the route features in winter are dominant over the whole year. Here the southwestern links may be related to the East Asian winter monsoon.

Figure 8 Distribution of transport paths in the network for spring (a), summer (b), autumn (c), and winter (d).

In complex networks, nodes that are closely related to each other are more likely to be grouped in the same cluster. Hence, cities are tightly bound to cities in the same cluster and uncorrelated to cities in other clusters. The pollution transport routes presented above indicate that curbing air pollution is more than just a local issue. In the following, we investigate the cluster features of our networks by utilizing the modularity algorithm described above. Considering a larger Q value means a more accurate community structure for network segmentation, we calculate the Q value at each division to obtain a better result. Here, 284 cities are divided into 9 clusters, where the Q value obtains the maximum value (0.56). The results present a strong regional character, as shown in Fig. 7b. Cities that have the same color represent the same cluster, which could be considered for collaborative governance. These nine regions include the above-mentioned three key regions: the BTH region, YRD region (containing Shanghai, Jiangsu, Anhui, and Zhejiang province), and PRD region (including Guangdong and Guangxi). The other interconnected areas are Heilongjiang and Jilin provinces; Jilin and Liaoning provinces (northeastern China); Hunan and Hubei provinces (central China); Jiangxi and Fujian provinces; Guizhou, Chongqing, and Sichuan provinces; and Shanxi, Shaanxi, Ningxia, and Gansu provinces.

4 Summary and discussion

By constructing PM_2.5 networks based on complex network approaches, it is found that the PDF of the degrees, weighted degrees, and edge lengths of PM_2.5 cities follow a power law distribution, which indicates that the variability of PM_2.5 concentrations in China is not random. Hence, it is reasonable to analyze the transmission and cooperation regions of PM_2.5 from the perspective of the whole national evolution over a long period of time. To quantify the relations of PM_2.5 among cities, the patterns of weighted degrees are investigated. Higher weighted degrees are observed in the BTH region, which is consistent with the patterns of high levels of PM_2.5 concentrations. Cites in the BTH region have stronger strength to export their PM_2.5 pollution to other cities. The distributions of weighted degrees exhibit significant differences in seasons, with the largest being in winter and the smallest being in summer.

Based on the PM_2.5 networks, the transport links and collaborative regions are analyzed. This showed that many links traveled from the Gobi Desert over southwestern parts of Mongolia and the Badain Jaran Desert to the BTH region. The other group extends southward from the BTH region to the YRD region and then further south to Fujian province with a 1 or 2 d time lag. This is consistent with previous studies obtained from the WRF-Chem model (Huang et al., 2014). Although we get a similar transmission pattern in winter, it possesses a stronger intensity. We demonstrate that this is possibly due to the influence of cold fronts, which exactly disperse the PM_2.5 accumulated in the North China Plain to the YRD region, which leads to the propagation of PM_2.5 from the BTH region to the YRD region. Hence, links between the BTH region and the YRD region obtained over the whole year are related to the cold front occurring in wintertime.

In addition, we created a community detection dataset based on the synchronicity of PM_2.5 concentrations and obtained nine clusters. Cities in the same regions should join together to control air pollution. This result provides theoretical support for the Joint Prevention and Control of Air Pollution (JPCAP) proposed by the national government. Regional cooperation should be promoted in these regions to implement regional policies to improve air quality.

A central implication of this study is that the transmission and collaborative regions can be explored via the complex network approach. For traditional model simulation, numerous parameters are needed in the simulation process. In contrast, complex network theory is performed based on time series of field observations, so the estimation process is faster and more economical. As our analysis is based on long-time PM_2.5 records in China rather than a particular region or period of air pollution, it may provide a reference and basis for the development of effective regulatory policies produced by the government to improve air quality. Previous researchers have demonstrated that the accumulated pollutants in the North China Plain can transport the pollution via strong winds to the YRD region based on traditional model simulation, which is similar to our study. We also observed that links that were transported from the BTHHS region to the YRD region show a 1 or 2 d time lag. The result is consistent with previous studies obtained from the WRF-Chem model. Hence, complex network methodologies are useful for the studies of the transport and clustering of air pollutants in faster and more economical ways. Furthermore, there is also potential for the study of other air pollutants, such as ozone and NO_x.

In addition, the study has some limitations. The relationships between PM_2.5 cities have been measured based on the lagged correlations, which have yielded useful results. However, the peak of cross-correlation in a correlogram may be spurious due to serial autocorrelation within each time series, which is another common feature in geophysical time series. Furthermore, the results cannot reveal causal relationships, which may suffer from problems related to interpretability.