Changing PM2.5 and related meteorology over India from 1950–2014: a new perspective from a chemistry-climate model ensemble

Surface PM2.5 concentrations in India have increased dramatically as emissions have risen in recent years. The role of meteorological factors in this increase is unclear, mainly due to a lack of long-term observations over the region. A 12-member ensemble of historical (1950–2014) simulations from the Community Earth System Model version 2-Whole Atmosphere Community Climate Model version 6 (CESM2-WACCM6) offers an unprecedented opportunity to examine simulated PM2.5 and meteorology for 20th century climates that can arise due to ‘climate noise’ under the same historical greenhouse gas and air pollutant emission trajectories. CESM2-WACCM6 includes interactive aerosol and gas-phase chemistry in the atmosphere coupled to ocean-sea ice-land models, and each ensemble member differs only in its initial conditions of the climate state. We systematically examine, decade-by-decade, the changes in PM2.5 and associated meteorology, including wind speed, surface temperature inversions, boundary layer height, precipitation, and relative humidity in four cities in India: Chennai, Kolkata, Mumbai, and New Delhi. Forced changes clearly emerge in meteorological variables from 1950 to 2014, including increases in both relative humidity and temperature inversion strength, and decreases in boundary layer height and average surface wind speed. The timing of these changes varies by city: boundary layer heights decrease most over New Delhi in the premonsoon season (ensemble average decrease of 400 m), but over Mumbai in the postmonsoon season (ensemble average decrease of 100 m). PM2.5 concentrations increase across India regardless of climate variability, with an almost threefold increase from 1950 to 2014 over New Delhi. Analysis of dimensionless variables shows that PM2.5 exhibits larger ensemble mean trends and smaller variability than the trends in the meteorological variables, enabling us to infer that the increase in PM2.5 is predominantly controlled by increases in anthropogenic emissions rather than climate variability. Overall, our simulations corroborate the dominant role of air pollutant emissions on poor air quality in India.


Introduction
In India, fine particle (PM 2.5 ) concentrations have been rapidly increasing; in New Delhi, the 2005 WHO daily mean PM 2.5 guideline of 25 µg m −3 [1] was exceeded on more than 90% of days from 2014 to 2019 [2], with 980 000 deaths in India attributed to ambient particulate matter pollution in 2019 [3]. In 2021, the 2. Methods

Model simulations
CESM2-WACCM6 is a climate model that includes full atmospheric chemistry in the troposphere and stratosphere [31,32] coupled to a modal aerosol scheme [33]. It is higher-resolution than prior atmospheric chemistry models coupled to full ocean-sea ice-land models, with grid cells of size 0.95 • latitude × 1.25 • longitude, and 70 vertical levels. We archived fields for PM 2.5 , aerosol components, and meteorological variables; some were archived daily, while others are available only monthly or annually. The model is driven by anthropogenic emissions of air pollutants and their precursors, and concentrations of long-lived greenhouse gases from 1950 to 2014. It represents the interactions between chemicals in the different layers of the atmosphere with each other, and the interactions of the atmosphere with the ocean, sea ice, and land vegetation to simulate weather and atmospheric composition in each grid cell. Note that a climate model is not expected to match the exact weather at any specific location or time, but rather, should simulate statistics similar to those observed in the actual atmosphere (and the broader climate system).
For each of the 12 individual ensemble members, the model ingests the same greenhouse gas and air pollutant emission scenarios but uses a unique set of initial conditions. Anthropogenic emissions are from the Community Emissions Data System (v2016-07-26) [34]. Biomass burning emissions are from the Biomass Burning for CMIP6 (BB4CMIP6) database [35], and greenhouse gases are from Meinshausen et al [36]. Dust and sea salt emissions are calculated online and depend on wind speed. This is described in Emmons et al [31]. Secondary aerosol formation from biogenic emissions involves a sophisticated approach that includes the volatility basis set scheme, whereby hundreds of intermediate semivolatile organic compounds are grouped by bins of volatility that are produced by oxidative chemistry from emitted precursors. These mechanisms are described in detail in Tilmes et al [32], which also explains that dimethyl sulfide and organic aerosol emitted from oceans are included as part of WACCM6's comprehensive chemistry.
We generate our initial condition ensemble from 1 January 1950 using three ensemble members that were launched from different years within a long pre-industrial control simulation that were extended using historical forcings (long-lived greenhouse gases and emissions of air pollutants and precursors). These ensemble members have different initial conditions for all Earth system model components (ocean, atmosphere, etc.). Additional ensemble members are generated from each of these three initial conditions by slightly perturbing the atmosphere initial conditions (more details about the ensemble can be found in Fiore et al [37]; our current ensemble is 12 members). These ensemble members represent 12 different possibilities for trends in atmospheric composition and meteorological variables from 1950 to 2014, with each ensemble member representing a climate system state that could have arisen under the same forcing scenario. While we lack historical observations for PM 2.5 and other meteorological indicators over India, these simulations provide a range of trends that might have been observed, at least under the assumption of minimal model error.

