Supplementary Information for Improving air quality assessment using physics-inspired deep graph learning

Existing methods for fine-scale air quality assessment have significant gaps in their reliability. Purely data-driven methods lack any physically-based mechanisms to simulate the interactive process of air pollution, potentially leading to physically inconsistent or implausible results. Here, we report a hybrid multilevel graph neural network that encodes fluid physics to capture spatial and temporal dynamic characteristics of air pollutants. On a multi-air pollutant test in China, our method consistently improved extrapolation accuracy by an average of 11–22% compared to several baseline machine learning methods, and generated physically consistent spatiotemporal trends of air pollutants at fine spatial and temporal scales.

where C is the observed concentration of air pollutant, x is the input vector, N is the number of samples for 1 e , M is the number of samples (including N supervised collocated samples) for 2 e , i d is the temporal index of day for sample i, Since the dependent variable (C) is log-transformed into log( ) C C   and so the continuity equation (Eq. 1) changes correspondingly and thus 2 e can be approximated using 2 e .We minimize the loss function, so 2 e will be close to 0 and we can make the predictions to maintain the continuity equation as much as possible.
Here, the time variable of day (d) and the horizontal directions (x and y) were used in the residual term, 2 e .The extraction of 2 e is shown in Algorithm 2.

Supplementary Note 3: Neighborhood heatmap and SHAP values
The output from the final graph convolution layer represents the spatial agglomeration of neighborhoods in the study area, and the linear activation was used for the output of the last graph convolution.For a daily grid of 1x1 km 2 over the study area, the outputs of the final graph convolution layers for all the pixels that represent the spread feature of air pollution at the longest spatial scale were merged into a heatmap, which provides a visualization of the contours of the neighborhood.As a game-theoretic approach, SHAP (Shapley Additive exPlanations) 3 measures the contribution of each feature to the prediction taking into account the interaction with other variables.It is calculated by comparing what a model predicts across all possible combinations of the other covariates with and without the particular feature under analysis 4 .With the national mean concentration of an air pollutant as the background value and 1000 selected observations as background samples, the SHAP value was estimated for each covariate of every prediction.Then the statistics of SHAP values were reported to evaluate spatiotemporal distribution of influential factors.We estimated the contribution of each covariate to PM2.5, especially those related to emissions.

Supplementary Note 4: AQI output
The air quality sub-index for each pollutant was calculated 5 : Then, the AQI was calculated: where IAQI is a sub-index of AQI for each air pollutant, n is the number of air pollutants.The primary air pollutant is defined as having the maximum sub-index of AQI.
In addition, we calculated the national averages of air pollutant concentrations and AQI weighted by the 1x1 km 2 population density.The population-weighted concentration or AQI was defined as: where ˆp M is the population-weighted air pollutant concentration, AQI sub-index or AQI, i is the cell index, i p  is the estimated population density for cell i, and i c  is the estimated pollutant concentration, AQI sub-index or AQI for cell i.We got the population density data from https://www.worldpop.org.Note: The variables including the MERRA2-GMI input of pollutants (shown in gray background) were not used in our sensitivity analysis.

Supplementary Table 2
Hyperparameters tested in the sensitivity analysis for the compared machine learning or statistical regression methods.

Algorithm
y coordinate for sample i,   , f  W b x is the output of the model, RSE  is the mean square error (MSE) loss function, ,  w b is the set of parameters (weight matrix and bias) and   ,   W b is the normalization of ,  w b .

