Impact of weather changes on air quality and related mortality in Spain over a 25 year period [1993–2017]

Climate change is a major public health concern. In addition to its direct impacts on temperature patterns and extreme weather events, climate change affects public health indirectly through its influence on air quality. Pollution trends are not only affected by emissions changes but also by weather changes. In this paper we analyze air quality trends in Spain of important air pollutants (C6H6, CO, NO2, NOx, O3, PM10, PM2.5, and SO2) recorded during the last 25 years, from 1993 to 2017. We found substantial reductions in ambient concentration levels for all the pollutants studied except for O3. To assess the influence of recent weather changes on air quality trends we applied generalized additive models (GAMs) using nonparametric smoothing; with and without adjusting for weather parameters including temperature, wind speed, humidity and precipitation frequency. The difference of annual slopes estimated by the models without and with adjusting for these meteorological variables represents the impact of weather changes on pollutant trends, i.e. the ‘weather penalty’. The analyses were seasonally and geographically stratified to account for temporal and regional differences across Spain. The results were meta-analyzed to estimate weather penalties on ambient concentration trends at a national level as well as the impact on mortality for the most relevant pollutants. We found significant penalties for most pollutants, implying that air quality would have improved even more during our study period if weather conditions had remained constant. The largest weather influences were found for PM10, with seasonal penalties up to 22 μg⋅m accumulated over the 25-year period in some regions. The national meta-analysis shows penalties of 0.060 μg⋅m per year (95% Confidence Interval, CI: 0.004, 0.116) in cold months and 0.127 μg⋅m per year (95% CI: 0.089, 0.164) in warm months. Penalties of this magnitude would correspond to 129 annual deaths (95% CI: 25, 233), i.e. approximately 3200 deaths over the 25-year period in Spain. According to our results, the health benefits of recent emission abatements for this pollutant in Spain would have been up to 10% greater if weather conditions had remained constant during the last 25 years.


Introduction
Climate change is widely recognized as a major public health threat (Watts et al., 2015) with a wide range of impacts including food and water security, flooding and sea level rise, spread of infectious diseases and extreme weather (WHO, 2018). Previous studies have reported increased mortality directly related to the modification of temperature patterns induced by climate change (Ren et al., 2011;Vardoulakis et al., 2014;Carmona et al., 2016;Lee et al., 2018). In addition, climate change influences air quality (Jacob and Winner, 2009) as ambient air pollutants are very sensitive to meteorological conditions (Elminir, 2005;de la Paz et al., 2016;Westervelt et al., 2016;Chen et al., 2018a) and exposure to airborne pollutants is also a leading contributor to global disease burden (Lelieveld et al., 2015;Cohen et al., 2017). Recent studies have suggested that two-way interactions between weather variables (e.g., temperature) and air pollution should be carefully considered to characterize synergistic effects on health (Stafoggia et al., 2008;Chen et al., 2018b), especially in the context of climate change Silva et al., 2017). Usually risk estimates from epidemiological studies model temperature as a confounder while few studies has investigated the role of temperature as an effect modifier to short-term exposure to pollutants such as ozone, whose dynamics https://doi.org/10.1016/j.envint.2019.105272 Received 22 January 2019; Received in revised form 2 October 2019; Accepted 14 October 2019 strongly depends on temperature too. Jhun et al. (2014) concluded that heat exposure might exacerbate adverse health impacts of ozone, but also found that interactions between air pollution and temperature are non-linear and very complex to translate into health effects, especially long-term effects. Nonetheless, since climate and air quality are closely intertwined, they should be addressed in an integrated manner in environmental policies. Previous research indicates that adverse health outcomes of climate change-induced air quality changes can be offset by future emission reductions (Stowell et al., 2017) and that benefits of these climate mitigations exceed their costs (Balbus et al., 2014;West et al., 2012;Xie et al., 2018). Up-to-date information on the extent of such effects is essential to the design of effective measures to abate harmful atmospheric emissions to protect public health.
In this paper, we investigate the influence of changes of meteorological conditions over the last 25 years on air quality trends in Spain. From here on, we will refer to the trends of meteorological conditions in this period as 'weather changes'. Querol et al. (2014) studied air quality trends in Spain for the period 2001-2012 and suggested that while air pollution reductions have been primarily driven by emission reductions, meteorological changes may also have influenced air quality trends throughout that period. Other time series studies on air pollution have also considered weather changes as an influential factor in Europe (Barmpadimos et al., 2012;Cusack et al., 2012;Garrido-Perez et al., 2018). While many regulated pollutants (SO 2 , NO x , CO, PM 10 and PM 2.5 ) have declined substantially over recent years, O 3 levels have remained constant or increased in some cases. The literature has proposed that O 3 has been difficult to control due to NO x reduction policies (Henschel et al., 2015;Jhun et al., 2015b) that may have induced unanticipated changes on the oxidant capacity of the atmosphere and the production of O 3 (Saiz-Lopez et al., 2017).
However, it has been demonstrated that tropospheric O 3 is very sensitive to meteorological influences both at synoptic (Santurtún et al., 2015;Ramos et al., 2018) and regional (Reche et al., 2018) scales, so the attributions to weather changes remain uncertain. Recent studies in the United States suggest that air quality would have improved more in recent years if weather conditions had remained constant, thus introducing the concept of 'weather penalty' (Jhun et al., 2015a). We apply similar methods to understand the influences of weather changes on air quality trends in Spain and the health impacts attributable to them. Of note, we do not intend to provide a holistic assessment of the impact of recent weather changes. Our goal is to isolate the effect that those changes have had on ambient pollution levels, and then the direct health implications of such influences only considering long-term effects of individual pollutants in all-cause mortality. In our analysis, we assess the impact on mortality of the pollutants most relevant to public health, namely particulate matter (Pope and Dockery, 2006;Kim et al., 2015;Burnett et al., 2018), O 3 (Malley et al., 2017;Díaz et al., 2018) and NO 2 (Samoli et al., 2006;Mills et al., 2015). However, we also estimate weather penalties for other major pollutants to gain a better understanding of the relationships between air pollution and weather variables.