Data analysis
We examine simulated trends of meteorological variables and PM 2.5 across India. We focus on meteorological variables identified by Schnell et al [38], who found that high PM 2.5 events over the Indo-Gangetic Plain (the most densely populated area of India along the northern border, referred to as the IGP) during winter are associated with high surface relative humidity, high temperature inversion strength, low planetary boundary layer height, and low surface wind speed. Temperature inversion strength is calculated as the difference between the 850 mb and surface temperature such that a less negative value indicates a stronger inversion. For calculating seasonal trends, we define seasons based on India's seasonal meteorology described by David et al [39]: winter (December-February), premonsoon (March-May), monsoon (June-September), and postmonsoon (October-November).
To infer the extent to which trends in meteorological variables, rather than emissions, might influence trends in PM 2.5 , we calculate dimensionless versions of PM 2.5 and each meteorological variable. Standardized versions of each variable allow us to directly compare trends and variability across the ensemble in a common frame of reference; if meteorological variables were dominating trends in PM 2.5 , we would expect to see similar variability in both. For each ensemble member and grid cell over India, we calculate the mean (µ n,g where n denotes the ensemble member and g denotes the grid cell location, a (lat,lon) pair) and standard deviation (σ n,g ) over the whole time series. We then subtract the temporal mean and divide by the standard deviation for each individual annual mean PM 2.5 value at year t (t = 1950, 1951, …, 2014): d n,g,t = x n,g,t − µ n,g σ n,g .
From each dimensionless time series (d n,g ) for each ensemble member n and grid cell g, we then compute σ ′ g,t (the standard deviation across the ensemble) for each year t in the time series, where σ ′ g,t is the standard deviation of {d 1,g,t , d 2,g,t , . . . , d 12,g,t }. We define the spread across the ensemble as the time-averaged mean standard deviation, µ ′ σ , the mean of . We can then assess how this standardized spread compares between PM 2.5 and meteorological variables, as well as whether trends in dimensionless variables emerge from the noise (µ ′ σ ). We focus some of our analysis on grid cells encompassing four cities in which high-quality U.S. Embassy and Consulate observations are available for model evaluation: Chennai, Kolkata, Mumbai, and New Delhi. Examining these cities found in four different regions of India allows us to analyze separate regional trends within India. Monthly observations of surface PM 2.5 between 2013 (before which data is unavailable) and 2015 were obtained from the U.S. Embassy in New Delhi and the U.S. Consulates in Chennai, Hyderabad, Kolkata, and Mumbai. Since the model simulations only extend to 2014, we compare these observations to the simulated years 2012-2014. The inclusion of three years, rather than simply the two overlapping years, is motivated by our interest in evaluating climatological patterns, particularly as a climate model will not generate weather patterns observed on a particular day but should represent the statistics of observed weather. The scale of a model grid cell is coarse (0.9 • × 1.25 • ) compared to the point measurement and should be interpreted as simulating regional average conditions. While comparing model grid cell values to in situ point measurements is not ideal due to the coarse model resolution and the localized nature of in situ measurements, it is common practice in global atmospheric chemistry model evaluation in light of the lack of regionally representative measurements [38,40]. We also compare modeled values to measurements of surface PM 2.5 in 74 cities averaged from 2012 to 2017 from the Central Pollution Control Board of India (CPCB) [41] and in observations in 2 additional cities obtained from a literature search [42,43].
For examining changes in annual average PM 2.5 over time, we use an online PM 2.5 model diagnostic (text S1). As we show below, the seasonality of the online model PM 2.5 diagnostic does not match observations, due to excessive dust and sea salt contributions to the smaller size bins. We thus construct a PM 2.5 estimate as the sum of individual PM 2.5 components in the Aitken, accumulation, and primary carbon modes [44], which includes black carbon (BC), particulate organic matter (POM), sulfate, sea salt, dust, and secondary organic aerosol (SOA) components. The component sum corresponds to particle diameter sizes of 0.02-0.08 µm for the Aitken mode and 0.08-1 µm for both the accumulation and primary carbon modes in MAM4 [44][45][46]. Monthly values are available for all components except for SOA, which was only saved annually in the model due to its detailed representation [32]. We use this component sum to examine the changing seasonal cycles of PM 2.5 by decade, assuming a uniform contribution from SOA throughout the year. For examining changes in individual components over time, we continue to examine only the Aitken, accumulation, and primary carbon modes to remove the outsized effect of coarser particles present in the model diagnostic. Neither the PM 2.5 diagnostic nor the component sum contain aerosol water content, which may contribute to underestimates.
Given the limited availability of high-quality surface PM 2.5 measurements, we also evaluate the model with monthly aerosol optical depth (AOD) retrieved from satellite from 2003 to 2014. Although AOD is often used to derive estimates of ground-level PM 2.5 [47][48][49], we compare modeled 550 nm AOD values directly to the AOD measurements from moderate resolution imaging spectroradiometer (MODIS) level 3, collection 6.1 instruments on the Aqua and Terra satellites. Satellite AOD values were obtained from Giovanni [50] by selecting the 1 • × 1 • grid cell that encompasses each city. We use the Dark Target algorithm which has been shown to compare better than Deep Blue to both ground based AERONET AOD and higher resolution MAIAC retrieval over India [51].  some cities is not observed. The difference is most notable in Mumbai, which has an observed winter peak and minimum in July-August, but a modeled March-April minimum and a July-August peak.