i
is the sub-index of air quality index for the air pollutant, p, Cp is the measured or estimated concentration for, p, are respectively the low and high values closest to the limit of Cp, H IAQI i is the sub-index of AQI corresponding to H BP i , and o L IAQI is the sub-index of AQI corresponding to L BP o .
Figure 3 Site-based test RMSE boxplot of seven air pollutants by eight methods.a, c, e, g, i, k and m show testing RMSE and b, d, f, h, j, l and n show site-based independent testing RMSE.Method: Deep graph modeling (DGM) proposed; Full residual encoder-decoder (FRE); GraphSAGE (SAGE); Random forest (RF); XGBoost (XGB); Kriging (Kri); Regression kriging (RKri); Generalized additive model (GAM).The whiskers (solid horizontal lines) represent the minimum and maximum RMSE, excluding outliers.The box bounds the interquartile range (IQR) from the 25 th to 75 th percentiles, with the 50 th percentile (median) as the bold line inside.The whiskers (dashed vertical lines) extend 1.5 times the IQR from the box, indicating the normal RMSE range.Data points (small circles) outside the whiskers are RMSE outliers.Supplementary Figure 5 Comparison of RMSE violin plots between our DGM method and various machine learning approaches, illustrating the predicted versus observed time series for all sites in the site-based independent test.) of the emission-related components to the PM2.5 predictions for four representative regions of China.a, China Mainland.b, the Jing-Jin-Ji metropolitan region.c, Yangtze River Delta.d, Pearl River Delta.e, Xinjiang Zhuang Autonomous Region.BC: Black carbon mass mixing ratio; DMSSMASS: DMS surface mass concentration; DUST: Dust surface mass concentration; GAE: Ground aerosol extinction coefficient; HNO3SMASS: Nitric acid surface mass concentration; Neighborhood: Neighborhood feature extracted using GC; NI: Nitrate Mass Mixing Ratio; NO: Nitric oxide; NO2: Nitrogen dioxide; O3: Ozone; OCS: Organic Carbon Surface Mass Concentration; SO2: SO2 Surface Mass Concentration.Supplementary Figure 7 Spatiotemporal evolution of PM10 for the sandstorm event in North China of April 15, 2015.a-f, national (the upper part) and local (the lower part of the Jingjinji metropolitan region) spatial distributions of estimated PM10 concentrations from April 11 to 16, 2015.The blue arrow is mean ground wind vector (m/s), the color gradient lines in each upper image show the contours of the graphical convolution heatmap, and the gray lines in each below enlarged image show major traffic roads.Supplementary Figure 8 Spatiotemporal evolution of PM2.5 for the sandstorm event in North China of April 15, 2015.a-f, national (the upper part) and local (the lower part of the Jingjinji metropolitan region) spatial distributions of estimated PM2.5 concentrations from April 11 to 16, 2015.The blue arrow is mean ground wind vector (m/s), the color gradient lines in each upper image show the contours of the graphical convolution heatmap, and the gray lines in each below enlarged image show major traffic roads.Supplementary Figure 12 Spatiotemporal evolution of NO2 for the haze event in Shanghai of December 19, 2018.a-f, national (the upper part) and local (the lower part of the Shanghai) spatial distributions of estimated NO2 concentrations from December 15 to 20, 2018.The blue arrow is mean ground wind vector (m/s), the color gradient lines in each upper image show the contours of the graphical convolution heatmap, and the gray lines in each below enlarged image show major traffic roads.

2 of different methods across seven air pollutants.
RFRandom forest: a method of averaging multiple deep decision trees, trained on different parts of the same training set, with the goal of reducing the variance.Number of trees (50-500; 50), criterion (mse), min samples split (1-10; 2), bootstrap (True, False), max feature ('auto') XGB XGBoost: an optimized distributed gradient boosting library designed to be highly efficient, flexible and portable.It is run under the gradient boosting framework.aOrdinary kriging: estimates a value at a point of a region for which a variogram is known, using data in the neighborhood of the estimation location.With local secondorder stationarity, it implicitly evaluates the mean in a moving neighborhood.aNo covariates used for ordinary kriging; all the same covariates used for all the other methods.

Table 7 Comparison of RMSE for the temporal trends: summarized monthly means of air pollutants for all sites in the site-based independent test.
Cor.: Pearson correlation with the ground observed values; b R 2 between estimates and observed values. a