Observational datasets
Air quality data were obtained from the Air Quality Area of the Spanish Ministry for the Ecological Transition. They provided a database of all the observations used for official air quality assessment in Spain during the last 25 years, from 1993 to 2017. This dataset contains hourly and/or daily (for some PM monitoring sites) concentrations of C 6 H 6 , CO, NO 2 , NO x , O 3 , PM 10 , PM 2.5 , and SO 2 from a total of 887 monitoring sites.
We selected monitoring sites that had at least 10 years of observations with a minimum validated data availability of 90% annually (following the European Union Air Quality Directive quality criterion of data capture for assessment purposes for fixed measurements). The daily average and maximum 1-h concentration (also the daily maximum 8 h-running average for O 3 ) were computed, excluding days with less than 18 valid hourly observations. That yielded a dataset of 7,266,609 daily records from 354 air quality monitoring sites. We present the spatial distribution of the air quality sites in Fig. 1 and a summary of each pollutant, including sample size in Table 1. Sites were  categorized into seventeen regions, corresponding to the Spanish administrative division (European Nomenclature of Territorial Units for Statistics; NUTS-2 level), as listed in Fig. 1. It should be noted that, the Canary Islands (around 1700 km south-west away from the center of the Iberian Peninsula) are also included in our analysis. Hourly meteorological data on temperature (°C), wind speed (m/s), water vapor pressure (hPa) and precipitation (mm) were derived from the observations of the Spanish Meteorological Agency (AEMET) weather monitoring network. The selection of these weather parameters was done accordingly to the outcomes of previous studies on the sensitivity of air pollution to meteorological factors (Jhun et al., 2015a and references therein). We computed daily averages to match the metrics used for air pollutants and excluded days with less than 18 valid hourly records. In the case of precipitation, accumulated daily rain was computed and a 0/1 value was assigned for each day since precipitation frequency is a more relevant metric than rainfall intensity for wet deposition processes (Jacob and Winner, 2009), especially for particulate matter. Similarly to that of air pollution, we carried out an analysis for the selection of the monitoring based on data availability. In this case, we selected only sites where a minimum of 90% of all the four variables were simultaneously available. This criteria produced a dataset of 348,288 daily records from 47 weather monitoring sites ( Fig. 1) with an average length of the available data of 20.3 years. Weather and air pollution monitoring sites were matched according to a minimum distance criteria, always within the same region. As a result, the median distance between air quality monitoring sites and weather stations was 26.5 km for O 3 and 24.2 km for PM 10 and 27.3 km for NO 2 .

