Association between NO₂ concentrations and spatial configuration: a study of the impacts of COVID-19 lockdowns in 54 US cities

The NO₂ concentration records before lockdown will be used to train historical patterns that contain seasonal and cyclic time-series information, and the SARIMAX model will be used to predict the NO₂ concentration. In this case, a prediction following historical patterns will not integrate with the new evolution after lockdown measures. Under the circumstance, two definitions are proposed to investigate the evolution:

• NO₂ difference denoted by Δd calculates the difference between the observed and predicted NO₂ concentration at a station after lockdown; and
• NO₂ change denoted by Δc calculates the difference of the mean NO₂ concentration before lockdown and that after lockdown based on observation or prediction.

Particularly, Δd can be significant as a result of the changes of urban mobility during the implementation of lockdown measures. Thus, we can investigate the relationship between Δd and the changes of urban mobility (Δm) subject to a baseline before lockdown.

The general ARIMA model consists of three parts, namely autoregression (AR), differencing (I) and moving average (MA), and the model can be presented as follows:

$$\begin{equation} \textrm{ARIMA}(p,d,q), \end{equation} \tag{ 1 }$$

where p denotes the order of auto-regressive model, d represents the order of difference, and q is the order of moving average model. Nonetheless, seasonal patterns cannot be established by introducing an ordinary ARIMA model. To address this problem, the seasonal auto-regressive integrated moving average (SARIMA) includes the seasonal parameters for the auto-regressive, differencing and moving average terms, and the time step for modelling the periodic time-series pattern:

$$\begin{equation} \textrm{ARIMA}(p,d,q)(P,D,Q)_t, \end{equation} \tag{ 2 }$$

where P denotes the order of seasonal auto-regressive model, D denotes the order of seasonal differencing, Q is the order of seasonal moving average model, and t refers to the time step parameter for seasonality. Furthermore, Seasonal Autoregressive Integrated Moving Average with eXogenous factors (SARIMAX) can be used to analyse the time-series data with additional variables in the regression operation (Valipour 2015):

$$\begin{equation} \textrm{ARIMAX}(p,d,q)(P,D,Q)_t. \end{equation} \tag{ 3 }$$

In this study, the time-series NO₂ measurements of 56 stations in the US were modelled using the SARIMAX model based on the exogenous weather variables including air temperature, air pressure, relative humidity, wind speed and percentage of cloud coverage. Since there might be more than one seasonal pattern for the time-series NO₂ data, six Fourier series with different periods were also included as exogenous variables (number of days = {15, 30, 90, 180, 365}), and 7 days (one week) was deemed to be the seasonal time step t in the SARIMAX model.

Before the SARIMAX modelling is performed, the randomness and normality were tested with Ljung–Box test and Jarque–Bera test. The Augmented Dickey–Fuller test and the Osborn, Chui, Smith, and Birchenhall tests were performed to determine the seasonally varying parameters (d and D) of the time series respectively. Canova-Hansen test was conducted to determine the time step parameter (t) for seasonality. Then, a stepwise approach (Hyndman et al 2008) was adopted for optimizing the other parameters (i.e. p, q, P, and Q) of the SARIMA model by minimizing the Akaike information criterion (AIC) value. The optimized parameters were introduced for the NO₂ prediction in this study. Python 3.7 and the statsmodels package (State Space Methods 2020) were used for training and prediction.

As traffic is one of the major sources of NO₂, in the present study, four indices are adopted to explain the characteristics of road networks in a 3 km-radius circular area centred at each air quality monitoring station (figure 1(a)), which are used for correlation analysis between NO₂ difference after lockdown. The four indices are enumerated as follows: (1) n_nd denotes the number of road intersections, (2) c_nd presents the number of the road segments at all the intersections, (3) s_rd stands for the total length of the road networks, and (4) n_rd is the number of the road segments. It is crucial for proposing n_nd since vehicles often stop at the intersections because of the traffic lights, which leads to the accumulation of a considerable amount of the emissions during idling (Minoura et al 2010). In comparison, c_nd shall be more representative than n_nd because the number of roads connecting at each intersection is counted that may indicate temporary stopping more confidently. For example, vehicles are more likely to be parked at an intersection of two roads than a single road with the same traffic volume. s_rd also plays a big part since it affects traffic flows fundamentally. The four indices are thus organized as $i = \{n_{nd}, c_{nd}, s_{rd}, n_{rd}\}$. In the study, G denote a topological graph of roads that edges E = {e} connect with each other by associating with the nodes O = {o}, getting G = {E, O}. Notably, n_nd is different from the number of nodes denotes by num(O) because only two edges connecting through a node means they are the same road essentially, which is caused by the data format. Thus, n_nd and c_nd are counted when a node at least associates with at least edges.

