Greenhouse gases stabilizing winter atmosphere in the Indo‐Gangetic plains may increase aerosol loading

The concentration of aerosols along the Indo‐Gangetic Plains (IGP) and their adverse effects on human health and the environment are increasing, particularly in winter, when the IGP is prone to high anthropogenic aerosol loading (i.e. particulate matter – PM). In this study, three reanalysis datasets including the MERRA aerosol reanalysis were analyzed to characterize the relationship between winter atmospheric stability and aerosols. Due to the lack of long‐term aerosol observations, an empirical relationship between PM and the atmospheric temperature lapse rate was derived. It is shown that PM and stability have a strong relationship at the lower troposphere. Analyses of Coupled Model Intercomparison Project (CMIP5) single‐forcing experiments indicate that the recent stabilization of the atmosphere in the IGP can be explained by both greenhouse gas (GHG) and aerosol forcings. Since a more stable atmosphere traps more PM, stabilization alone can increase aerosols even in the absence of changes in emission sources. Regional simulation with fixed anthropogenic aerosol loading also supports this finding. Thus, enhanced stabilization caused by both aerosols and GHG in the atmosphere can further increase winter aerosol loading in the IGP.


Introduction
In recent decades, the Indo-Gangetic Plains (IGP) has become a hotspot of rising anthropogenic aerosol emissions, as evidenced by measurements of moderate resolution imaging spectroradiometer (MODIS) aerosol optical depth (AOD) - Figure 1(a), and remarked upon by various studies (Vijayakumar et al., 2007;Kar et al., 2010). The increased concentration of anthropogenic aerosols, coupled to their bearing on climate change and effect on human health, has resulted in a marked proliferation of research into atmospheric aerosols. Past studies (Gautam et al., 2011;Srivastava et al., 2012) have found that the concentration of anthropogenic aerosols in the IGP is not as straightforward as a function of source, but is complicated by topography and seasonality; this is akin to what occurs with particulate matter (PM) (e.g. PM 2.5 ) in many terrain basins (Gillies et al., 2010;Wang et al., 2012).
The winter months within the IGP feature a stable boundary layer due to cool temperatures and weak ambient wind flow [November-January, Figure 1(c)], hence the potential for high-aerosol loading; this being contrary to the summer monsoon season where high temperatures accompany strong monsoonal winds [May-July, Figure 1(b)]. Aerosol measurements, using the angstrom exponent of aerosols that increases with decreasing size ( Figure S1), reveal that the make-up of atmospheric aerosols in the IGP do indeed vary seasonally -i.e. smaller aerosols (e.g. black carbon) are predominant in winter and larger ones (e.g. dust) in summer. The large extent of human induced (anthropogenic) aerosols in winter therefore underscores the importance of studying any corresponding changes in the local climate. Moreover, when it comes to understanding the IGP's aerosol sensitivity to the changing atmospheric thermal structure driven by anthropogenic global warming, a lack of systematic research coupled to the paucity of a long-term observational aerosol data record is problematic; this motivated the research undertaken and presented here.
In this study, we utilized atmospheric and aerosol reanalysis data in conjunction with climate model outputs to characterize atmospheric stability and examine the relationship between atmospheric stability conditions and anthropogenic aerosols. The diagnostic analyses were followed by the derivation of an empirical relationship between atmospheric aerosol loading and stability, which subsequently was used to quantify the concentration of aerosols in the IGP, and how it interacts with changing atmospheric stabilization.  University of Delaware (UDEL) precipitation datasets (Willmott and Matsuura, 2001), NASA's MERRA (Rienecker et al., 2011), and the NCEP/NCAR 40-year reanalysis datasets (Kalnay et al., 1996). For attribution, 11 models that participated in the Coupled Model Intercomparison Project (CMIP5) (Taylor et al., 2012) were used. These single-forcing experiments only go up to 2005. Table S1 details the acronym, full name, and full description of each model. Data from the space borne MODIS instrument (Hsu et al., 2004) and the Aerosol Robotic network (Holben et al., 1998) were used for the depiction of aerosol loading and properties. We also examined PM derived from the MERRA Aerosol Reanalysis (MERRAero) (Rienecker et al., 2011) from 2003 to 2013. In MERRAero, PM can be obtained from the surface mass concentrations of dust, sulfate (SO), black carbon (BC), organic carbon and sea salt. Over the IGP, SO and BC emissions are significantly higher in the winter than in the summer; this is due to differences in weather and industrial activity. This article focuses on the relationship between 'anthropogenic aerosol' emissions and stability for the period defined as 'winter', therefore dust was excluded from the analysis to minimize possible biases.