Trends and weather penalty calculation
We applied generalized additive models (GAMs) using nonparametric smoothing (Hastie and Tibshirani, 1990) to explore the relationship among ambient concentration levels and relevant weather parameters in Spain. This statistical approach has been used in several studies to (i) adjust for inter-annual meteorological variation using smoothing spline functions (Cox and Chu, 1996;Pearce et al., 2011;Barmpadimos et al., 2012); (ii) isolate pollution effects (Wilson et al., 2004;Pope and Dockery, 2006); or (iii) assess the efficacy of emission mitigation efforts (Zheng et al., 2007;Jhun et al., 2013). GAMs are a powerful tool due to their flexibility to fit different time scale trends (seasonal, long-term, etc.) as well as nonlinear associations among variables, making them ideal for our research purpose. We have followed the methodology used in Jhun et al. (2015a), that consists in comparing two closely related GAMs fit to the same long-term pollution trend; one without adjusting for weather parameters (Eq. (1)) and other adjusting for weather parameters (Eq. (2)): represent the linear pollutant annual trend (weather unadjusted and adjusted respectively), expressed in μg⋅m −3 per year; γ and δ are vectors of coefficients that explain monthly and weekday variability within the time series, respectively; and s s s s , , , 1 2 3 4 denote the smoothing spline functions that take into account the nonlinear relationships between the air quality index scrutinized and temperature (°C), wind speed (m/s), water vapor pressure (hPa) and precipitation (frequency) respectively in the weather adjusted model. The weather adjusted model (Eq. (2)) removes the inter-annual influence of meteorological changes, i.e., β adjusted 1, represents the trends that pollution would have followed if weather parameters would have remained constant during the 25 year period analyzed. Therefore, any difference between β unadjusted 1, (that represents the actual trend, including the influence from weather) and β adjusted 1, (weather-adjusted trend) denotes the change in air quality attributable to changes in meteorological conditions or 'weather penalty' (μg⋅m −3 per year) (Eq. (3)).
A positive penalty ( > β β unadjusted adjusted 1, 1, ) implies that weather changes led to poorer air quality (either greater increases or smaller reductions in ambient concentration). Negative values indicate that changes in meteorology dampened trends in poor air quality, implying a weaker increasing trend or stronger decreasing trend. All the GAMs discussed in this paper were computed through the gam function in mgcv package (Wood, 2011) in R (R Core Team, 2018).
The procedure described above provides the central estimate of weather penalties. In order to provide a measure of the uncertainty, we also computed the standard errors of both trends and the penalties. We applied a block bootstrap procedure that consists of creating randomized subsets of the actual pollution data, or pseudo-datasets (Politis, 2003). Following the methodology of Jhun et al. (2015a), a 20-day block size was considered to generate 100 pseudo-datasets for each region, season and pollutant. The corresponding adjusted and unadjusted betas and thus, the weather penalties, were computed for each of the 100 pseudo-datasets and the standard deviation of these 100 realizations was used to derive standard errors (i.e., the confidence intervals for the central estimates), that accounts for site-to-site heterogeneity within each region. Finally, region-specific betas and their differences (penalties) were meta-analyzed through a random-effects model (Berkey et al., 1998) to combine within-region and betweenregion variability into a national level estimate in order to provide an aggregated assessment of weather-related changes in pollution in Spain during the 1993-2017 period.
In addition to GAMs for pollution series, a general linear regression model (Chambers, 1992) was applied to each meteorological time series to estimate the trends of temperature, wind speed, water vapor pressure and precipitation frequency (binomial regression) to help identifying the causal relationships between the weather penalties obtained and the trends of the main meteorological variables. Similarly, regional trends were meta-analyzed to obtain aggregated national trends, consistent with the pollution trends results.
It should be noted that all the analyses were seasonally and geographically stratified to account for temporal and geographical differences in our study. This is relevant considering the spatio-temporal variability of the weather parameters reported for Spain in the literature: temperature ( (Gallego et al., 2011;Herrero et al., 2011;Cortesi et al., 2012). For consistency with previous studies (Jhun et al., 2015b), we stratified our analyses by warm (May-October) and cold (November-April) season, although the monthly weather variability is also discussed in Section 3.1. All trends, penalties and health effects were specifically assessed for the seventeen regions NUTS-2 regions shown in Fig. 1.

Health impact estimation
Our interest is in the mortality impact of the differences in air pollution that occurred due to these weather changes, both as absolute values and in comparison with the mortality reductions that have occurred due to improvements in air pollution since 1993.
If the mortality observed in year y (M y ) is assumed to reflect the pollution experienced in that year for a given pollutant, C y , under a relative risk model , where M 0 is the mortality that would be expected if pollution levels were below threshold levels and RR C ( ) y is the relative risk corresponding to C y . Although it may well be true that the effects of exposure to air pollution may be experienced several years after exposure, it is not uncommon for air pollution risk assessments to make this simplifying assumption. If it is further assumed that the concentration-response function relating air pollution C to mortality risk is essentially linear, with a slope β (fractional increase in RR per µg/m 3 ), then = − RR e β C T ·( ) for > C T and 1 if ≤ C T, where T is the threshold. Under these circumstances, the mortality observed in year y can be expressed as: If the pollution level in year y had been different, say C y a , , then under this model the observed level of mortality would also have been different: The difference between the mortality observed at C y , and under the alternative conditions, C y a , , is the quantity of interest in this study. As long as both concentrations are above the threshold, then: ) . Once this is substituted into Eq. (6) we obtain: Which can be simplified to yield: Using this expression, the difference between the number of deaths attributable to air pollution under any two scenarios of interest ( M Δ ) can readily be computed by Eq. (9): where the difference of concentration between these two scenarios ( − C C y ya , ) is noted as δ. The estimate is valid under the assumptions that -(i) a relative risk model holds; (ii) either (a) the relationship between mortality risk and pollution exposure is linear, without a threshold; or (b) both C y and C y a , are above any threshold; and (iii) the full impact of air pollution exposure is experienced in the year of exposure (i.e., no lag).
This approach, which is consistent with that of the US Environmental Protection Agency's BENMAP tool (US EPA, 2018) and the World Health Organization (WHO, 2018), has been used to compute -(i) the change in mortality attributed to air pollution due to the weather penalties derived above ( = δ Penalty), and (ii) the change in mortality attributed to air pollution due to the improvements in air quality since 1993 ( = δ β unadjusted 1, ). Although our study includes other pollutants, we limited our health impact analyses to PM 10 , O 3 and NO 2 . There is a breath of β coefficients (fractional increase in RR) in the literature for the three pollutants considered. Although some studies on short-term effects have been published for Spain , we are mainly interested on long-term effects of air pollution and no country-specific cohort studies are available for this country. In particular, we focus on mortality data for all causes, since they tend to be more reliable than cause-specific mortalities (WHO, 2013a) and provide a better overview of the overall impact of pollution trends or the weather penalties for this particular study.
For consistency, we considered the RR proposed by the multicenter European Study of Cohorts for Air Pollution Effects (ESCAPE) (Beelen et al., 2014) for PM 10 and NO 2 . This study investigated the association between natural-cause mortality (A00-R99 according to the International Classification of diseases-10th revision) and long-term exposure to several air pollutants using data from 22 European cohort studies, which created a total study population of 367 251 participants. They reported a 4.0% increase (95% CI: 0.0%, 9.0%) in mortality risk for an increase of 10 μg⋅m −3 in average PM 10 levels, i.e. a RR of 1.040 (1.000, 1.090). The corresponding RR increase for an increase of 10 μg⋅m −3 in average NO 2 is 1.010 (0.990, 1.030) according to this study. Although it is very hard to strictly allocate observed effects to single pollutants (Künzli et al., 2000), many studies have proposed lower RR for NO 2 , in comparison with those of PM 10 . Nonetheless, NO 2 has been positively associated with all-cause mortality increases in European cities (Carey et al., 2013;Giulia et al., 2013). In addition, the reduction of NO 2 concentrations has been the primary target of urban air quality plans in cities such as Madrid (Borge et al., 2014;Borge et al., 2018). It should be noted that in our analysis we limited NO 2 RR to values ≥1 that are more coherent and plausible from the health perspective and prevent misinterpretations of the resulting trends and weather penalties. For ozone, the American Cancer Society Cancer Prevention Study (CPS) (Jerrett et al., 2009) reported a 4.0% increase (95% CI: 1.0%, 6.7%) in the risk of death from respiratory causes per 10 ppb increase in daily maximum O 3 concentrations (in the April -September period) among people age 30 and older in the US, although no significant association with all-cause mortality was found (RR of 1.001 (0.996, 1.007) in that study. However, this estimate was recently updated under the CPS-II epidemiological cohort study (Turner et al., 2016) which reported a mortality risk increment for the same population range (all causes, including external) associated with long-term O 3 exposure of 2% (95% confidence interval: 1-4%) for a 10 ppb increase in the mean annual daily 8-hour maximum. Since this is a more recent reference and it may provide a better view of potential O 3 across different seasons (Malley et al., 2017), we have considered this RR in our assessment.
Previous studies (Cooke et al., 2007) support that mortality from long-term exposure risks from global studies (Cohen et al., 2005) may be extrapolated to populations in other in other geographic regions. In fact, the RR proposed by Jerrett et al. (2009) for respiratory causes are those recommended by the Health risks of air pollution in Europe (HRAPIE) project (Héroux et al., 2015) experts to assess long-term mortality due to tropospheric O 3 in Europe (WHO, 2013b).
We used official data from the Spanish Statistical Office (INE, 2018) to incorporate mortality data, (M y a , ) integrated in 6 month-periods for specific age ranges (the figures are provided in the supplementary material; Table S1 and Tables S2a and S2b). This allowed us to make use of season-specific trends and apply Eq. (9) by season. This approach complies with the ethical and legal requirements according to the data protection law. Although ESCAPE cohorts differed in the mean baseline age, we consistently considered the population age 30 and older for the health impact assessment.
To provide some characterization of the uncertainty in our estimates of the mortality impacts of these weather changes we apply Gauss' Law of Error propagation (known by statisticians as the delta method) as was done by Agresti (2012). Gauss law suggests that the error in a function may be estimated as the sum of products of the partial derivative of the function with respect to each uncertain variable times an estimate of the variance of each variable. If the uncertain variables are correlated, then cross-products reflecting their covariance must be added.
Using this approach we combined the standard errors of β (from the confidence intervals of the original concentration-response function reference) and δ (obtained from the bootstrapping procedure discussed in the previous section) to estimate the total uncertainty in estimates of mortality impacts. Assuming that both variables are unrelated and thus, the error covariance is zero (according to Eq. (10)).
This approach provides a valid estimate of the impact of parameter uncertainty as long as the errors in both beta and delta are small. It does not address epistemic or model uncertainty, which in some cases may be much larger than parameter uncertainty.

Meteorological trends
The trends of weather parameters included in this study exhibited a substantial spatio-temporal variability that may have influenced air pollution trends, as discussed below. 3.1.1. Temperature Temperature changes during the cold season ranged from −0.082°C per year in the Basque Country (PV) to 0.050°C per year in Aragon (AR) (Fig. 2a). Regional trends are much more consistent during the warm period. Except for PV, temperatures increased all over the country, reaching values up to 0.083°C per year in Madrid (MD). That indicates that average temperature during the warmer months may have increased as much as 2°C in the center of the Iberian Peninsula in 25 years.
A national meta-analysis of regional temperature trends yielded a statistically significant variation of −0.07°C⋅yr −1 in the cold months and 0.040°C⋅yr −1 in the warm months (i.e. 0.4°C per decade). That represents an accumulated change of −1.7% and 4.9% relative to the respective national average values in the period of interest. These results ( Fig. 2b) are consistent with those from previous studies that have reported annual temperature increases around 0.1-0.2°C per decade with greater increases in the spring and summer (del Río et al., 2011).

Wind speed
According to our weather dataset, wind speed trends are rather variable across Spain during our study period. Most regions show a reduction in wind speed throughout the year (as much as −0.08 m⋅s −1 per year in Cantabria, CA) although an intense increment is observed in the Canary Islands (IC) (0.06 m⋅s −1 per year) (Fig. 2c). Despite a lack of evident spatial patterns, aggregated nation-wide results (Fig. 2d) indicate a dominant negative trend, especially during winter months, in agreement with the findings of Azorin- Molina et al. (2016). A global variation on wind speed of −0.3 m⋅s −1 over the 25-year period is observed, that corresponds to an 11.2% reduction relative to the average value, that is in line with the general trends in the northern hemisphere (Vautard et al., 2010).

Humidity
Water vapor pressure trends show a large variability, both spatially (Fig. 2e) and temporally (Fig. 2f). No statistically significant changes on water vapor pressure were observed at national level during the 1993-2017 period. This is consistent with previous findings suggesting that while relative humidity has significantly decreased in Spain over the last 50 years due to a rise in temperature, specific humidity had not changed significantly during that period . However, water vapor pressure may have decreased substantially in some inland areas such as Castille-La Mancha (CL) or Madrid (MD), especially during the warm season, with total variations in the annual mean pressure water pressure of −1.6 and −0.8 hPa over the 25 year period, respectively. On the other hand, a substantial increment is found in other regions such as the Balearic Islands (IB) (1.7 hPa). These results suggest that average annual water vapor pressure may have changed from −15.0 to 10.5% (relative to the respective regional mean values) across Spain during our study period.

Precipitation
Our meta-analysis of linear regression models shows a statistically significant increase in precipitation frequency during the cold season and a decrease in the warm season nationally (Fig. 2g). According to our results, precipitation frequency changed by 1.8 and −6.2% for the cold and warm season, respectively, over the period analyzed, with an average number of 59.8 and 41.5 rainy days (accumulated precipitation > 0) per year in the cold and warm season respectively. However, a monthly analysis at national level (Fig. 2h) does not show a consistent seasonal trend, due to a strong variability at regional level (Fig. 2g). The abundant literature on precipitation trends in the Spain (Vicente-Serrano et al., 2017) does not provide consistent conclusions due to differences on datasets, precipitation indexes used and methodological approaches. Gallego et al. (2011) and Herrero et al. (2011), i.a., found that both total precipitation amount and the number of rainy days exhibit large spatial variations as they are strongly affected by local geographical features. These studies have also identified an increasing variability of precipitation regimes in time as well as non-linear trends. This may explain the wide confidence intervals obtained in our metaanalysis (Fig. 2g).

Air quality trends and weather penalties
Unadjusted air pollution trends showed an overall decrease in pollution levels in Spain for both cold and warm seasons for all the pollutants analyzed except for O 3 . Although the methods and datasets differ, our results are in agreement with a recent analysis by Querol et al. (2014) for the 2001-2012 period. As illustrated in Fig. 3, NO 2 , PM 10 and PM 2.5 ambient levels have decreased remarkably during our study period. Concentration reductions are observed regardless of season but more prominently during the winter, with nationally-averaged accumulated changes up to −15 μg⋅m −3 for NO 2 (February) (Fig. 3a), more than −30 μg⋅m −3 for PM 10 (March) (Fig. 3c) and approximately −8 μg⋅m −3 for PM 2.5 (in February and March; in this case over a 17-year period) (Fig. 3d). In contrast, ambient O 3 levels (daily 8hour maximum O 3 moving average concentrations) have increased by nearly 10 μg⋅m −3 as an average. Increases have been higher during the cold seasons, up to 19 μg⋅m −3 in February (Fig. 3c). Increases for the summer months (June to August) were non-significant. In order to understand to what extent these air pollution trends were influenced by recent weather changes, we applied the statistical models discussed in Section 1. The unadjusted trends we estimated are the product of both emission and meteorological changes, while the weather-adjusted trends remove the influence of weather changes on air quality. Consequently, the difference between the unadjusted and weather-adjusted trends reflect the impact of long-term meteorological changes or weather penalties. Such penalties account for both direct effects (e.g., pollution transport and transformation phenomena) and indirect effects (e.g., changes in energy consumption) of weather changes. Meteorology affects a series of emission sources, both anthropogenic (e.g. those related to heating/cooling) and natural, such as mineral dust or biogenic VOC. One of the advantages of our methodology is that the influence of emissions is implicitly accounted for in the exact same way in both GAM models (unadjusted and weatheradjusted), so our penalties are not affected by that factor.
The nationally-aggregated weather-unadjusted and weather-adjusted trends as well as the corresponding weather penalty obtained from the meta-analysis of pollution trends at the regional level are summarized in Fig. 4 and Table 2. Of note, the effects of meteorological variables vary depending on the pollutant. In general, weather penalties for CO, a much less reactive specie, are smaller in relative terms (relative to the respective average concentration), suggesting that recent weather changes may have affected not only transport phenomena (e.g., advection), but also atmospheric chemistry or wet deposition processes. Below we report the results at the regional level for each of the pollutants considered in our health impact analysis (Fig. 4).
In Fig. 5 we compare unadjusted and weather-adjusted trends in a bubble plot. The magnitude of the penalty is given by the distance to the x = y line or identity line. Points close to that line imply that β unadjusted 1, and β adjusted 1, are similar in value and consequently, the weather penalty is smaller. In other words, it reflects a smaller influence of weather changes on recent air pollution trends in a given region. The size of the bubbles represent the uncertainty of the penalty estimate. A similar graph for the rest of pollutants can be found in Fig.  S1. Fig. 5a shows weather penalties for NO 2 by region. Consistent with Fig. 4d, the magnitude of the penalties (distance to the identity line) is larger for the warm period and usually positive > β β ( ) unadjusted adjusted 1, 1, , i.e., most of the points are in the lower-right half of the graph. However, the influence of weather changes on this specie differ significantly by region. Aragon (AR) had the highest penalty both in the cold and warm seasons (0.40 and 0.46 μg⋅m −3 ⋅year −1 respectively). This may be related to warmer and more stagnant conditions identified in this region over the study period. Temperature plays a fundamental role on NO 2 atmospheric chemistry (Atkinson, 2000). Higher temperatures would increase the oxidation rate of NO while weaker winds would limit dispersion. The latter may be more influential since the relative penalties are very similar to those found for NO x that can be collectively considered as an inert pollutant. Although CA (Cantabria) region presents a very similar variation on wind speed (Fig. 2c), NO 2 penalties are considerably lower (0.19 and 0.09 μg⋅m −3 ⋅year −1 for the cold and warm seasons, respectively). In addition, reduction of precipitation frequency in the summer would explain the higher relative penalty in Aragon (AR) during the warm season, given the sensitivity of NO 2 to wet deposition processes (Martin, 1984;Yoo et al., 2014). Observed changes on humidity may also play a role since it affects the relative abundance of the hydroxyl radical (OH), which in turn has a fundamental influence on the chemical processes that control the concentration of gaseous pollutants, including NO 2 (Beirle et al., 2011). This, in combination with moderate decreases or even increases in wind speed may be responsible for the positive influences of weather changes on NO 2 trends (negative penalty) observed in areas such as Castille-La Mancha (CM) or Madrid (MD). Given the interactions between weather parameters included in the weather-adjusted model, directly attributing changes in pollution trends to individual weather parameters is very challenging.

Ozone
Regional level penalties on O 3 are often the opposite of weather penalties on NO 2 in most regions (Fig. 5b), including Castille-La Mancha (CM) and Madrid (MD) regions, discussed above. Even though O 3 levels have increased more evidently during the cold season (Fig. 3b), weather penalties are usually negative for that period which is consistent with the result shown in Fig. 4g. In contrast, our results indicate that weather changes had favored higher O 3 concentrations in the warm season. The greatest regional penalties were found in Extremadura (EX) during the warm season (0.50 μg⋅m −3 ⋅year −1 ). Summers in EX have become substantially warmer, drier, and stagnant over our study period, which is consistent with previous studies that have looked into the sensitivities of O 3 to these weather parameters (Tu et al., 2007;Murphy et al., 2007;Kumar et al., 2015). However, a penalty of 0.18 μg⋅m −3 ⋅year −1 was identified for Madrid (MD) during the cold season, even when temperature has decreased and wind speed has increased during this period. The reason for the O 3 penalty in this case, is likely related to indirect effects on atmospheric chemistry. We found negative penalties in this region both for NO x and C 6 H 6 , one of the most abundant aromatic hydrocarbons in urban areas (Cocker III et al., 2001) and thus, a good proxy for VOCs (Wang et al., 2010). Such weather-induced changes may have modified the oxidizing capacity of the atmosphere, in addition to the changes driven by recent emission reductions (Saiz-Lopez et al., 2017), ultimately leading to increases of O 3 . Other factors, not explicitly included in our analysis, such as the recent increase of solar radiation in the Iberian Peninsula (Sanchez-Lorenzo et al., 2013) may have had an impact on atmospheric photochemistry and thus, the O 3 yield. Other regions with significant weather penalties during the warm season are Castille and León (CL) (0.33 μg⋅m −3 ⋅year −1 ) and Castille-La Mancha (CM) (0.32 μg⋅m −3 ⋅year −1 ), also located in the central area of the Iberian Peninsula where summers are generally becoming increasingly dry and warm according to our weather trend analysis.

Particulate matter
According to the results summarized in Table 2, weather changes in Spain had the greatest influence on PM 10 among all the pollutants we analyzed. This is also reflected in Fig. 3c, with most regions plotted in the right-lower half of the graph reflecting weather-related increases in PM 10 for both the cold and warm seasons. Consistently with the results for NO 2 , the Aragon (AR) region presents strong penalties in both seasons (0.88 and 0.59 μg⋅m −3 ⋅year −1 for the cold and warm seasons respectively). In other words, if weather had remained constant, average PM 10 ambient concentration levels in this region would have had an additional decrease of around 18 μg⋅m −3 over the 25 year period considered in this study (22 μg⋅m −3 as an average for the cold season). Despite generalized air quality improvements in all regions regarding this pollutant, this has strong policy implications since the effect of air quality plans may have been substantially affected by R. Borge, et al. Environment International 133 (2019) 105272 weather trends. That being said, it is worth noting that the effect of weather changes on PM 10 levels in Asturias (AS), the region with the largest improvements regarding this pollutant (a total reduction of 36 μg⋅m −3 in the PM 10 annual mean from 1993 to 2017), was negligible (penalties of −0.03 and 0.04 μg⋅m −3 ⋅year −1 for the cold and warm periods respectively), implying that the successful reductions of PM 10 levels was primarily related to emission abatement measures. In contrast, weaker winds and reduced precipitation frequency may be related to the penalties found for Cantabria (CA) in the cold season (0.16 μg⋅m −3 ⋅year −1 ) or Castille-La Mancha (CM) in the warm season (0.42 μg⋅m −3 ⋅year −1 ). Previous studies have identified stagnation as one of the main causes for pollution increases in Europe (Garrido-Perez et al., 2018). However, other factors not explicitly considered in this analysis, such as long-range transport patterns, may have a significant effect on observed trends (Cusack et al., 2012;Santurtún et al., 2015). Changes on such factors may be the cause for the weather penalties observed in the Canary Islands (IC), where despite a substantial increase in wind speed, a positive penalty was observed both seasons (0.06 and 0.08 μg⋅m −3 ⋅year −1 for the cold and warm periods, respectively). The results for PM 2.5 indicate that the influences on this pollutant are less pronounced and probably less affected by modification in deposition processes, and more affected by the weather impact on the precursors of secondary components, relatively more important for this fraction (Querol et al., 2004;Salvador et al., 2012). In agreement with the results for PM 10 , larger penalties are found in the warm season, with values up to 0.14 μg⋅m −3 ⋅year −1 in Galicia (GA) and 0.10 μg⋅m −3 ⋅year −1 in Madrid (MD).

Health impacts
The estimates of the national-level impacts of these weather penalties on mortality are summarized in Fig. 6 and Table 3.
Recent reductions in PM 10 levels have had a significant impact on mortality (Fig. 6). The number of deaths attributable to PM 10 has  ) and penalties at the national level meta-analysis. The figures in brackets represent the 95% confidence interval.  The results for O 3 indicate that recent increases of the annual average daily maximum 8-h O 3 concentration, reference metric used by Turner et al. (2016), were associated with 145 deaths (95% CI: 63, 227) per year. Nationally aggregated penalties according to our regional trend meta-analysis suggest that mortality linked to weather penalties for this pollutant is non-significant (−6, 95% CI: −24, 12 deaths per year). Nonetheless, it should be considered that this is the result of aggregating opposite seasonal trends (typically negative penalties in winter and positive in summer), that in some regions have the same order of similar magnitude than the unadjusted trend. This points out that O 3 -related mortality have been strongly affected by recent meteorological changes and indicates the need to perform regional and seasonal specific analysis to better understand the influences of meteorological factors on this secondary pollutant.