Model evaluation
We attribute the modeled monsoon season peak to excess larger particles (dust and sea salt) included in the PM 2.5 diagnostic, which will be corrected in newer model versions. When we instead calculate the sum of available PM 2.5 components in the Aitken, accumulation, and primary carbon modes [44] (see section 2 for details), the seasonal cycle aligns more closely with observations in Kolkata, Mumbai, and New Delhi (blue curves in the top row of figure 1), although the component sum worsens the underestimate of PM 2.5 in all cities (likely because of the omission of coarse particles, as well as nitrate and ammonium, which have been shown to be important components of PM 2.5 in India [38]). In Chennai, both the component sum and the model diagnostic have peaks in May and October that are not present in the observations. Black carbon, dust, sea salt, particulate matter, and sulfate all show these peaks in Chennai (figure S1). Note that this component sum does not capture the seasonality of SOA, since SOA is applied uniformly throughout the year. SOA has been shown to peak during winter over India [52][53][54][55], which matches the overall PM 2.5 component sum seasonal cycle, so the assumption of constant SOA likely decreases the amplitude of the overall seasonal cycle, with an underestimate of winter SOA and an overestimate of monsoon season SOA. However, given the magnitude of SOA concentrations relative to the total PM 2.5 sum, we do not expect SOA seasonality to substantially influence the alignment of the component sum with observations.
Although we do not expect AOD (column light extinction by all aerosols, both fine and coarse) and surface PM 2.5 seasonal cycles to match since their relationship varies with aerosol composition and meteorology [56,57], AOD retrieved from satellite offers another observational test of the CESM2-WACCM6 representation of aerosols. The second row of figure 1 shows average monthly AOD from 2003 to 2014 in the model (red line for ensemble average and shading for the range) and from the MODIS instruments aboard the Aqua (black line, morning overpass) and Terra (green line, afternoon overpass) satellites. Aqua and Terra values are similar in all cities excluding Chennai, in which the AOD retrieved from Terra is slightly higher than from Aqua in all months except September, October, and November. In all cities, modeled AOD falls close to the values retrieved from MODIS. In Mumbai, the simulated and observed AOD seasonal cycles match almost exactly, with a July peak of between 0.8 and 1 and a winter minimum. In New Delhi, the satellite and model AOD also share a July peak, but the satellite products show a higher winter AOD (model AOD <0.4 in January, while satellite AOD >0.7). In Chennai and Kolkata, satellite values exceed modeled values. In Kolkata, the model and satellite AOD both have a June peak, but the satellite AOD shows a second peak in August that is not captured by the model. In Chennai, the model peaks in May and The bottom row of figure 1 shows the spatial distribution of AOD across India, averaged from 2003 to 2014. The left panel shows AOD from Aqua (we omit Terra for this analysis because values between Aqua and Terra are closely correlated, leading to the same conclusions), while the middle panel shows the CESM2-WACCM6 ensemble average. CESM2-WACCM6 values are slightly lower than the satellite values, especially over the IGP where AOD is highest. However, the spatial patterns of AOD generally align, with high values over the IGP, lower values in central India, and lowest values in Southern India. The most notable anomaly is over Western India near the Gujarat region, where the model shows elevated AOD concentrations that are not present in the observations. This discrepancy may arise from a combination of excessive sea salt and dust. Tilmes et al [32] noted a global discrepancy between CESM2-WACCM6 and observed AOD which they hypothesize is from excess sea salt, and dust is transported to this coastal region during the pre-monsoon season [58]. The right panel shows the spatial correlation between Aqua and the ensemble average grid cell values over India, averaged from 2003 to 2014. While the coefficient of determination is relatively low (r 2 = 0.29), the root mean square error (RMSE) and mean bias error are low (0.125 and −0.025 respectively). Over the IGP, Aqua retrieves AOD values ranging from 0.5 to 0.8, while the model rarely simulates an annual mean AOD value greater than 0.5.
Precipitation is a major sink for PM, so we evaluate its seasonal representation in CESM2-WACCM6. The spatial variability and seasonality of precipitation align (spatial RMSE of 2.7 mm d −1 ) with data from the Global Precipitation Climatology Project, as shown in figure S2.
We further evaluate the model using near-surface CPCB measurements of PM 2.5 in 74 cities and 2 additional sets of PM 2.5 measurements obtained with a literature search ( figure 2). We compare 2012-2017 averages to CESM2-WACCM6 values for the grid cells encompassing each city, averaged from 2009 to 2014 (blue points). We also compare measurements from Murari et al [42] and Joseph et al [43] Figure 3 shows the spatial patterns and trends from 1950-2014 in the anthropogenic plus biomass burning emissions that are used to drive the CESM2-WACCM6 ensemble simulations. We lump together all primary PM emissions, defined as the sum of the surface flux of primary sulfate, black carbon, and POM. We use SO 2 , the precursor gas to sulfate aerosol, to represent changes in the emissions most relevant to anthropogenic

