Long-term observations and modeling of aerosol loading over the Indo-Gangetic plains (IGP), the Indian desert region and Himalayan slopes are analyzed in the present study. The Box-Jenkins popular ARIMA (AutoRegressive Integrated Moving Average) model was applied to simulate the monthly-mean Terra MODIS (MODerate Resolution Imaging Spectroradiometer) Aerosol Optical Depth (AOD550 nm) over eight sites in the region covering a period of about 13years (March 2000-May 2012). The autocorrelation structure has been analyzed indicating a deterministic pattern in the time series that it regains its structure every 24month period. The ARIMA models namely ARIMA(2,0,12), ARIMA(1,0,6), ARIMA(3,0,0), ARIMA(2,0,13) ARIMA(0,0,12), ARIMA(2,0,2), ARIMA(1,0,12) and ARIMA(0,0,1) have been developed as the most suitable for simulating and forecasting the monthly-mean AOD over the eight selected locations. The Stationary R-squared, R-squared, Root Mean Square Error (RMSE) and Normalized BIC (Bayesian Information Criterion) are used to test the validity and applicability of the developed ARIMA models revealing adequate accuracy in the model performance. The values of Hurst Exponent, Fractal Dimension and Predictability Index for AODs are about 0.5, 1.5 and 0, respectively, suggesting that the AODs in all sites follow the Brownian time-series motion (true random walk). High AOD values (>0.7) are observed over the industrialized and densely-populated IGP sites associated with low ones over the foothills/slopes of the Himalayas. The trends in AOD during the ~13-year period differentiates depending on season and site. During post-monsoon and winter accumulation of aerosols and increasing trends are shown over IGP sites, which are neutralized in pre-monsoon and become slightly negative in monsoon. The AOD over the Himalayan sites does not exhibit any significant trend and seems to be practically unaffected by the aerosol built-up over IGP. © 2014 Elsevier B.V.