Zoom In Zoom Out Reset image size

However, based on the four indices, it is estimated that the dispersion of vehicle emissions at road networks has homogeneous impacts on stations' air quality when the spatial distance between them is taken into consideration. To make a better estimation, in the present study, it is assumed that the dispersion of each index follows a Gaussian distribution from each element of the index to the station (figure 1(b), equation (4)):

$$\begin{equation} y = y_0 + Ae^{-\frac{(x-xc)^2}{2w^2}}. \end{equation} \tag{ 4 }$$

In the function, y₀ denotes the offset and xc is the centre, both equalling 0 to present a normalized Gaussian. A is the amplitude to denote the magnitude of the element, w implies the width to control the speed of the dispersion, and x stands for the distance from the element to the station (figure 1(c)). When computing n_rd and s_rd enriched with the Gaussian function, x is the distance from the middle point of the road segment to the station, as an approximation. Then, the total impacts of all the elements of an index at a station can be accumulated (equation (5)):

$$\begin{equation} Y = \sum y = \sum_{m = 1}^{n} g(A_m, x_m, w). \end{equation} \tag{ 5 }$$

A set of abbreviations and their definitions are summarized in table 2 for clear presentation of this study.

Daily NO₂ observations at 56 stations in 54 cities across the US are derived from the Air Quality Open Data Platform during a period from June 2014 to May 2020 (AQODP 2020), which are complemented by the geographical coordinates. The weather data taken as exogenous factors in the SARIMAX model are collected from OpenWeatherMap (Open Weather 2020) with the coordinates of the 56 NO₂ stations. The COVID-19 Community Mobility Reports from Google provides daily mobility changes in percentage subject to the baseline on 15 February 2020 (Google 2020), indicating user mobility in Google Maps during the COVID-19 pandemic. Mobility is aggregated and anonymized as six categorical purposes, including retail and recreation, groceries and pharmacies, parks, transit stations, workplaces, and residences. The originally collected NO₂ data and the mobility data have the same temporal resolution on one day, so that preprocessing is not needed. The study is based on the data collection between 15 February 2020 and 6 May 2020, which is the last day of the NO₂ data. Road networks in the US are acquired from OpenStreetMap (2020).

Figure 2 presents 56 NO₂ stations in the US, covering a large geographical space, and table 1 summarizes six climate types that prevail at the stations according to National Weather Service (2020). The stations are located in various micro-environments with different road networks, and they are located in urban or rural areas of those cities that implement lockdown when the COVID-19 pandemic becomes quite severe. The training has been done at each station since June 2014. Then, the SARIMAX model is used to predict daily NO₂ at each station from the start of lockdown to the end of lockdown.

Zoom In Zoom Out Reset image size

Table 2. Definitions of the abbreviations used in this study.

No.	Abbre.	Definition
1	Δd	NO2 difference: the predicted NO2 concentration minus the observed NO2 concentration at a station after lockdown
2	e	relative error: Δd divided by the observed NO2 concentration
3	Δc	NO2 change: the mean NO2 concentration after lockdown minus that before lockdown
4	Δcp	NO2 change in percentage: Δc divided by the mean NO2 concentration before lockdown
5	Δm	the mean change of urban mobility subject to a specific day before lockdown; {Δmw , Δmr , Δmt } are the changes of mobility for workplace, recreation, and public transit
6	nnd	the number of road intersections
7	cnd	the number of the road segments connecting at all the intersections
8	srd	the total length of the road networks
9	nrd	the number of the road segments
10	Δmrd	$\Delta m_{rd} = s_{rd}\cdot \Delta m_w$

Two predictions are made, with one prediction period starting 60 days (n = 60) before the lockdown and ending a number of days after the lockdown (until 6 May 2020), and the other prediction period is starting 30 days (n = 30) before the lockdown and ending a number of days after the lockdown. Figure 3 visualizes prediction results of two selected sites that the red curves (observation) present a seasonal and periodical pattern, which are lower than the blue curves (prediction). The difference between the blue and red curves (i.e. NO₂ difference) confirms our hypothesis that the prediction has not incorporated the disruption of lockdown and can be explained by the reduction of mobility.

Zoom In Zoom Out Reset image size