Pollutant
Finally, our NO 2 findings suggest that weather changes may have lessened the impact of NO x emission abatement policies in Spain, with a penalty of 14 deaths per year (95% CI −3, 31). Nonetheless, our assessments point out that the reductions of NO 2 levels may have had a significant effect in terms of long-term mortality, around −164 cases a year (95% CI: −300, −28).
Of note, the results at the regional level may differ substantially from the aggregated trends at the national level. As illustrated in Fig. 7a (number of deaths standardized by population), the PM 10 weather penalty ranged widely from more than 30 deaths annually (per 1 million people) in Aragon (AR) to a negligible influence in regions such as Madrid (MD), Murcia (MU) or Navarra (NF).
Regional differences are even more evident for O 3 (Fig. 7b), due to contrasting seasonal trends. While O 3 related mortality substantially decreased in areas such as Aragon (AR), Catalonia (CT) or Madrid (MD), it increased prevalence in Galicia (GA) and Navarre (NF). These regional differences in weather penalties highlight the need to take into account weather influences for each region specifically. As for NO 2 , Fig. 7c shows that relative penalties are in all cases small in comparison with those obtained for PM 10 . Consistent with those for this pollutant, the results for the actual (unadjusted) NO 2 trend indicate that the highest relative health benefits from recent pollution trends are found in the Asturias (AS) region.

