Can the weakening of Indian Monsoon be attributed to anthropogenic aerosols?

Literature found that the increasing concentration of Anthropogenic Aerosol (AA) is the key reason behind the weakening trend of Indian Summer Monsoon Rainfall (ISMR), based on the Coupled Model Intercomparison Project Phase5 (CMIP5) simulations with AA-only forcing. Here, we re-examine and find that AA-only simulations show country-wide drying, in contrast to the observed east-west asymmetry in the recent ISMR trend. For further evaluation, we decompose the changes in moisture convergence during summer monsoon into dynamic and thermodynamic components. We find that multi-run ensemble averages for individual CMIP5 models do not capture the observed dominance of the changes in dynamic component over the thermodynamic one. An optimal fingerprinting technique for detection and attribution also fail to attribute the changes in ISMR to AA, either because of large internal variability and/or intermodel spread. This implies the need for more careful assessment of AA-only simulations for the ISMR before attributing the changes to AA.


Introduction
The spatial and temporal variability of Indian Summer Monsoon Rainfall (ISMR) impacts the water availability, agricultural productivity and the Gross Domestic Product (GDP) of the country. Therefore, a reliable projection of ISMR with its spatio-temporal variability in global warming scenarios is important for planning adaptation and mitigation strategies. Currently the climate models participating in the Coupled Model Intercomparison Project Phase5 (CMIP5) are the best tools available for long term projection of climate variables. However, the CMIP5 models are known for their biases and limitations in simulating climatology and seasonality (Jourdain et al (2013)), inter-annual and intra-seasonal variability (Sperber et al 2013) and multi-decadal trend (Saha et al 2014) of ISMR. Saha et al (2014) found that the ISMR exhibits a negative trend in the latter half of the 20th century; but the majority of the CMIP5 models fail to capture it. The reason for decreasing trend of ISMR and failure of the climate models are still a subject of active research. The greenhouse gas (GHG)-induced warming of the Indian Ocean, as well as the direct and indirect effects of anthropogenic aerosols have been documented as possible reasons for weakening of Indian monsoon (Bollasina et al 2011, Roxy et al 2015. Saha et al (2014) demonstrate that the CMIP5 models are unable to properly capture important geophysical processes associated with monsoon. Sabeerali et al (2015) identify that improving the convective parameterization scheme in coupled models is necessary to capture the observed characteristics of monsoon. Salzmann et al (2014) and Salzmann and Cherian (2015) point out that the internal climate variability plays an important role in determining the trend of Indian monsoon and probability of a CMIP5 model realisation to capture the trend based on internal variability alone is small. Guo et al (2015) show that the CMIP5 models, that are simulating both direct and indirect effects of aerosol, are able to capture the observed decreasing trend of ISMR better.
Aerosols can alter the spatiotemporal pattern of precipitation by scattering and absorbing solar radiation (direct effect) and by changing cloud properties such as albedo, droplet size etc (indirect effect). The impacts of aerosols on Indian summer monsoon are heavily debated in scientific literature.  proposed the 'Elevated Heat Pump' (EHP) hypothesis, where the direct effects of aerosols are responsible strengthening of Indian monsoon. Studies have shown evidences in favour of the EHP hypothesis using both observational datasets  and model simulations (Wang et al 2009). On the other hand, Ramanathan et al (2005) concluded that aerosols are associated with weakening of meridional thermal contrast, Hadley circulation and Indian monsoon. Using model simulations, Bollasina et al (2011) showed that anthropogenic aerosols (AA) play a key role behind the observed decreasing trend of ISMR in the second half of the 20th century. According to Guo et al (2015), a multi-model mean of CMIP5 AA-only experiments can capture the observed decreasing trend of rainfall, where GHG-only simulations show a trend in the opposite direction. Based on the CMIP5 outcomes, they have attributed the weakening of monsoon to anthropogenic aerosols and suggested a better representation of both direct and indirect effects in CMIP5 models will improve the rainfall projection. To explore this idea further, in this study we have analyzed AA-only, GHG-only and All forcing Historical simulations from 9 CMIP5 models and compare various characteristics of rainfall simulated by them with that of observed gridded rainfall product acquired from the India Meteorological department (IMD). More details about the climate models used in this study are given in table S1. In figure 1(A), we have compared the time evolution of ISMR obtained from AA-only, GHG-only and All forcing Historical CMIP5 simulations, as well as observed rainfall, throughout the 20th century. It should be noted that the observed ISMR from IMD shows an increasing trend in the first half of the century, but a slightly decreasing trend, which is not statistically significant unlike other observed rainfall datasets (Ramesh and Goswami 2014), in the second half. On the other hand, The AA-only trend is decreasing throughout the century; it matches with the observed trend, but only for a specific window of time, not the entire century. The GHG-only trend is in opposite direction to that of AA, and the All forcing trend lies in between.
Even though ISMR, the summer monsoon rainfall averaged over the landmass of India, is a widely used metric for quantifying temporal variability and trend of Indian monsoon, it doesn't capture the spatial heterogeneity associated with it. In figure 1(B), we have shown the time evolution of summer monsoon rainfall over India, as well as two homogenous rainfall zones, Central India and North-East India, throughout the 20th century. Central India, which mostly comprises of the core monsoon zone, shows a decreasing trend of rainfall in the second half of the century. However, the rainfall in North-East India has an increasing trend during the same time period. Even though the area with positive rainfall trend in this region is a small one, the magnitude is much higher than that of Central India region (figure S1 is available online at stacks.iop.org/ERC/1/061006/ mmedia). The highly heterogeneous nature of Indian monsoon implies that different zones can have different trend of rainfall and trends of two zones in opposite direction can compensate each other, leaving a false impression of no significant change. However, majority of the scientific studies do not take the observed eastwest asymmetry of rainfall trend into account, while examining the impact of AA. Some of them only consider the trend of average rainfall over the whole country (Ramanathan et al 2005, Guo et al 2015, while others show a country-wide drying effect due to aerosols (Bollasina et al 2011, Sajani et al 2012.
Our goal in this study is to seek answer to this question: can weakening of Indian monsoon in recent decades be attributed to aerosols, based on CMIP5 model outcomes?In order to determine if anthropogenic aerosol forcing can account for the spatial heterogeneity associated with monsoon, we examine how the frequency of different spatial structures of monsoon rainfall changes with time and how the AA-only experiments from CMIP5 models simulate it. We also calculate the dynamic and thermodynamic components of changes in moisture convergence simulated by CMIP5 models for the entire 20th century and compare them with the observed changes. Lastly we use an optimal fingerprinting based attribution-detection technique to figure out if the observed multi-decadal trend of Indian monsoon can be attributed to anthropogenic aerosol.

Study area
The chosen area for this study is the landmass enclosed by the political boundary of India, which is divided into 7 rainfall homogenous zones by IMD (Parthasarathy et al 1996). We have further modified the rainfall zones: introduced a new zone at the Western-Ghat region, and ignored Jammu and Kashmir and North-East Hilly regions for lack of enough rain gauge stations in those zones (figure 1(c)). The Central India and the North-East India are two most notable zones considered for our study, because of their asymmetrical rainfall trends.

Data
Based on the availability of all required datasets, 9 CMIP5 models are chosen for this study, and their details are presented in table S1. 7 out of 9 of these models include parameterizations for indirect effects of aerosol. Multiple ensembles from AA-only, GHG-only and All forcing historical simulations for the following variables: precipitation, zonal and meridional wind speed, and specific humidity, are downloaded at monthly scale from Earth System Grid Federation (ESGF) gateways. Additionally monthly precipitation data for NAT-only and Pre-Industrial Control (CTL) simulations are downloaded from the same source. The long-term  daily gridded observed rainfall product with high spatial resolution (0.25°×0.25°) is acquired from India Meteorological Department (IMD). Monthly meridional and zonal wind speed and specific humidity data from ERA-20C dataset are collected from European Centre for Medium-Range Weather Forecasts (ECMWF) data repository.

Trend detection
Prior to trend detection, we standardize each time-series with respect to its long-term mean and standard deviation, and pass it through an 11-year moving average filter to remove high frequency variability. We employ a modified Mann-Kendall test (Hamed and Rao, 1998) to detect presence of a linear trend at 95% significance level. If a trend is detected, we calculate it by linear regression.

Spatial structure-based classification
We segregate the total Indian summer monsoon rainfall into 6 zonal rainfall time-series, each of which is standardized with respect to its own long-term mean and standard deviation, and categorized into one of 7 classes based on supplementary table S2. In order to identify the years with similar spatial pattern of rainfall on Indian landmass, we apply k-means clustering technique with Euclidean distance to those 7 classes. Each cluster represents a unique spatial pattern of ISMR. We evaluate the clusters with Silhouette index (Rousseeuw, 1987) and find that 9 is a reasonably optimal number of clusters. We apply the clustering methodology on the model simulated data by reusing the centroid obtained from observed data classification. More details on the classification methodology can be found in Sahastrabuddhe et al (2018).

Identification of dynamic and thermodynamic components
We follow the methodology proposed by Seager et al (2010) to decompose the changes in moisture convergence into its dynamic, thermodynamic and transient eddy component.. In this study we use long term mean of monthly mean data and ignore the contribution of transient eddy moisture convergence. We only focus on mean flow convergence and its components (dynamic and thermodynamic). Basically the dynamic component is a result of changes in wind circulation, while the specific humidity is constant. On the other hand, the thermodynamic component is due to changes in specific humidity without any change in wind circulation. More details on the derivation of these components are provided in the supplementary text.

Attribution and detection
We have applied ordinary least square approach for optimal fingerprinting method (Allen and Tett 1999) in order to assess the contribution of external forcings (AA, GHG and NAT) in the observed changes in rainfall. We use linear regression to estimate the regression coefficients (also known as scaling factors, β), that should be applied to the amplitude of the simulated responses (X) to the forcings in order to best fit the observations (y). The internal climate variability (ε) of models is estimated from CTL simulations.
The upper and lower limit of β is estimated from a student's t-distribution, for 95% confidence interval. A signal is considered detected when the associated scaling factor is not consistent with zero. Attribution additionally requires the scaling factor to be consistent with unity. More details on the optimal fingerprinting methodology are provided in the supplementary text.

Results
We divide the whole 20th century into four epochs of 25 years each and combine them to form three overlapping 50-years time periods (1901-1950, 1926-1975 and 1951-2000), for which we calculate rainfall trends simulated by AA-only, GHG-only and All forcing historical simulations from 9 CMIP5 models for the All India, Central India and North-East India region. We compare these calculated trends with observed trends and show them in figure S2. We find that the simulated trend by models may reasonably match the observed trend at certain time period, region and for certain forcing, but in majority of the cases the models fail to capture it.
To examine the spatial structures associated with Indian monsoon rainfall trend (figure S1), we classify each year of the 20th century into 9 clusters, based on the rainfall anomaly in 6 homogenous rainfall zones in India. We apply the location of the centroid of the clusters obtained from the observed data, to classify the model simulated data as well. In figure S3, we show the spatial pattern of June-September (JJAS) precipitation associated with each cluster for AA-only, GHG-only and All forcing historical simulations from CMIP5, as well as the observed dataset from IMD. The cluster 1 is the most notable one, as it captures a negative anomaly in the Central India rainfall and a positive anomaly in rainfall at north-eastern region. It is expected that an increase in cluster 1 frequency is associated with the observed rainfall trend at these regions. In figure 2(a), we have plotted the frequency of each cluster as simulated by the AA-only CMIP5 simulations in each epoch and compared them with the observed frequencies. We find that indeed the number of occurrence of cluster 1 has increased highly between epoch 3 and 4. But in the AA-only model simulated rainfall the frequency of cluster 1 is decreasing ( figure 2(a)). The main reasons for the weakening of monsoon in AA-only simulations are increase in occurrence of cluster 4, 5 and 7 and decrease in occurrence of cluster 6, which is not the case for observed data except for cluster 5. None of the aforementioned clusters are associated with east-west asymmetry of rainfall anomaly and most of them correspond to either positive or negative anomalous rainfall all over India. We further classify each cluster into three categories: dry, normal and wet, for each homogenous zone, as well as all India; and depending on the zone being considered, a cluster can shift from one category to another (table S3). In figure 2(b) we show the frequency of dry, normal and wet years, as simulated by AA-only CMIP5 simulations for every epoch, for Central India, North-East India and All India region and compared them with observed frequencies. The slope of the line joining two epochs is an indicator of the trend of rainfall in the corresponding 50-years period. As expected, we find an increase in wet years and decrease in dry years in North-East region, and the opposite in the central India between epoch 3 and 4 in the observed data. At the all India level, the number of dry years is constant, and there is a slight decrease in the number of wet years. The AA-only CMIP5 simulations are unable to capture these changes regionally. We have also shown the frequencies of each cluster and each category for GHG-only and All forcing simulations in figure S4 and figure S5 respectively. To ensure that the results obtained from this analysis are not heavily dependent on the chosen time period, we performed clustering with the same dataset, but only with more recent years . A comparison of the outcomes of these two analyses is shown in the table S4. We find that the majority of the overlapped years do belong to the same cluster, which suggests robustness of the result. Overall, the AA-only simulations have a tendency to increase the number of dry years and decrease the number of wet years throughout the century. The GHG-only simulations show a reverse effect and the All forcing simulations correspond to somewhere in between. None of the forcing scenarios match the observed changes for the whole century.
The failure of CMIP5 models to capture the change of frequencies of spatial structure associated with Indian monsoon rainfall suggests that they may not be simulating the underlying geophysical processes, which are affecting those changes, properly. Saha et al (2014) investigated some of those geophysical processes, such as southern Indian Ocean warming and anomalous cyclonic formation at tropical western Pacific, and found that CMIP5 models have limitations in capturing these processes. Ultimately, inability to capture these processes will be reflected in the moisture convergence change simulated by the models. In this study we examine the observed change in moisture convergence over Indian region and how CMIP5 models are simulating it. Changes in the long term mean of moisture convergence are driven by changes in moisture holding capacity of the atmosphere due to warming (thermodynamic component) and changes in the wind circulation (dynamic component). We calculate the observed and model simulated dynamic component of moisture convergence change, by only considering the changes in wind circulation due to warming, but keeping the specific humidity constant to its historical level; and vice versa for the thermodynamic component. We have shown the observed changes in moisture convergence between epoch 3 (1951-75) and epoch 4 (1976-2000), as captured by the ERA20C dataset in figure 3(A). The foothills of Himalaya and the north-eastern India receive an increased moisture convergence in later decades, while Western and Central India, as well as Arabian Sea and Bay of Bengal receive decreased moisture convergence. We decompose the moisture convergence into its components and show it in the figures 3(B) and (C). The dynamic component dominates and closely matches the spatial pattern of the change in moisture convergence. In case of the thermodynamic component, the western-ghat region and some part of central India shows a slight decrease in moisture convergence, and the rest of the country shows slight increase; overall the magnitude of these changes are much less compared to that of dynamic component. We also show the contribution of precipitation and evapotranspiration change to moisture convergence change separately in figure S6 and find that change in precipitation dominates over that of evapotranspiration. The dominant dynamic component of moisture convergence is correlated with the precipitation change. The spatial pattern of the thermodynamic component has some resemblance with the evaporation change, but not everywhere. We show the thermodynamic and dynamic component of moisture convergence change simulated by 9 CMIP5 models with AA-only forcing in figures 3(D) and (E) respectively. The ensemble-averaged model simulated dynamic and thermodynamic components neither match the ERA20C outcomes, nor consistent with each other, in terms of spatial pattern. The magnitudes of both the components in majority of the models are comparable, and don't reflect the dominance of dynamic component as shown in observed changes. We also compare the GHG-only and All forcing model simulated convergence and their components and in none of the cases CMIP5 models have a consensus among themselves. Performing the analysis for other time periods give a consistent result.
However, analysing the spatial pattern of climate changes solely based on ensemble means of multiple model realisations can be misleading. As Salzmann et al (2014) and Salzmann and Cherian (2015) pointed out, internal climate variability plays a significant role in determining the Indian monsoon trend along with aerosol and GHG forcings. If the observation is considered to be one possible realisation of a non-linear system, among multiple, it cannot be expected to match the ensemble mean of model realisations. We repeat the analysis for each realization separately and find that the spatial patterns simulated by them largely vary from each other due to high internal variability. For example we presented the dynamic and thermodynamic components of moisture convergence simulated by all realisations from GISS E2 H model's AA experiment in figure S7. The magnitudes of the simulated changes from each realization are much higher than that of the ensemble mean, because some of them capture changes in opposite direction and compensate each other. We didn't find a single realisation among all the models and forcing experiments that matches both the dynamic and thermodynamic components of moisture convergance change reasonably. However, some realisations in figure S7 (r3i1p107 and r5i1p107) are closer to observations than others.
Even though the CMIP5 models are clearly unable to capture the magnitude and spatial pattern of change in rainfall and moisture convergence, it may still be possible to statistically detect the responses to different forcings simulated by them into the observed precipitation. We apply an ordinary least square (OLS) based optimal fingerprinting technique to determine if AA, GHG or natural (NAT) forcing can be attributed to the changes in the ISMR in 20th century. We present the calculated scaling factors with their limits in table S5 and plot them in figure 4. To attribute the changes to a certain forcing, the corresponding scaling factor with its upper and lower bound, needs to intersect with 1 and not overlap with 0 at the same time. However none of the forcings match these criteria. Best estimates of the scaling factors are smaller than unity, which indicates that the responses to the forcings are overestimated by the models. Wide limits of estimated β suggest high internal climate variability in the models and/or high inter-model spread.

Conclusions
In this study we examine multiple single forcing experiments from CMIP5 model simulations in order to find if the weakening of Indian monsoon in recent decades can be attributed to any of those forcings. There is some uncertainty regarding the trend of Indian monsoon, as multiple sources of observations differ with each other (Ramesh and Goswami 2014). Apart from the IMD gridded dataset, we also examined the trend using 2 other observation dataset (APHRODITE and IITM) for 1951-2005 and presented a comparison in figure S8. The time evolution of ISMR and Central India rainfall closely matches among different dataset. However, the North-East India rainfall captured by IMD differs significantly from the other two datasets. Overall the positive trend of North-East Inda rainfall is also captured by other datasets, but with less magnitude than that of IMD, and sometimes not significant at 95% significance level. A similar trend is also present in the precipitation captured by ERA20C reanalysis dataset ( figure S6(b)). In case of CMIP5 simulations, there is a general tendency of AAonly simulations to decrease the strength of monsoon, while the GHG-only simulations increase the strength over time. Considering the similarity between AA-only simulated and observed ISMR trend it could be hypothesized that the anthropogenic aerosols are responsible for the observed weakening. However analyzing . Scaling factors (β) for AA, GHG and NAT forcings, estimated by OLS-based approach for optimal fingerprinting technique of attribution-detection. The dots represent the mean β and the errorbars represent the limits, assuming a student's t-distribution with 95% confidence interval. the spatial and temporal pattern of the changes in Indian monsoon rainfall reveals that no particular forcing can explain the observed pattern of rainfall at every rainfall zone and at all time period throughout the century. For example, a multi-model mean of changes between 1951-1975 and 1976-2000 rainfall captured by CMIP5 AA simulations show a mild positive trend all over the country, which is far from the observed scenario ( figure S9). The similarity between model-simulated and observed trend for a specific time period or region is possibly a result of internal climate varibility. The weakening or strengthening of monsoon associated with internal climate variability can align with different time periods in different realisations of model simulations, which makes attribution of the observed changes to a particular forcing difficult (Salzmann and Cherian (2015)). Apart from the uncertainty associated with internal climate variability, the impacts of AA and GHG forcings on Indian monsoon trend are also not conclusive. Even though it is generally considered that AA is associated with decreasing ISMR and GHG is associated with increasing ISMR based on CMIP5 model simulations, elevated heat pump hypothesis ) and studies based on Indian Ocean warming differ with that (Roxy et al 2015). We fail to detect the signals from the responses of AA, GHG and NAT forcings by employing an OLSbased optimal fingerprinting technique for attribution-detection, either because of large internal variability and/or intermodel spread. Examining the frequencies of spatial structures associated with monsoon rainfall and the dynamic and thermodynamic components of moisture convergence change shows that ensemble averages for individual CMIP5 models are unable to capture the underlying geophysical processes. Our findings suggest that even though the AA-only simulations of CMIP5 models are capturing the trend of weakening Indian monsoon, the observed changes in monsoon cannot be attributed to AA based on CMIP5 model outcomes alone.