## 2 Materials and Methods
### Spatial and Temporal Analysis  
Thermal analysis was carried out using the Landsat 8 satellite sensor image, and the land surface temperature (LST) was calculated from band 10 and band 11 by using land  
Figure 3: Land-use map of the study area and location of air sampling.  
surface emissivity (\(\varepsilon\)). In addition, the normalized difference vegetation index (NDVI) was calculated using band 4 and band 5 to determine the proportion of vegetation (\(PV\) and \(\varepsilon\)) to estimate the LST [32; 33; 34; 35], as shown in Figure 4. Three images from May 2014, 2017, and 2020, with a ground resolution of 30 m, were selected.  
At the same time, the study involved the measurement of particulate matter (PM10, PM4, PM2.5, and PM1) and gases (CO, CO\({}_{2}\), NO, NO\({}_{2}\), SO\({}_{2}\), O\({}_{3}\), H\({}_{2}\)S, and Volatile organic compounds VOCs) using portable instruments at 18 points from 18 to 25 May 2014. The measuring was according to Presidency of Meteorology and Environment Protection standards (PME) in Riyadh and performed by Geotechnical and Environmental Company (GECO), as shown in Table S2. The distance between sample sites was 2-3 km, as shown in Figure 3.  
The study used descriptive analysis for air sampling, to picture and summarize the data [36], and used correlation coefficients to assess the strength of linear relationship between variables. Principal component analysis (PCA) was carried out to transform the original set of variables into a smaller set of linear combinations [37].  
Data analysis was carried out geostatistically, using ArcGIS 10.1, while ordinary kriging was done by creating a semivariogram. The study estimated the central tendency, dispersion, and shape to determine the type of distribution of numerical variables. A semivariogram was created by building the spatial variability structure of each attribute, as in shown in Equation (1).  
\[\gamma(h)=\frac{1}{2N(h)}\sum_{i=1}^{N(h)}\left[\left(Z(x)-Z(x+h)\right)^{2}\right. \tag{1}\]  
where \(\gamma(h)\) is the semi-variance, \(Z(x)\) is the value of initial potential at site \(x\), \(Z(x+h)\) is the value of potential at site (\(h\)) distance apart from (\(x\)), and \(N\) is the number of sample pairs [38].  
Ordinary kriging maps were established to predict values at a non-sampled point that assumes a constant unknown meaning, Equations (2a) and (2b):  
\[Z(X_{0})=\mu+\varepsilon(X_{0}) \tag{2a}\]  
Figure 4: Flowchart to determine land surface temperature.  
\[Z(X_{0})=\sum\lambda_{i}\,\gamma(x_{i})\,\ \ \ \sum\lambda_{i}\ =\ 1 \tag{2b}\]  
where \(\mu\) is an unknown constant, \(\varepsilon(X_{0})\) is the the error associated with an unknown location \(X_{0}\), \(Z(X_{0})\) is the estimated value of \(Z\) at \(X_{0}\), and \(\lambda_{i}\) is the the weight that gives the best-possible estimation from the surrounding points.  
Finally, the trend of the monthly average of NO\({}_{2}\) and PM10 concentrations was determined for the Southern Ring Road air station from January 2016 to September 2018 and the Al-Khaldiyah air station from September 2019 to May 2020. These two air stations were selected as they are the closest to the study area and have air quality conditions similar to those in the study area. Moreover, a continuous data set for the time period mentioned above was available from these two air-quality-monitoring stations.