## 2 Materials and Methods
### Long-Term Trends and Meteorological Adjustment  
In the following, trends in dry-season average visibility at all seven TMD stations (CB, LC, RY, ST, CT, SK, and PB) were estimated using linear regression. A non-parametric Mann-Kendall test was employed in conjunction with block bootstrapping (to account for serial correlation in a time series) to test the significance of trend (against the null hypothesis of no trend) [58]. At CB, ST, CT, and PB, long-term data (since 1981) are available; it was then possible to perform trend estimation over two consecutive periods: P1: pre-2000 (1981-2000); and P2: post-2001 (2001-2015). For LC, SK, and RY, the trend analysis was limited to the recent period (starting with the years 2001, 2001, and 2007, respectively) only.  
In order to account for the effect of meteorology on trends in visibility, a generalized linear model (GLM) framework, used by [PERSON] et al. (2007) [59], was adopted. GLM is an advanced statistical model capable of accounting for relationship between response (i.e., predictand) and explanatory (i.e., predictor) variables via a link function [60; 61]. Visibility at each of the CB, LC, and RY stations is the interested response variable, while meteorological variables are the predictors, whose GLM formulation is given by:  
\[g(\mu_{i})\sim f_{1}(x_{i,1})+f_{2}(x_{i,2})+\ldots f_{j}\big{(}x_{i,j}\big{)} +\ldots+f_{p}\big{(}x_{i,p}\big{)}+Y_{k}, \tag{4}\]  
where subscript \(i\) denotes the \(i\)th observation and subscript \(j\) (1, 2, 3, \(\ldots\), \(p\)) denotes the \(j\)th predictor variable, \(f_{j}\big{(}x_{i,j}\big{)}\) is the smoothing function (here cubic spline) of the \(j\)th predictor at the \(i\)th observation, \(\mu_{i}\) is the \(i\)th value of the response variable, \(Y_{k}\) is the effect of the \(k\)th calendar year, and \(g(\mu_{i})\) is the link function relating the mean of the response variable to the predictors. Here, GLM was implemented in the R software [56]. Several options of error distribution and link function were tried, and Gaussian distribution and identity link were finally chosen since they tend to give an optimal fit. The adjusted coefficient of determination (shortly, adjusted \(R^{2}\)), representing variance of data explained (%) by the model, was used as the goodness-of-fit measure for each GLM model. Residual diagnostics (here, independent error, normality of error, homoscedasticity of error, and outlier leverage) were inspected and found to be fairly acceptable [55]. A parsimonious approach was adopted to achieve the final model. Eight predictors were initially chosen for the GLM, namely temperature, relative humidity, cloud cover, mixing height, wind speed, recirculation factor, back-trajectory cluster, and persistence (Supplementary Materials, Table S2). Persistence was here assigned as the daily value from the previous day and included to represent inter-day impact or imprint generally due to atmospheric stagnation or lack of atmospheric ventilation. It is noted that certain variables were not included in this initial list due to their correlation with one of the initial variables, e.g., global radiation and cloud cover. Similarly, the entire GLM procedure was applied to PM\({}_{10}\) but with log transformation. Given Equation (2), relative humidity with and without log transformation was tried but gave comparable results in both visibility and PM\({}_{10}\) cases. Thus, no log transformation was used here. In the GLM processing, stepwise backward elimination was carried out to obtain a parsimonious fit. Given the final GLM fit achieved, a meteorologically adjusted trend (adjusted trend) was computed from the partial response with respect to the effect of the calendar year (\(Y_{k}\)).