Methodology
The study area along the IGP is outlined in Figure 1(a); this area features particularly large aerosol loadings due to its high-population density. Given the maximum PM in the seasonality of the IGP, we focused on the winter months of November and December (ND).
To depict atmospheric stability, we used lapse rate, as it is directly relevant to vertical motion and, subsequently, the dispersion and concentration of aerosols in the atmosphere. Both the actual (T) and potential temperature ( ) were used to compute lapse rate (i.e. -d [T]/dp). To make compositing possible, all utilized reanalysis datasets were regridded onto the same 2.5 ∘ grid spacing using a third order Bessel interpolation. Attribution involved repetition of the lapse rate analysis with different forcing composites of CMIP5 model outputs. To test the evolution of PM concentrations in the IGP without time-varying anthropogenic aerosol emissions, the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) version 3.8.1 (Grell et al., 2005) was used. The experimental setup, main physics and chemical schemes used for the simulations are listed in Table S2.

Results and discussion
3.1. Instability Figure 2 shows the changes in the atmospheric stability displayed as vertical profiles of lapse rate anomaly as a function of time. Linear trends at each pressure level, computed as the least square best fit to the lapse rate over the entire underlying period of 35 years, are overlaid in contours. The climatological mean of the lapse rate was removed prior to computations to depict year-to-year changes. As a consequence of removing the long-term mean, one is able to observe two distinct periods in variability, i.e. one in the past centered around the 1990s that is indicative of a period that is less atmospherically stable, compared to the other centered around 2010 where the situation is reversed. All the three reanalysis datasets (Figures 2(a)-(c)) indicate the development of increasingly atmospheric stable conditions from the surface to about 700 hPa. Also evident is a similar occurrence in the middle-to-upper troposphere with the exception of NCEP1 (Figure 2(a)). Regardless of marginal distinctions between the three observational datasets, a consistent picture arises from the composite plot ( Figure 2(d)), one that signals a stability modification for the entire troposphere. In other words, a shift has occurred, one from a climatologically stable winter atmosphere environment to one of greater stability; this implies extended or strengthened winter temperature inversions over the IGP.
To attribute the shift in the observed lapse rates, we analyzed a composite of 11 CMIP5 models. Due to the limited temporal extent of CMIP5 data (up to 2005), only a comparable 35-years time span (i.e. 1970-2005) with respect to the reanalysis datasets was used: Starting with the natural runs (Figure 2(e)), there is no clear-cut shift toward increased atmospheric stability as the lapse rate anomalies tend toward neutrality -this is not surprising due to the cancelation of natural variability among the different simulations over time. The greenhouse gas (GHG) forcing (Figure 2(f)) runs do, however, reveal a similar interdecadal shift in comparison with that in Figure 2(d) toward a more pronounced stable atmosphere in the middle-to-upper troposphere around 2005. The modeled GHG forcing is not particularly synchronous with the observations, as a more stable atmosphere develops much earlier than that observed despite the same direction. This is expected, however, since the GHG-only forcing would induce more pronounced warming than the observation, enhancing the stabilization effect earlier.
In the aerosol run (Figure 2(g)), an increasingly stable boundary layer is also present, with corresponding positive (and somewhat accelerated) trends in most of the troposphere. The combination of all three forcings (Figure 2(h)), i.e. the ensemble (ALL or historical experiment) scenarios, results in uniform stabilization throughout the tropospheric column, but once again, it reflects a disparity with the reanalysis composite (Figure 2(d)) in the timing of the shift in stability, i.e. about a decade earlier. Apparently, the models tend to also overestimate stabilization in the middle troposphere in comparison to the reanalysis. Nevertheless, the observations are clear as to the change that has occurred and the CMIP5 results suggest that there is an effect to enhanced stabilization that partly lies in anthropogenic GHG's and aerosols.
Caution is necessary here in interpreting the model results. The CMIP5 free runs cannot be compared to calendar years, and aerosol effects are only partially represented in the majority of the models (e.g. Kiehl, 2007). Even models that include atmospheric chemistry modules have shown large uncertainties in aerosols (e.g. Kinne et al., 2006).