Changes from 1950 to 2014
In examining the spatial distribution of PM 2.5 and related meteorology over India, we find some general trends emerge across the entire ensemble. Figure 4 shows the change of surface PM 2.5 , relative humidity at 2 m reference height, temperature inversion strength, planetary boundary layer height, and average surface wind speed from 1950-1959 to 2005-2014. For each ensemble member, for each variable, in each grid cell, we subtract the 1950-1959 average from the 2005-2014 average. The minimum, average, and maximum values taken over all the ensemble members were selected for each grid cell to demonstrate the variability in these trends between the ensemble members. We use the same criteria for significance testing as Fiore et al [60], whereby changes are considered significant if the ensemble mean change exceeds the 95th percentile of a sample constructed from differences in 1950-1959 means across pairs of ensemble members.
The PM 2.5 increases in figure 4 generally follow the spatial patterns of the precursor emissions shown in figure 3, with the highest emissions-and increases in PM 2.5 -across the IGP. The spatial distribution of Overall, on average, relative humidity and temperature inversion strength increase across India while boundary layer height and average surface wind speed decrease. The magnitudes of these trends, however, vary across ensemble members. While forced changes in meteorology over India occur (particularly in boundary layer height), climate variability has also played a role in shaping trends in these meteorological variables from 1950 to 2014. We demonstrate this further in the following section.