Figure 4(a) provides the distribution of the mean relative errors (e) of the prediction at all the stations. Evidently, the errors before lockdown (b) are smaller than after lockdown (a), with the medians at e = 0.13 for n = 60 and e = 0.12 for n = 30, respectively. It suggests that satisfactory prediction accuracy is achieved prior to the lockdown accordingly. When the lockdown was activated, the medians become larger at e = 0.22 for n = 60 and e = 0.20 for n = 30, which indicates that it is challenging to obtain accurate prediction caused by disruptive lockdown measures. Besides, both errors for n = 30 are slightly smaller than that for n = 60, probably because n = 60 has a longer prediction time while n = 30 allows an additional 30 day training that better incorporates the most recently seasonal variation into the prediction.

Zoom In Zoom Out Reset image size

To probe into the impact of the lockdown measures, in the present study, the NO₂ changes before and after lockdown for the 56 stations based on prediction (p) and observation (o) are calculated, respectively. It is found that many stations have undergone a substantial reduction after lockdown in view of the changes of absolute values (Δc) (figure 4(b)), which become more significant when conversion into percentages that Δc_p equals Δc divided by the mean concentration before lockdown (figure 4(c)). Specifically, the observation suggests that half of the NO₂ measurements decreased at least with −1.45 mg m⁻³ for n = 30 and −2.41 mg m⁻³ for n = 60 (figure 4(b)), which is equivalent to a 23.08% and 30.38% reduction of the NO₂ concentration before lockdown. Meanwhile, the predicted number is slightly larger than the observation.

According to the COVID-19 Community Mobility Reports from Google that provides relative urban mobility subject to a baseline before lockdown in each city, all the involved cities have shown a huge decrease in the mobility for workplace with $|\Delta m_w|\gt20\%$, in which Seattle, Ashburn, Burlington, and Minneapolis have witnessed the largest reduction with $|\Delta m_w|$ equalling 41%, 38%, 36%, and 33%, respectively (figure 5(a)). By inspecting each station, the observation suggests that a considerable reduction of NO₂ after lockdown is realized in most cities, while a tiny increase is found in four cities quite unexpectedly (figure 5(b)). Notably, the changing trend of $|\Delta m_w|$ is unclear with the decrease of Δc, which stimulates our motivation to explore the influential factors.

Zoom In Zoom Out Reset image size

Considering the vehicle emission is one of the largest resources of NO₂, a significant reduction of urban mobility during the pandemic is supposed to have impacts on the changes of NO₂ in cities, which focuses on a detailed investigation by contrast with other studies at regional and national scales. Based on the COVID-19 Community Mobility Reports, the Pearson correlation analysis between the NO₂ difference (Δd) and the averaged mobility changes (Δm) is made, in categories of n = 30 (figure 6(a)) and n = 60 (figure 6(b)). It shows that the mobility for workplace decreases significantly in all cities ($-40\% \leqslant \Delta m_{w}\leqslant -15\%$), which validates a moderate and negative correlation with Δd in both groups (the Pearson correlation coefficient R =−0.47, $p\lt$0.05). Similarly, the mobility for recreation has shown a considerable decrease in a vast majority cities ($-30\% \leqslant \Delta m_{c}\leqslant 1\%$), having a weak and negative correlation with Δd (R =−0.37 for n = 30 and R =−0.40 for n = 60). In comparison, the greatest decrease in transit mobility is observed in almost all the cities ($-50\% \leqslant \Delta m_{t}\leqslant 5\%$), which contributes to a weak and negative correlation (R =−0.32 for n = 30 and R =−0.34 for n = 60) (figure 6).

Zoom In Zoom Out Reset image size

As an expectation, a small increase of Δm_c = 0.94% and Δm_c = 5.95% is revealed in the recreation mobility for the Cedar Bluff State Park in Kansas and the transit mobility for Murphy Ridge in Wyoming, respectively. It is reasonable since some populations may deliberately choose not to travel to crowded locations and move to rural areas as temporary residences such as parks. Three implications can also be induced by the results. Firstly, the lockdown policy helps reduce urban mobility for recreation, transit, and workplace greatly during the pandemic, leading to a considerable reduction of NO₂. Secondly, the correlations between Δm and Δd, even though not robustly significant, suggest the validation of the hypothesis that the difference of NO₂ has resulted from the changes in urban mobility. Last but not least, other factors such as road characteristics or local climates may also affect the changes of NO₂ since the current correlations are not significantly strong.