Empirical assessment of aerosols
Since most CMIP5 models do not have built-in interactive or dynamic aerosol modules, and given that the time span of available aerosol data is limited, an empirical modeling approach serves as a means to estimate long-term regional aerosol levels. Hence, we derived a statistical relationship between PM and lapse rate using linear regression, following Wang et al., 2015. The choice of an optimal lapse rate level for temperature was achieved by averaging between the 925 and 800 hPa levels. Several tests with different instability indices ( Figure S2) revealed that any correspondence of PM with lapse rate diminishes somewhat beyond 700 hPa, consistent with the finding of Mishra and Shibata (2012). It is worth noting that the use of either actual or potential temperature does not change the relationship between PM and instability at lower altitudes. Figure 3(a) shows a scatter plot between PM anomaly and lapse rate at the optimal level for both actual (ΔT) and potential (Δ ) temperatures over time. Long-term trends were first removed to minimize anthropogenic sources and interdecadal variability. From the regression in Figure 3a, it is clear that PM [y(t)] and lapse rate [x(t)] exhibit a good linear relationship. Overall, the inter-comparison shows a close agreement between the two variables as indicated by the moderate-to-strong correlation. In fact, the large spatial extent of the datasets used, coupled to the relatively large size of the domain lend confidence to the strong statistical association.
Using the linear regression model (Figure 3(a)), it is possible to reconstruct ΔPM [y(t)] as a function of lapse rate [x(t)], where lapse rate is the proxy for inferring stability or otherwise. To do this, we used lapse rate data at selected pressure levels, extracted from the four CMIP5 forcing composites (Figures 3(b)-(e)). From the reconstruction, we can infer past PM levels and changes in PM concentration under different climate forcing experiments. In the HISTORICAL-forcing (Figure 3b) scenario, PM increases with time and appears driven by both increasing GHG (Figure 3(d)) and aerosols (Figure 3(e)). For the natural conditions (Figure 3(c)), there is little to no change. The most striking trend revealed is in Figure 3(d), which suggests that under GHG forcing, the concentration of aerosols in the IGP would continue to rise assuming PM emission sources remained constant -this is a possibility that was hitherto unknown.

Direct simulation of aerosols
One might recall that, apart from aerosols themselves, GHG's can also induce stabilization in the IGP The optimal level is 925-800 hPa -shown in red for potential temperature and blue for actual temperature, overlaid with their respective linear trends in the same colors.
level while GHG concentrations and warming continued. The simulated PM, which is shown in Figure 4, exhibits a clear increase even without the time-varying anthropogenic aerosols. The marked interannual variability also indicates that PM does undergo climatic modulations in the region. Both the climate effect and the GHG-induced warming can be seen in the simulated temporal evolution of temperature lapse rate anomaly in the IGP, plotted in Figure 4 alongside the simulated PM. The WRF-Chem simulation lends support to the notion that increasing GHG emissions as a singular factor can lead to a consequent increase in PM concentrations. In the same way, PM rise (Figure 3(e)) results in aerosol-induced stabilization (Figure 2(g)), which induces a positive feedback and so, traps more aerosols in the boundary layer -a known phenomenon (Lau et al., 2005). Another implication from the long-term trends is that GHG's have a longer atmospheric residence time and so, prolongs atmospheric stability conditions enhancing aerosol loading.