Limitations
There are several limitations to our study that should be considered for a correct interpretation of the results and for the design of further research. Since our main goal is to identify weather penalties on pollution series at regional and national level in Spain, we build our GAM models at the regional level (NUTS-2). Despite the computational implications of calculating individual trends for each monitoring site, it should be noted that quality-controlled, comparable meteorological data are not available for all of them. That is why we matched pollution series from 354 air quality monitoring sites to 47 AEMET weather stations based on proximity criteria. Although we used a very dense network of monitoring sites, there may be considerable differences within a given region. The bootstrap analysis approach used allows estimation of confidence intervals that reflect site-to-site heterogeneity. Our results show that on a regional scale, the confidence intervals of weather penalties were smaller than air pollution trends and penalties. This indicates that while site-to-site variation of air pollution trends may be larger, the impacts of weather changes on trends are less variable between sites, and subsequently, within regions. As for the national trends, the random-effects model used to perform the meta- Fig. 6. Annual mortality changes in 1993-2017 as a result of unadjusted and weather-adjusted trends in PM 10 , O 3 and NO 2 and corresponding weather penalties from national meta-analysis. The 95% confidence intervals are shown. Table 3 Summary of average mortality changes (cases⋅yr −1 ) according to pollution trends and weather penalties at national level (meta-analysis) over the period . The figures in brackets represent the 95% confidence interval.  Fig. 7. Population-standardized (deaths per 1 million people) annual mortality changes in 1993-2017 as a result of unadjusted and weather-adjusted trends in PM 10 , O 3 and NO 2 and corresponding weather penalties at regional level. The 95% confidence intervals are shown.
R. Borge, et al. Environment International 133 (2019) 105272 analysis assumes that average effect size in the population varies randomly from region to region. Since we found strong differences between regions, we carried out some tests to detect potential biases. The results, included in the supplementary material, indicate that most of the variance in the meta-analysis comes from heterogeneity between regional trends, which may result in biased results for the meta-analysis of weather penalties at the national level in some cases (SO 2 both cold and warm season and PM 10 warm season). This highlights the need to carry out more detailed studies at the regional level. The methodology presented in this contribution demonstrates, however, the ability to support such studies. There are some limitations regarding the health assessment related to the trends and weather penalties too. We assumed that the entire population (over 30) according to the official census (Table S1) was exposed to air pollution, i.e. regional trends are representative for the whole population within that region. Although this cannot be taken for granted, we deem this hypothesis acceptable given the high spatial density of observations used and the fact that 74% of the monitoring sites used in this work are located in urban or suburban areas. Therefore, our trends are mainly related to urban trends, were population concentrates, so we think they may be representative in terms of exposure. This is also consistent with the relative risks derived from epidemiological studies that do not account for intra-urban variation of pollution levels (Hoek et al., 2013). A very recent study (Picornell et al., 2019) carried out in Madrid metropolitan area compared the exposure estimates according to the classic methodology that associates population and pollution according to their residence (static approach) with those obtained from massive cell phone data (CDRs) (dynamic approach) in the region. Although significant differences exist at district or neighbor level, the aggregated results for the metropolitan region were very similar, confirming that this approach is valid for the scale of analysis in this study.
As for the methodology used to estimate mortality changes, it is widely recognized that worldwide PM 2.5 exposure comprises the majority of air pollution mortality (Pope and Dockery, 2006;Kim et al., 2015;Burnett et al., 2018) and has been taken as a reference in previous air pollution health impact assessments in Spain (Boldo et al., 2011). However, historic data on fine particles is limited in Spain (as shown in Table 1) since routine measurements of PM 2.5 did not begin until 2001. Although representative PM 2.5 datasets are available for some regions (such as Madrid), our analysis relies on PM 10 trends to improve data availability encompassing all the Spanish regions. Finally, the relative risks used in the concentration-response functions chosen are not specific for Spain due to the lack of national cohort studies that could support the analysis of long-term pollution effects. Nonetheless, we clearly presented our methods and input data so any further study can contrast the results by changing the underlying assumptions (Malmqvist et al., 2018).