Emission versus meteorological drivers of PM 2.5 trends
In an effort to place different variables into a common framework for comparison, we non-dimensionalize all variables and examine their standard deviations across ensemble members, calculated as described in section 2.2 ( figure 5). If the meteorological drivers were dominating the trend in PM 2.5 , we would expect to see similar variability in meteorology and PM 2.5 . The spread of dimensionless PM 2.5 across the 15-member ensemble is low (<0.3) compared to the other variables (>0.7) throughout India. We interpret the much smaller variability in PM 2.5 , along with larger ensemble mean trends (slopes; see also section 6.1), as indicating that PM 2.5 is predominantly controlled by the rise in anthropogenic emissions (figure 3), which are identical in each ensemble member, leading to a lower spread in the simulated PM 2.5 concentrations. Conversely, the higher spread of the meteorological variables demonstrates that they are influenced by climate variability to a much larger extent than PM 2.5 , even though we are able to detect some forced changes in meteorology ( figure 4). This confirms the dominant role of rising anthropogenic emissions, rather than natural climate variability, in raising PM 2.5 concentrations over India over the last 65 years.
The PM 2.5 spread is highest over Western India in the same location that modeled AOD values are unexpectedly high. This supports the idea that the excessive AOD could reflect the influence of dust or sea salt, which vary strongly with meteorology and thus across the ensemble, in contrast to PM 2.5 components dominated by anthropogenic emissions which are identical across the ensemble. The PM 2.5 spread is especially low over the IGP, where anthropogenic emissions are highest (figure 3).
Relative humidity and boundary layer height also show a lower spread over the IGP. This indicates that while meteorology is more affected than PM 2.5 by climate variability, in areas where emissions are high, meteorology may also be playing a stronger role in contributing to the forced trends in PM 2.5 . Boundary layer height especially shows a lower spread (as low as 0.3) over the IGP, which may imply a direct response to the rise in local aerosol emissions [14,61]. Wind speed shows a more uniform spread across India, with slightly lower standard deviations further from coastlines. The spread of inversion strength is generally uniform and high (0.9-1.0) as well but shows some lower standard deviations (0.7-0.8) over the IGP and in Southern India, also possibly suggesting a role for a local meteorological response to anthropogenic emissions. While modeled high aerosol loading has been shown to impact meterology in India [62,63], additional simulations would be needed to conclusively identify a causal relationship here.

City-specific changes
In this section, we zoom into four selected cities. We first probe more deeply the changes in individual components. Then, we explore changes in the seasonality, and apply our non-dimensional analysis to infer the relative roles of anthropogenic emissions versus climate variability on the PM 2.5 changes from 1950-2014.  We also examined the role of meteorology on PM 2.5 trends over India within each city by calculating trends in dimensionless meteorological variables (figure S2). The overall trends are the same as those shown in figure 4 (decreasing wind speed and boundary layer height, increasing strength and relative humidity in all cities), but even when these changes occur in all ensemble members, the variability across the ensemble for meteorological variables is much larger (2-3 times larger average maximum-minimum values) compared to that of PM 2.

