Spatial Patterns in the Extreme Dependence of Ozone Pollution between Cities in China’s BTH Region

Abstract

Ozone pollution in China has become increasingly severe in recent years. Considering the damage that extreme ozone pollution may cause and the fact that the occurrence of extreme ozone pollution among different locations may be related, this paper uses hourly ozone concentration data from national monitoring stations to investigate the co-movement of extreme ozone pollution in the Beijing-Tianjin-Hebei (BTH) Region. The extreme dependence analysis is adopted to assess such extreme co-movements between different cities. The co-occurrences of extreme ozone pollution at the same time or with certain time differences in the region are analyzed. City groups suffering simultaneous extreme pollution and those where the pollution occurs with certain time differences are identified under certain criteria. Furthermore, the order in which cities experience extreme ozone pollution is determined. With the publication of the New Three-year Action Plan for Winning the Blue Sky War, our results may be important for improving the joint early-warning and emergency response mechanism at city levels in the BTH Region.

1. Introduction

In recent years, ozone pollution (which refers to ground-level ozone in this paper) has become a growing problem worldwide. According to the State of Global Air 2020 [1], the global exposure to ozone increased from about 47 ppb in 2010 to nearly 50 ppb in 2019 on average, which is 30% to 70% higher than it was 100 years ago. This increase may be due to the rising emissions of the chemicals that form ozone and the increases in global and local temperatures that accelerate ozone formation. Although the national average ozone pollution in the USA decreased by 29% from 1980 to 2021 (Ozone trends by the U.S. Environmental Protection Agency, https://www.epa.gov/air-trends/ozone-trends, accessed on 4 January 2023) because of stringent emission control measures, significant increases in ozone (0.20–0.59 µg m${}^{-3}$ yr${}^{-1}$ for the annual means) are found over the European suburban and urban stations from 1995 to 2012 [2]. Meanwhile, South and East Asia have experienced significant increases in ozone and its precursor emissions and hence have become global ozone pollution hot spots in recent years [3]. Taking China as an example, in 2019, the annual average concentration of ozone in the major 337 cities of China increased by 6.5% compared with the last year [4].

Extreme ozone levels are the most serious and bring great harms that must be paid attention to. A large number of studies have shown that high ozone concentrations can cause a variety of health problems [5,6,7,8] and cause damage to the ecosystem [9,10,11,12,13,14]. In particular, the WHO air quality guidelines [15] report that as ozone concentrations increase above the guideline value, health effects become increasingly numerous and severe. More importantly, due to the transport of pollutants from one location to nearby areas under certain atmospheric conditions [16,17,18,19], the occurrence of ozone pollution among different locations is related. For example, a number of works have studied the regional and long-range transport of ozone and its precursors from one area to the other, showing that the ozone pollution in the area that receives the pollutants is usually significantly increased [20,21]. Under these circumstances, local emission-reduction strategies are not necessarily the most effective [22,23]. Therefore, the co-occurrences or co-movements of extreme ozone pollution at different locations, which are very useful for pollution control and policy-making, are worthy of further study.

Among the few works that address co-movement rather than the extreme co-movement of air pollution [24,25,26,27], spatial autocorrelation analysis has been widely used for centralized pollution data [28,29,30]. However, extreme co-movement, which refers to the correlation between extreme values, is much more complicated to analyze. To describe extreme co-movements among multiple stations, a model builder has to consider many locations and many spatial parameters, which can be hard to implement in a context of multivariate or spatial extreme values and time series [31,32]. Some recent studies have dealt with spatial extremes and processes [33,34,35]. They show that, when modeling extreme co-movement, a fundamental question is whether two series are asymptotically dependent, since for many environmental processes, the spatial dependence weakens as events become more extreme. Therefore, in this paper, the potential spatial co-movement of extreme ozone pollution between different locations was analyzed in the presence of extreme dependence (also called asymptotic dependence or tail dependence).

It is worth noting that the occurrence of ozone pollution in China shows extremes in its wide coverage, high concentration, and long duration [3,36,37,38,39], especially in the Beijing–Tianjin–Hebei (BTH) Region. Pollution in 2020 can be taken as an example. In 28 April 2020, there was a wide range of ozone pollution in the BTH Region and its surrounding areas (28 cities in all). The daily maximum 8-h mean of these areas exceeded 160 µg m${}^{-3}$. In particular, the moving average of the daily maximum 8-h means in Beijing, Tangshan, and Langfang exceeded 240 µg m${}^{-3}$, while the hourly concentration in Beijing even approached 300 µg m${}^{-3}$ on 1 May. Until 2 May, 196 of the 337 largest cities in China demonstrated an excess of 160 µg m${}^{-3}$ of their daily maximum 8-h means. (The WHO air quality guidelines for ozone and its related interim target (daily maximum 8-h mean) are 100 and 160 µg m${}^{-3}$, respectively, [15]. The first limit ensures the public are adequately protected, although some health effects may occur below this level. At the second limit, ozone endangers public health and may have an important impact on health.)

Considering the severity of extreme ozone pollution in the BTH Region and the economic importance of this region in China, this paper investigates the extreme dependencies inside this region. The hourly station-level ozone concentration data from 2014 to 2019 were used. Both simultaneous extreme dependence and extreme dependence with a time difference between each station pair in the BTH Region were considered. Then, city groups suffering simultaneous extreme pollution and those where the pollution occurred with a certain time difference were identified. Furthermore, the order in which cities experienced extreme ozone pollution was determined. According to these results, city-level early-warning and emergency response strategies could be established. The close relationships and time paths among these cities could also be used to enhance the prediction accuracy of an atmospheric pollution model. With the publication of the New Three-year Action Plan for Winning the Blue Sky War in 2018, the above strategies will be important for solving the problem of extreme ozone pollution in China.

2. Extreme Value Selection and Spatio-Temporal Variation of Extreme Ozone Pollution

2.1. Data Source and Approach for Selecting Extreme Values

The hourly station-level ozone concentration data used in this paper were obtained from the China National Environmental Monitoring Center. The numbers of national air quality monitoring stations each year (including stations in use and those still under test) from 2014 to 2019 were 945, 1496, 1498, 1644, 1644 and 1745, respectively. Until 2019, the stations evenly covered the whole country, except for a few provinces in Northwest China, where there were fewer stations because of the low populations in those areas.

Next, we identified the extreme values for each station as those values above a threshold based on a station-specific quantile. We could have used the same threshold for all stations. However, the size of the sample of extreme values would differ too much among stations. Moreover, the extreme values obtained by the approach with a station-specific threshold can capture the seasonality and clustering features of ozone pollution better than the traditional block maxima approach. To ensure that the fitting of the extreme distribution for each station was robust, we used the typical choice of the 90% quantile to obtain extreme values. Using an average of 30 days per month, this threshold represents the minimum of the most severely polluted 72 h of each month for each station. Note that the higher the threshold for a station is, the more serious the extreme ozone pollution for this station is.

2.2. Variation in Extreme Ozone Pollution in China and in the BTH Region

The threshold is the minimum of the extreme values. Figure 1 shows monthly contour maps of the thresholds for stations across China in 2019, which is taken as an example to illustrate the temporal (seasonal) and spatial variation of extreme ozone levels.

From Figure 1, extreme ozone pollution is more severe in the second and third seasons than that during the other two seasons, which is related with the strong ultraviolet radiation and high temperature during this time.

More specifically, from January to March in 2019, the thresholds for all stations across China were below 160 µg m${}^{-3}$. During April to September, which were the worst months of the year, the area with thresholds above 160 µg m${}^{-3}$ had expanded significantly, reaching a maximum in June, as depicted by the areas in orange and red in Figure 1. In October, the areas with higher thresholds gradually moved to lower latitudes, especially coastal areas in the southeast, where temperatures remained high and ultraviolet radiation was still strong. The trends for 2014–2018 are exactly the same.

Regarding the potential co-movements of extreme ozone pollution between different locations, on one hand, it is more meaningful to focus on the second and third seasons, which experience more severe ozone pollution. On the other hand, if two locations are too far from each other, even though there may be a relationship between extreme ozone events, the link is usually too weak to have any practical application [35,40]. Therefore, it is more reasonable to analyze a region with a relatively suitable size. Due to the above points, the BTH Region (as marked with a red circle in June, Figure 1), which experienced a high threshold from April to September, is the focus of this paper. Furthermore, this region includes China’s capital (Beijing) and the international shipping center of northern China (Tianjin). There are very close economic and administrative relations among cities in this region. Exploring the co-movement of extreme ozone pollution in this region could be very important for improving the joint prevention and control strategy for air pollution. The total numbers of stations in this region from 2014 to 2019 were 80, 79, 80, 83, 81, and 87, respectively. These stations belong to 13 cities. Still, the hourly station-level data are used. For simplicity, the number and distribution of the stations in each city are provided in our Supplementary File.

It is worth noting that some values are missing in the original data. On one hand, this is because some of the stations were still being tested at that time. On the other hand, this may be due to temporary hardware errors and data update failures at the time of recording. To deal with this issue, stations with 50% or more missing data points in a month were ignored. Taking 2019 as example, the numbers of stations deleted from the data list were nine stations each in January, February, March, and July, one station each in April and May, two stations in June, and eight stations each from August to December This 50% threshold was a trade-off between overcorrection and undercorrection. A much lower threshold would mean too many stations were ignored, whereas a much higher threshold would mean that too many stations with a large percentage of missing data would be analyzed. In fact, with the 50% threshold, the few stations removed were mostly those with no data points. Removing these stations is obviously better than keeping them, because retaining stations with many missing values would lead to biased estimates. As to the stations left in the data list, their ratios of missing data are mostly lower than 0.1 in most months, except that one station each in August and September has those ratios slightly larger than 0.3. Even so, since this paper focuses on the extreme dependence of ozone pollution between 80 stations on average, the influence of certain missing data at one station in a month is relatively limited, which will be further investigated in our future project.

The boxplots of thresholds among stations in the region were made from April to September in 2015, 2017, and 2019; see Figure 2. The median and mean of the thresholds in each month of 2017 and 2019 were mostly higher than those values in the same month of 2015, showing more severe extreme pollution in 2015. The variances in the thresholds in the months with higher means or medians are larger than those in the months with lower means or medians. This shows that the more severe the extreme pollution in a month is, the larger the difference in thresholds among locations is.

3. The Extreme Dependence Approach

This paper investigates the extreme co-movement of ozone pollution in the presence of extreme dependence. The tail quotient correlation coefficient (TQCC) was calculated to test the extreme dependence. The details of the approach are as follows.

3.1. The Extreme Dependence and the Test Statistics

Two identically distributed random variables X and Y with the distribution function F are extreme-independent (or tail-independent) if

$$\lambda =\underset{u\to x_F}{lim}P(Y>u\mid X>u)$$

(1)

is 0, where $x_F=sup\left\{x\in \mathbb{R}:F\left(x\right)<1\right\}.$ $\lambda$ is called the bivariate tail-dependence index. This measures the dependence of the bivariate upper tails. If $\lambda >0,$ then X and Y have tail dependence; i.e., there are extreme co-movements between X and Y.

The tail dependence given in expression (1) implies the conditional probability of observing extreme ozone pollution at one station given that extreme ozone pollution has already been observed at another station within a time window. Modeling tail dependence has been challenging for years. In [35], the authors introduced an empirically efficient tail quotient correlation coefficient (TQCC) as a statistic to test the null hypothesis of tail independence, namely $H_0:\lambda =0$, against the alternative hypothesis, $H_1:\lambda >0.$ This approach has been proved to be consistent under a much more general condition, which is more applicable in practice.

Given that ${\left\{\left(X_i,Y_i\right)\right\}}_{i=1}^n$ are the series of two Fréchet distributed random variables, the TQCC statistic is

$$q_{u_n}=\frac{\underset{1\le i\le n}{max}\left\{\frac{max\left(X_i,u_n\right)}{max\left(Y_i,u_n\right)}\right\}+\underset{1\le i\le n}{max}\left\{\frac{max\left(Y_i,u_n\right)}{max\left(X_i,u_n\right)}\right\}-2}{\underset{1\le i\le n}{max}\left\{\frac{max\left(X_i,u_n\right)}{max\left(Y_i,u_n\right)}\right\}\times \underset{1\le i\le n}{max}\left\{\frac{max\left(Y_i,u_n\right)}{max\left(X_i,u_n\right)}\right\}-1}$$

(2)

where $u_n$ represents the threshold of $X_i$ and $Y_i$. It is shown in [35] that the TQCC is less sensitive to the choice of $u_n$. The authors also recommended using the 95th percentile (as also applied in this paper).

Since the extreme dependencies are not always strong in all the months in a year, the extreme dependencies for a year are usually relatively small. Therefore, in this paper, TQCC statistics were computed by month to assess the monthly extreme dependence, which can also better describe the seasonality. Specifically, Function (2) was calculated for each station pair in the BTH Region.

The TQCC test statistic, as shown in [35], follows the distribution of a chi-squared random variable with four degrees of freedom under the null hypothesis of tail independence:

$$2n\left\{2-exp\left(-\frac{1}{u_n}\right)\right\}{q}_{u_n}\stackrel{L}{\to }{\chi }_4^2$$

(3)

The null hypothesis is rejected if the test statistic is larger than the critical value for a given level a (the upper a percentile of the ${\chi }_4^2$ distribution); i.e., the tail dependence between X and Y does exist, and its estimate is $q_{u_n}.$

It is important to note that in Equation (2), X and Y are assumed to be unit Fréchet random variables. However, this does not always hold in practice, such that X and Y must first be transformed separately into a unit Fréchet distribution. One way to achieve this is to fit the GEV distribution (as detailed in Section 3.2) for each station first, obtain the estimate function $\widehat{F}\left(x\right)$, and then use the transformation $-1/log\left(\widehat{F}\left(x\right)\right),$ which is the unit Fréchet scale. Ref. [35] showed that the TQCC calculated using the above transformation with the estimated parameters of the GEV distribution converges into the same asymptotic distribution with the true parameters.

3.2. Fitting the Generalized Extreme Value Distribution Based on a Point Process Approach

The function of the GEV distribution is as follows:

$$F\left(x\right)=exp\left\{-{\left(1+\xi \frac{x-\mu }{\sigma }\right)}^{-1/\xi }\right\}$$

(4)

where $1+\xi (x-\mu )/\sigma >0$. Here, $\mu$ is the location parameter and $\sigma$ is the scale parameter. $\xi$ is the shape parameter. When $\xi >0,$ $F\left(x\right)$ is a Fréchet distribution function, its density has a lower bound $-\sigma /\xi +\mu$, and it is long-tailed and right-skewed. When $\xi <0,$ $F\left(x\right)$ is a Weibull distribution function, its density has an upper bound $-\sigma /\xi +\mu ,$ and it is short-tailed and left-skewed. In the case of $\xi \to 0,$ the distribution function is interpreted as the limit of Function (4), leading to a Gumbel distribution function (omitted for simplicity), with no bounds for its density.

The point process approach is used to fit Function (4), which is more robust to time dependence and potential non-stationarity [41,42]. In this approach, the extreme values above a specified threshold, rather than the traditional block maximas, are used. Then, the times at which exceedances over a certain threshold occur and the excess values over the threshold are combined into one process. This shows that with a suitable normalization, this process behaves like a nonhomogeneous Poisson process, the likelihood function of which can be expressed as follows:

$$\begin{array}{c}\hfill L=exp\left\{-n{\left[1+\xi \left(\frac{u-\mu }{\sigma }\right)\right]}^{-1/\xi }\right\}\\ \hfill {\displaystyle •}{\displaystyle \prod _{i=1}^N}\left\{{\sigma }^{-1}{\left[1+\xi \left(\frac{y_i-\mu }{\sigma }\right)\right]}^{-1/\xi -1}\right\}\end{array}$$

(5)

where $\sigma$, $\xi$, and $\mu$ are the parameters of the GEV distribution, Function (4). n is the desired number of extreme values for the GEV distribution. For example, in the case that n is the number of days in a month, the estimated parameters correspond immediately to the GEV parameters of the daily maximum distribution of the observed process. $y_i$ $(i=1,2,\dots ,N)$ are the hourly extreme values (above threshold u) in each month, which are used in the estimation, where N is the total number of extreme values.

To match the calculation of monthly TQCC, the GEV distribution should be fitted by month for each station. A typical 90% quantile of each month, as in Section 2.1, was used to select the extreme values which can also guarantee a sufficient sample size and increase the fitting accuracy. Note that this threshold is different with $u_n$ in Function (2). While the threshold in this section is for fitting the GEV distribution via the point process approach, $u_n$ determines the extremes of co-movement in Function (2). These values are not necessarily the same.

The parameters were estimated by maximizing the log of Function (5). In practice, these maximum likelihood estimators are local maxima found by nonlinear optimization. In [43], the authors recommend quasi-Newton iterations for optimizing and obtaining the starting values by fitting the Gumbel distribution. All of the estimators are asymptotically normal. Since the GEV distribution is fitted for each station separately, the significance test has to be conducted in the same way. However, since the total number of stations is very large, these single tests may result in a significant increase in the total number of false rejections of the null hypothesis. Thus, a practical approach, i.e., multiple testing [44,45], is used to derive a more accurate inference, in which the false discovery rate is controlled at 0.05.

4. Co-Movement of Extreme Ozone Pollution in the BTH Region

4.1. TQCC Approach with Time Differences

When calculating extreme dependencies between stations, there may exist a certain time difference d (hours) between the start times of the two datasets. For two stations, $s_i$ and $s_j$, if the first extreme value for station $s_i$ is earlier than that for station $s_j$, then extreme ozone pollution occurred earlier at station $s_i$ than at station $s_j$ (denoted by $d>0$), and vice versa. We set the maximum time difference as $d_{max}$ ($d_{max}>0$); then, considering the integers in $[-d_{max},d_{max}]$, there are $2d_{max}+1$ cases altogether. $d=0$ is the case for simultaneous extreme dependence, and TQCC, calculated by Equation (2), is denoted by $q_0$. It is also important to know the largest extreme dependence among all cases (the related TQCC is represented by $q_{max}$) and the related d (represented by $d_{biggest}$). Here, $d_{biggest}$ can also be interpreted as the time difference between when one station experiences extreme ozone pollution and when such pollution would subsequently occur at another station with high probability.

4.2. Monthly Variation in TQCC

With regard to the monthly TQCCs calculated in this paper, the maximum time difference is set as 48 h. This is because the extreme values used to calculate TQCC are above the $90\%$ quantile, which means that we consider only $10\%$ of the time in a month to have experienced extreme pollution. Thus, the longest period of extreme pollution is naturally 3 days at most. If the maximum time difference is set to 48 h, then d is an integral number of hours ranging from $-48$ to 48 (spanning 4 days), which covers the longest duration of an extreme pollution event (3 days).

The TQCC between every two stations in the BTH Region (around 80 stations) was computed for different values of d. Since TQCC is symmetric for station $s_i$ and station $s_j$, as shown in Equation (2), TQCC for station pairs can be represented by a triangular matrix. Figure 3 and Figure 4 show triangular matrices for $q_0$ and $q_{max}$ for each month in 2019, respectively.

Figure 3 and Figure 4 show that both $q_0$ and $q_{max}$ increase in each month until October and then decrease a little. These monthly trends for TQCC do not match the seasonality shown in Figure 1. Extreme ozone pollution was relatively severe from April to September due to the strong ultraviolet irradiation. However, the strong TQCC between stations could be due to the transmission of pollutants under certain atmospheric and meteorological conditions, which are not the causes of the extreme pollution itself.

The results for the two months are interesting. First, in August, the extreme dependence in Figure 3 is relatively high, since the median and mean of $q_0$ are close to $0.5$. Furthermore, the values of $q_0$ in Figure 3 and $q_{max}$ in Figure 4 are similar in this month for each pair of stations. This shows that even when considering the time difference, the extreme dependence in August is mainly the strongest for $d=0$. Therefore, August is an important month for analyzing the simultaneous extreme dependence, especially since the threshold shown in Figure 1 is also high.

Furthermore, although $q_0$ in September is relatively small (the median and mean of $q_0$ are both below $0.4$; see Figure 3), $q_{max}$ is relatively high (the median and mean of $q_{max}$ are both around $0.6$; see Figure 4). The large difference between $q_0$ and $q_{max}$ in September shows that there is a time difference in the occurrence of the extreme dependence between stations. Therefore, September, which also experiences severe extreme pollution, as shown in Figure 1, is an important month for analyzing the time difference. These monthly variations were similar in other years.

4.3. Cities with Simultaneous Extreme Dependence of Ozone Pollution in the BTH Region

We computed the simultaneous extreme dependence between every two stations for August. However, the relationships are not clear since there are too many station pairs in the BTH Region. Therefore, to assist the government with its policy making, we identified cities with strong simultaneously extreme dependence with the following steps:

Step 1: TQCC selection. Select all values of $q_0$ greater than or equal to a threshold $q_{Tsh}$.

Step 2: City identification. For each selected value of $q_0$, identify whether the two stations belong to the same city or two different cities (a city pair).

Step 3: Ratio calculation. For each city pair, calculate the ratio of the number of selected $q_0$ values to the total number of station pairs inside the pair.

Step 4: City selection. Select all city pairs that have a ratio greater than or equal to a threshold $r_{Tsh}$. These selected city pairs have enough station pairs suffering extreme ozone pollution simultaneously with very large probability.

Table 1. Results of selected city pairs for different $q_{Tsh}$ and $r_{Tsh}$.

	${\mathit{q}}_{\mathit{T}\mathit{s}\mathit{h}}$	${\mathit{r}}_{\mathit{T}\mathit{s}\mathit{h}}$	Results
Number of selected city pairs (number of cities involved)	0.4	0.5	61 (13)
		0.6	57 (13)
		0.7	54 (13)
Number of selected city pairs (number of cities involved)	0.5	0.5	31 (12)
		0.6	26 (12)
		0.7	23 (11)
Number of selected city pairs (number of cities involved)	0.6	0.5	8 (8)
		0.6	3 (5)
		0.7	2 (3)

Note: There are 13 cities in the BTH Region. Since TQCCs ≥ 0.7 is very rare between stations, the maximum value of $r_{Tsh}$ is set to 0.7.

shows the numbers of city pairs and cities involved for different values of $q_{Tsh}$ and $r_{Tsh}$.

The number of selected city pairs and the number of cities involved mainly decrease with $q_{Tsh}$ but not with $r_{Tsh}$. The strong extreme-dependence case ($q_{Tsh}=0.6$, $r_{Tsh}=0.5$), which has eight selected city pairs, is analyzed to give insights into the simultaneous extreme dependence between cities.

There are two interesting city groups for which any two cities are a selected city pair: Langfang–Baoding–Shijiazhuang and Langfang–Tianjin–Shijiazhuang. Considering that the ratio of selected $q_0$ between Tianjin and Baoding was 0.33 (which was below our threshold $r_{Tsh}=0.5$), it could be concluded that, in general, any two cities in Langfang, Tianjin, Baoding, and Shijiazhuang (C8–C2–C7–C9 as shown in Figure 5) satisfy the hypothesis that if any city in the group experiences extreme ozone pollution, then the other three cities will also experience the same simultaneously with high probability. These four cities are called a simultaneous city group.

These four cities are all in the middle of the BTH Region. Their close geographical location is one of the key reasons for the co-occurrence of extreme ozone pollution.

4.4. The Time Difference of Extreme Dependence of Ozone Pollution in BTH Region

This section focuses on the extreme dependence of ozone pollution with a time difference for cities in BTH Region. As in the previous section, the following steps are used to select city pairs with extreme dependence among stations.

Step 1: $d_{biggest}$ selection. Select all values of $d_{biggest}$ that fall into the intervals given by $d_{Tsh}$.

Step 2: City identification. Identify the city pair that the stations for each selected value of $d_{biggest}$ belong to for each interval of $d_{Tsh}$.

Step 3: Ratio calculation. Calculate the ratio of the number of selected values of $d_{biggest}$ to the total number of station pairs inside every city pair for each interval of $d_{Tsh}$. The ratio is the probability of the occurrence of extreme pollution within the time difference.

Step 4: City selection. Select all selected city pairs that have a ratio greater than or equal to the threshold $r_{Tsh}$. These city pairs have a large probability that extreme pollution will occur within a certain time difference.

Since the values of $d_{biggest}$ are mostly clustered within $[-24,24]$ out of $[-48,48]$ and considering the importance of an inter-day specific analysis for an emergency response, $d_{Tsh}$ is set to $[-12,12]$, $[-24,-12)\cup (12,24]$, and $[-48,-24)\cup (24,48]$. These ranges ensure that enough $d_{biggest}$ values fall into each time interval for analysis.

Table 2. Results of selected single cites and city pairs for different $d_{Tsh}$ and $r_{Tsh}$.

	${\mathit{d}}_{\mathit{T}\mathit{s}\mathit{h}}$	${\mathit{r}}_{\mathit{T}\mathit{s}\mathit{h}}$	Results
Number of selected city pairs (number of cities involved)	$[-12,12]$	0.6	11 (10)
		0.7	10 (9)
		0.8	7 (7)
Number of selected city pairs (number of cities involved)	$[-24,-12)$$\cup (12,24]$	0.6	6 (11)
		0.7	3 (5)
		0.8	2 (3)
Number of selected city pairs (number of cities involved)	$[-48,-24)$$\cup (24,48]$	0.6	34 (13)
		0.7	22 (13)
		0.8	16 (12)

Note: There are 13 cities in the BTH Region.

gives the results for different combinations of $d_{Tsh}$ and $r_{Tsh}$ in September 2019. The higher the value of $r_{Tsh}$, the higher the probability that extreme pollution will occur within a certain time difference. We use $r_{Tsh}=0.8$ to illustrate the results.

For $d_{Tsh}=[-12,12]$, any two cities among Tangshan, Baoding, and Qinhuangdao (C6–C7–C5, as connected by the blue lines in Figure 5) are the selected city pairs. Therefore, those three cities are called the 12-h city group, since if extreme pollution occurs in any one city in the group, it will also happen in the other two cities within 12 h.

Note that interval $[-12,12]$ includes $d=0$. However, the simultaneous city group is not the same as the 12-h city group. This is because for the simultaneous city group, extreme ozone pollution happens in the cities at exactly the same time, in which case $r_{Tsh}=0.5$. However, for the 12-h city group, $r_{Tsh}$ is much higher (0.8), which ensures a high ratio of station pairs falling into this time interval; i.e., the 12-h city group is more precise.

The number of selected city pairs decreased considerably when $d_{Tsh}=[-24,-12)\cup (12,24]$. There is no city group in this case, unlike in the 12-h city group. When $d_{Tsh}=[-48,-24)\cup (24,48]$, the selected city pairs include almost all cities in the BTH Region, showing universal extreme dependence among cities during an extreme pollution event. For some cities, if extreme pollution happens in any one city, it will spread to the other cities within 24 to 48 h. We name these cities the 24-to-48-h city groups. A typical example is the Hengshui–Langfang–Zhangjiakou group (C10–C8–C3, as connected by the orange lines in Figure 5), which relates to the highest ratio among all the other city pairs.

Note that there are six different cities in the simultaneous city group and 12-h city group (C2–C5–C6–C7–C8–C9 in Figure 5), which are aligned down the middle of the BTH Region from west to east. Since there is one overlapping city in the two groups, C7, if any one of these six cities experiences extreme pollution, no matter which group it belongs to, the other five cities will experience it very soon (within 12 h or even simultaneously).

In contrast, the typical 24-to-48-h city group is arranged from south to north (C3–C8–C10 in Figure 5), which has a longer time difference for the co-occurrence of extreme pollution than for the six cities arranged from west to east. This is probably because the Taihang mountain acts as a barrier between C3 and C8–C10 in the south, which may potentially delay the transmission of pollution.

It is interesting that C8 falls into both the west–east cities and the south–north cities. Since the west–east cities usually experience extreme pollution within 12 h, once any one of these cities experiences extreme pollution, C8 may suffer very soon after, and then the same will happen to the south–north cities within 48 h.

4.5. Sequential Order for the Occurrence of Extreme Ozone Pollution

As mentioned in Section 4.1, $d>0$ means that the occurrence of extreme ozone pollution at station $s_i$ is d hours earlier than at station $s_j$, whereas $d<0$ means the opposite and $d=0$ means it occurs simultaneously. Therefore, the sign of $d_{biggest}$ can be used to determine the sequential order for the occurrence of extreme pollution in the two cities with the highest probability.

We focus on the monthly $d_{biggest}$ in the second and third seasons. These are the worst two seasons, so it is more meaningful to study them. All the values of $d_{biggest}$ for each city pair were calculated. Then, the ratios for three cases were obtained: $d_{biggest}>0$, $d_{biggest}<0$, and $d_{biggest}=0$. Finally, the sequential order of each pair of cities was found from the sign for the largest value among those three ratios.

There are some special results for Beijing and Cangzhou. The occurrence of extreme ozone pollution in Beijing usually appears after the cities to its south during April and May, but is more likely to happen otherwise (earlier than most of these cities), especially in September. Pollution is more likely to occur later in Cangzhou than the cities around it from June to August. The possible reasons for these results and their implications are as follows.

During April and May, the level of extreme ozone pollution for all cities has not yet reached its maximum and is relatively low overall. This is shown in Figure 2, in which the thresholds for the cities during April and May are lower than in June. Meanwhile, the differences in the ozone concentrations among cities are probably not large enough (evidenced by the smaller variance of the threshold, as shown in Figure 2) to drive the strong transmission of pollutants from high to low concentration areas. The direction of transmission and diffusion of pollutants is more likely to be affected by the southerly wind in spring. Thus, the occurrence of extreme ozone pollution in Beijing is relatively later than in most of the cities to its south (C2, C7–C13; Figure 6).

June is the start of summer, and the wind speeds drop everywhere in the BTH region. Therefore, the diffusion and transmission of pollutants due to the wind becomes weaker. Moreover, pollution accumulates rapidly due to the high temperatures and low wind speeds in summer. This results in an increase in the extreme ozone pollution level and, more importantly, an increase in the differences in levels among cities, which is evidenced by the larger variance in thresholds shown in Figure 2. Specifically, in June, the monthly threshold of extreme pollution in Cangzhou is 201.3 µg m${}^{-3}$, which is far lower than the threshold in cities C1, C2, C7–C10, and C12–C13 (for the location of these cities, please refer to Figure 6), whose thresholds are 214.8, 210.3, 223.7, 203.3, 231, 213, 230, and 229 µg m${}^{-3}$, respectively. The potential concentration difference among these cities probably causes the pollutants to spread from the cities with high levels (those mentioned above) to cities with lower levels, such as Cangzhou. It seems that pollution levels in Cangzhou are more likely to continue rising for some time due to the transmission, although other cities have experienced their own extreme values. Therefore, the occurrence of extreme pollution in Cangzhou tends to be later than in other cities from June to August.

Since September is in fall, both the illumination intensity and temperature decrease, such that the meteorological conditions are no longer optimum for the formation of ground-level ozone. Therefore, the extreme pollution levels decrease widely in the region. The wind becomes stronger than in summer, which may facilitate the transmission and diffusion of pollutants. Therefore, the prevailing north wind probably results in the earlier appearance of extreme pollution in northerly Beijing compared with other cities (C2, C7, and C10–C13; the locations of these cities are as shown in Figure 6).

5. Discussion and Conclusions

Considering the severity of extreme ozone pollution in the BTH Region, as well as the fact that the occurrences of extreme ozone pollution in different locations are usually related events rather than independent, this paper adopted the extreme dependence approach to study the extreme co-movements of ozone pollution between cities in this region. The hourly ozone concentration data from national monitoring stations from 2014 to 2019 were used. City groups suffering simultaneous extreme pollution and those where the pollution occurs with certain time differences were identified under certain criteria. Furthermore, the order in which cities experience extreme ozone pollution was determined. Possible reasons resulting in the above spatial pattern were also given, which are the diffusion and transmission of ozone pollution and the concentration differences among different locations.

Similar to the co-movement of ozone pollution focused in this paper, previous works have investigated the regional interaction of ozone in the BTH Region caused by the transport of pollutants. For example, some studies show that more than 30% of the contribution to ozone mass burdens in Beijing and Tianjin comes from the emission sources in Shandong province and Hebei province [46]. Others show that atmospheric transport can largely modulate regional interactions of ozone pollution in China [20]. In particular, ozone transport from Central and Eastern China contributes 36% of the enhanced ozone concentrations in North China (including BTH Region) [36]. The above works support the existence of the co-movement of ozone pollution in our paper. Furthermore, [47] states that the highest ozone concentrations appear along the Taihang Mountains, which is consistent with our results on the reason for the longer time difference of the 24-48 city groups than the simultaneous and 12-h city group.

However, in contrast to the previous work, firstly, this paper focuses on extreme ozone pollution rather than the mean ozone concentration [28,30,37,48]. Since a mean is the average level of ozone concentration over a period of time, information about extreme values, which cause the most damage, is lost. Consequently, an analysis using average values may underestimate the harm and impact of ozone pollution. Therefore, the extreme co-movement of pollution, rather than just the co-movement, was investigated in this paper. Secondly, this paper focuses on extreme dependence, i.e., the probability of the occurrence of extreme ozone pollution at one location given that extreme ozone pollution has been observed at other locations, which focuses on statistical results, rather than meteorological reasons, etc. Of course, the factors that influence the co-occurrence of ozone pollution are important, but they are usually complicated to model, and thus will be further investigated in our future work. The analyses of the extreme dependence in our paper are also powerful and can be used for early warning, pollution prediction, joint prevention and control, etc.

Specifically, considering the close co-movement among cities in the simultaneous city group and 12-h city group, we suggest that these cities should adopt the joint early warning and response mechanism. Once any one of these cities is going to experience a certain high ozone concentration, all of the cities can activate their response measures at the same time to prevent the occurrence of extreme pollution over a wide area.

In contrast to the above cities, arranged from west to east, the 24-to-48-h city groups, arranged from north to south, have longer time differences. It is possible that these time and path differences among these cities could be used to build an atmospheric pollution model to enhance prediction accuracy.

With regard to the special sequential order of Beijing and Cangzhou in the occurrence of extreme pollution, we suggest using Cangzhou as an early-warning city from June to August and using Beijing as an early-warning city in September, considering that, in these two cities, extreme ozone pollution often occurs earlier than in other cities during these periods.

Despite the above descriptions, theoretical or empirical analysis on how extreme dependence can be applied to early warning, prediction, and the joint prevention and control of pollution is still required. Moreover, modeling the extreme relationships among multiple locations, especially with covariates, is still a challenge. We will explore these issues in our future research.