Conclusions
Ambient air quality in Spain has substantially improved over the last 25 years , with the exception of tropospheric ozone. For the same period, significant trends in meteorological variables, referred to as weather changes, have been observed. In this study, we report a general trend of increasing temperatures, especially in summer, as well as substantial decreases in wind speed. Humidity, evaluated through water vapor pressure, presents considerable changes in some regions although no significant trend was found at national level. In addition, we observed a generalized decrease in the number of rainy days in summer. We assessed the combined effect of these weather changes on trends in air quality by applying general generalized additive models (GAMs). While it is difficult to discern the influence of individual weather parameters on pollution trends, our models provide a robust combined estimate of the non-linear relationships among weather parameters and pollution levels. We found that air quality would have improved even more during our study period if meteorological conditions had remained constant. We also found that the impact of weather changes differ greatly by pollutant, suggesting that the influence of weather changes may vary widely from transport phenomena such as advection or diffusion to atmospheric chemical reactions governing secondary pollutant formation. We also found significant variability in weather penalties by region.
Our results suggest that despite weather penalties, recent emission abatement efforts in Spain have been successful in reducing air pollution-related mortality. We found that reductions in ambient PM 10 concentration have reduced the number of deaths attributable to air pollution over the 25-year study period by 33,047 (95% CI: −13,584, 52,510). According to our results, the benefits from air quality improvements in Spain regarding this pollutant would have been ap-proximately10% larger if weather conditions had not changed during this time period. We also found non-negligible weather penalties for O 3 and NO 2 .
Our study suggests that when policies to further improve air quality are designed they should consider not only the current meteorology, but also the changes in weather likely to result from climate change. If current trends continue, larger emissions reductions may be needed to achieve the same air quality goals in the face of weather penalties.
The regional differences found in this study stress the need for specific regional (and seasonal) studies were more detailed consideration about weather trends and emission source apportionments for different pollutants can be taken into consideration. To provide a better basis for policy developments, it will be advisable to complement such studies with air quality simulations based on deterministic models that can explicitly consider the relationships between air quality and weather parameters in a climate change context.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.