Shifting PM 2.5 and meteorology seasonal cycles
Seasonal cycles of PM 2.5 differ in the four cities. In all grid cells, seasonal cycles shift as winter PM 2.5 levels more than double from 1950 to 2014. Figure 7 shows the average monthly PM 2.5 component sums for every other decade from the 1950s to 2010-2014 for Chennai, Kolkata, Mumbai, and New Delhi. Each line represents the average over all ensemble members while the shading represents the range across the ensemble members for a given decade (note the last period is only half a decade, 2010-2014). Kolkata, Mumbai, and New Delhi all show a small May-August peak in PM 2.5 with a slightly larger winter peak in the 1950s. By the 2010s, the winter peak eclipses the initial May-August peak in these cities. Chennai, however, experiences peaks in both March-May and October in the 1950s that increase along with winter PM 2.5 levels. Monsoon season PM 2.5 levels in all cities increase each decade, with overlap of ensemble members between some decades but a clear forced increase overall, while winter PM 2.5 levels increase, with little overlap between decades, indicating significant increases between decades for all ensemble members. The range across the ensemble also increases with each decade, implying that PM 2.5 decadal averages are more influenced by climate variability now than in the mid-19th century. This is confirmed by dimensionless analysis: the dimensionless spread across the ensemble increases by 60.17%, 56.73%, 19.42%, and 20.04% from 1950-1955 to 2010-2014 in Chennai, Kolkata, Mumbai, and New Delhi, respectively (table S1). Higher PM 2.5 concentrations during these more recent periods may be more susceptible to meteorological influence, as in Chen et al [6]. For 2010-2014, the use of a five year rather than ten year average may be inflating the range, but even ignoring 2010-2014, the range increases across the ensemble over time.
Seasonal cycles of meteorological variables also vary by city and by decade. Figure 8 shows the seasonal cycles in the first and last decades of the simulation (1950-1959 and 2005-2014)  Forced increases are also detected in inversion strength, with no overlap between ensemble members in parts of the monsoon season in Mumbai and in the premonsoon season in New Delhi. Kolkata and Chennai show similar seasonal cycles in inversion strength (winter low, premonsoon season high) without a clear forced change between decades. As expected, precipitation is at its maximum in all cities during the monsoon season. The ensemble average precipitation during the monsoon season increases in all cities from 1950-1959 to 2005-2014, but there is overlap between the ensemble members between decades in all months, indicating that differences in precipitation reflect climate variability rather than a forced change. All cities show an ensemble average increase in relative humidity between the first and last decade, with no overlap between ensemble members between decades during many months. This finding indicates a forced increase in relative humidity. The increase is largest in New Delhi, where average relative humidity increases from 20% to 35% in May. Overall, significant changes occur in all meterological variables examined except precipitation with different seasonal timing between cities: with no ensemble overlap between decades, decreases in average surface wind speed occur in New Delhi during the premonsoon and monsoon seasons, boundary layer height decreases in all cities in different seasons, inversion strength increases in Mumbai and New Delhi in different months, and relative humidity increases in all cities.
The forced changes in meteorology, also shown in section 5.1, could be responses to the local aerosol increases, or to global increases in aerosols and/or greenhouse gases, or other forcings during 1950-2014 to which meteorology responds [9]. The forced changes detected in figure 8 are consistent with a response to local aerosol increases, which can stabilize the lower atmosphere [64,65], although again, this does not conclusively identify a causal relationship. Differentiating between different mechanisms (globally forced changes in meteorology impacting aerosols versus local aerosols impacting meteorology) could be done in the future through a parallel set of model simulations that keep aerosol emissions over India at 1950 levels but allow everything else to change. Future air quality planning would benefit if the underlying drivers were identified.
Regardless of their origin, these meteorological changes exacerbate air pollution even in the absence of emissions growth. Moreover, it should be noted that the timing of the forced changes, and thus the seasons that air pollution worsens in response to these changes, differs between cities and meteorological variables. For example, boundary layer heights decrease most in New Delhi during the premonsoon season, while the largest decreases in Mumbai occur postmonsoon. Consideration of these variations can help inform air quality planning.

Conclusions
A newly available set of CESM2-WACCM6 model ensemble simulations with interactive tropospheric chemistry and aerosols show that PM 2.5 levels rise across India from 1950 to 2014 regardless of meteorology differences, corroborating the dominant role of air pollutant emissions, rather than climate variability, on poor air quality in India. To the best of our knowledge, previous work has not considered the role of internal variability on trends in PM 2.5 and meteorology over India. The newly available historical ensemble simulations analyzed above enable us to rule out climate variability as the driver of the increasing PM 2.5 trend. Trends in dimensionless PM 2.5 and meteorological variables strengthen this conclusion: steep PM 2.5 increases with low variability (inter-ensemble spread <0.3) relative to meteorological variables (inter-ensemble spread >0.7) indicate that PM 2.5 is predominantly controlled by the rise in anthropogenic emissions, which are identical across all ensemble members.
Our selection of four cities in different regions of India revealed that the dominant PM 2.5 component shifts from primary organic matter in 1950 to sulfate in Chennai and Mumbai, and SOA in Kolkata. While CESM2-WACCM6 captures observed inter-city variability and winter peaks in PM 2.5 , it underestimates winter PM 2.5 concentrations.
Meteorological variables related to the accumulation and ventilation of pollution vary across the ensemble members due to natural (internal) climate variability, but long-term trends forced by the same combination of rising greenhouse gases, anthropogenic and biomass burning pollutant emissions in each ensemble member clearly emerge. For example, relative humidity and temperature inversion strength increase while boundary layer heights and average surface wind speed decrease across India from 1950 to 2014. The stronger trends with smaller inter-ensemble ranges over high-emission areas suggest a possible meteorological response to increases in local aerosol emissions. New model sensitivity simulations that keep aerosol emissions at 1950 levels in India but allow other emissions to change as prescribed in these simulations would be necessary to separately quantify the roles of globally forced meteorological changes versus local meteorological responses to rising aerosols.