Other measures
Several variables were analyzed to show seasonal synoptic changes in the IGP, which primarily serve as validation for the outcome of the earlier lapse rate computations. Figure S3 shows the seasonal long-term changes (slope) in 800 hPa potential temperature. The atmosphere appears to have warmed up significantly in winter ( Figure S3(a)) in contrast to a net cooling in summer ( Figure S3(c)). Furthermore, we analyzed the slopes of winter precipitation and 2-m temperatures (T 2m ) within 2 different periods of time: precipitation is observed to have declined over the last 34 and 20 years ( Figure  S4(a)), along with surface cooling (Figure S4(b)).  -3 1980 1985 1990 1995 2000 2005 2010 1980 1985 1990 1995 2000 2005 2010 Table S2.
Concomitant to the rainfall decline is the hypothesis that aerosols increase the lifetime of clouds inducing an increase in the concentration of smaller droplets, which in turn lead to decreased drizzle production and reduced precipitation efficiency (Bollasina et al., 2008;Ganguly et al., 2012). The impact of aerosols on precipitation is by no means a straightforward conclusion. It is important to note studies such as Kim et al. (2016) that have associated the high concentration of aerosols over northern India to enhanced rainfall in the late spring and early monsoon months, followed by a suppression of monsoon rainfall over all of India. This is consistent with the elevated heat pump hypothesis by Lau et al. (2005). The vast extent of the deduction such as we have discussed illustrates the complexity of aerosol forcing on climate, and underscores the need for more research to be conducted in this area.
Also apparent is a reduction in T 2m ( Figure S4(b)), which together with the decreasing IGP precipitation are aligned with increased aerosol concentrations and enhanced stability conditions (Figures 2(a)-(h)), especially at lower altitudes. The surface cooling ( Figure  S4(b)) coupled to a warming trend aloft ( Figure S3(a)) is not consistent with a normal environmental lapse rate of decreasing temperature with height, and therefore reinforces the inferences of a more stable atmosphere. The results presented here suggest that the trapping of light absorbing aerosols (e.g. BC) through the enhancement of a thermal inversion layer suppresses boundary layer convection, and is consequential as it constrains vertical motion. Consequently, the light absorbing aerosols aloft serve to strengthen the inversion cap on the boundary layer through warming, and may feed back to cool the surface (Markowicz, 2003), increasing the stability of the troposphere. The upshot is the trapping of more aerosols especially so in basin-type terrain.

Conclusion
The research presented here looked at aerosol-climate interactions in the context of regional climate change. The study involved the interactions between winter aerosols in the IGP and atmospheric stability using several reanalysis datasets and CMIP5 models outputs. The lapse rate analysis across the IGP shows that the atmosphere is becoming more stable, and both GHG's and aerosols have a role in this stabilization. Through an empirical model approach and WRF-Chem simulation, it is shown that aerosols/PM in the IGP can increase in the absence of a change in their emission sources -a process caused by increased anthropogenic GHG that was previously not realized. Moreover, a positive feedback exits whereupon aerosol-induced stabilization increases the accumulation of PM. Collectively the consequences of a warming atmosphere is in the modulation of aerosol concentrations compounded by the IGP's basin-like orography. In the face of increasing anthropogenic GHG's, aerosol loading in the region will worsen whereupon exacerbation of pollution and winter smog should be a future expectation. Figure S2. Scatter plot of changes in PM (MERRAero) and lapse rate (ERA-I) at different pressure levels for November and December months, from 2002 to 2013. Both potential and actual temperature were used, denoted by P and A on the legend, respectively. The long-term trend in each dataset was removed. In parentheses are the correlation coefficients (R), between PM and each lapse rate proxy. The dotted lines in corresponding colors are the least squares best fit to each set, with the optimal level (925-800 hPa) emboldened. Figure S3. Seasonal potential temperature (MERRA) slope at the 800 hPa level . The Tibetan plateau and altitudes above1000 m have been masked out. Units have been converted to per 34 years. Figure S4: Slope of precipitation (UDEL) and 2-m temperature (ERA-I), for different time periods. The Tibetan plateau and altitudes above 1500 m have been masked out. Units have been converted to per 34 years  and per 20 years (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). Table S1: CMIP5 specifics as used in the study. The CMIP5 models we used that utilize interactive aerosol and prognostic cloud microphysics are emboldened.