As the urban mobility restraint by road networks, in the present studies, it is proposed that the road characteristics, such as the density of road networks, may have notable impacts on emission. Therefore, three static indices are put forth to represent road characteristics $i = \{n_{nd},c_{nd},s_{rd}\}$: (1) the number of road intersections n_nd , (2) the total number of the road segments connecting at all the intersections c_nd , and (3) the total length of the road networks s_rd . Figure 7 demonstrates the correlations between i and Δd, categorizing into two groups of n = 30 and n = 60. It is pointed out that both groups obtain somewhat strong and positive correlations that the coefficients R = {0.54, 0.54, 0.60} for n = 30, resulting a better performance than (Δm, Δd) in figure 6. It can be explained that a larger indicator in i makes a more significant effect on the changes of NO₂ after lockdown, suggesting that road networks' characteristics can greatly affect Δd. However, the three indicators treat intersections and roads with homogeneous impacts on the air quality stations, regardless of the spatial distance between the roads and stations.

Zoom In Zoom Out Reset image size

To make the indices more representative, we used a Gaussian function y = ∑ g(A, x, w) to estimate the dispersion of NO₂ from the vehicle emission, in which A represents the magnitude of the road property, x denotes the distance from roads to the NO₂ station, and w indicates the width of the function. In addition, another road index named n_rd is initiated, which means that the dispersion associated with the road follows the Gaussian distribution but ignoring the characteristics of the road length to compare with the performance of s_rd . Figure 8 shows the coefficient curves between $i = \{n_{nd},c_{nd},s_{rd},n_{rd}\}$ and Δd, comparing between n = 30 and n = 60. Overall, all the curves grow steadily and approach the upper bounds gradually with an increase of w. Particularly, c_nd and n_nd share the same growing trend with the increase of w that their $R_{\textrm{max}}\lt0.55$ (figures 8(a) and (b)). Also, c_nd is insignificantly larger than n_nd with the same w and the same group n = 30 or n = 60. Comparatively, s_rd and n_rd perform better than c_nd and n_nd regarding the correlation coefficients R. Meanwhile, s_rd overtakes n_rd notably for both n = 30 and n = 60, and s_rd shows the most prominent correlation with $R_{\textrm{max}}\thickapprox 0.60$ when w = 3 km (figures 8(c) and (d). Two phenomena can be observed from the figure. On the one hand, s_rd obtains the strongest impact on Δd, which is convincing since s_rd incorporates the NO₂ source from roads and NO₂ dispersion following the Gaussian distribution into the correlation. For instance, a long road is likely to produce more NO₂ while a short distance to the station is related to a greater amount of NO₂. On the other hand, w = 2.25 km is an empirical distance that road networks start to have a rather weak impact on Δd.

Zoom In Zoom Out Reset image size

The present study considers that regional climates can also affect the dispersion of NO₂. It is based on the evidence that air temperatures were used to describe urban heat islands and climate largely determines the magnitude of urban heat islands (Zhao et al 2014, Manoli et al 2019). As both air temperatures and NO₂ are air properties essentially, NO₂ concentrations may also be influenced by climate. To explore the combined impacts of urban mobility and road networks in partition of different climates, the correlation analysis is performed in different climates, namely, arid and humid climates, hot and cold climates, and coastal and inland climates (figure 9(a)). The correlation analysis is performed between Δd and $\Delta m_{rd} = s_{rd}\cdot \Delta m_w$, where s_rd follows the Gaussian function with w = 3 km that obtains the largest R and Δm_w is the change of workplace mobility as discussed in figure 6. Generally, results obtained from n = 30 (figures 9(b)–(d)) and n = 60 (figures 9(e)–(g)) are almost the same that (Δd, Δm_rd ) are strongly correlated. Particularly, the correlations are more prominent in the arid, hot, and coastal climates with $R = \{-0.84,-0.65,-0.77\}$ than in the humid, cold, and inland climates, namely $R = \{-0.60,-0.64,-0.64\}$ (figures 9(b)–(d)). In contrast, correlations based on n = 60 are slightly weaker than n = 30. That is to say, the arid, hot, and coastal climates tend to facilitate the dispersion of NO₂ based on a given urban mobility and road networks. It is also found that Δm_rd shows the most prominent correlation with Δd, and the significance of the correlations decreases from Δm_rd with the Gaussian distribution, s_rd with the Gaussian distribution (figure 8(d)), s_rd (figure 7(c)), to Δm_w (figure 6) when they are in the same condition. It suggests that the dynamic mobility constraint by static road networks significantly affects the changes of NO₂ in a large geographical space.

Zoom In Zoom Out Reset image size