Assessment of regional aerosol radiative effects under the SWAAMI campaign – Part 1: Quality-enhanced estimation of columnar aerosol extinction and absorption over the Indian subcontinent

Improving the accuracy of regional aerosol climate impact assessment calls for improvement in the accuracy of regional aerosol radiative effect (ARE) estimation. One of the most important means of achieving this is to use spatially homogeneous and temporally continuous datasets of critical aerosol properties, such as spectral aerosol optical depth (AOD) and single scattering albedo (SSA), which are the most important parameters for estimating aerosol radiative effects. However, observations do not provide the above; the space-borne observations though provide wide spatial coverage, are temporal snapshots and suffer from possible sensor degradation over extended periods. On the other hand, the ground-based measurements provide more accurate and temporally continuous data but are spatially near-point observations. Realizing the need for spatially homogeneous and temporally continuous datasets on one hand and the near non-existence of such data over the south Asian region (which is one of the regions where aerosols show large heterogeneity in most of their properties), construction of accurate gridded aerosol products by synthesizing the long-term space-borne and ground-based data has been taken up as an important objective of the South West Asian Aerosol Monsoon Interactions (SWAAMI), a joint Indo-UK field campaign, aiming at characterizing aerosol–monsoon links and their variabilities over the Indian region.

In Part 1 of this two-part paper, we present spatially homogeneous gridded datasets of AOD and absorption aerosol optical depth (AAOD), generated for the first time over this region. These data products are developed by merging the highly accurate aerosol measurements from the dense networks of 44 (for AOD) and 34 (for AAOD) ground-based observatories of Aerosol Radiative Forcing over India NETwork (ARFINET) and AErosol RObotic NETwork (AERONET) spread across the Indian region, with satellite-retrieved AOD and AAOD, following statistical assimilation schemes. The satellite data used for AOD assimilation include AODs retrieved from MODerate Imaging Spectroradiometer (MODIS) and Multiangle Imaging SpectroRadiometer (MISR) over the same domain. For AAOD, the ground-based black carbon (BC) mass concentration measurements from the network of 34 ARFINET observatories and satellite-based (Kalpana-1, INSAT-3A) infrared (IR) radiance measurements are blended with gridded AAODs (500 nm, monthly mean) derived from Ozone Monitoring Instrument (OMI)-retrieved AAODs (at 354 and 388 nm). The details of the assimilation methods and the gridded datasets generated are presented in this paper.

The merged gridded AOD and AAOD products thus generated are validated against the data from independent ground-based observatories, which were not used for the assimilation process but are representative of different subregions of the complex domain. This validation exercise revealed that the independent ground-based measurements are better confirmed by merged datasets than the respective satellite products. As ensured by assimilation techniques employed, the uncertainties in merged AODs and AAODs are significantly less than those in corresponding satellite products. These merged products also all exhibit important large-scale spatial and temporal features which are already reported for this region. Nonetheless, the merged AODs and AAODs are significantly different in magnitude from the respective satellite products. On the background of above-mentioned quality enhancements demonstrated by merged products, we have employed them for deriving the columnar SSA and analysed its spatiotemporal characteristics. The columnar SSA thus derived has demonstrated distinct seasonal variation over various representative subregions of the study domain. The uncertainties in the derived SSA are observed to be substantially less than those in OMI SSA. On the backdrop of these benefits, the merged datasets are employed for the estimation of regional aerosol radiative effects (direct), the results of which would be presented in a companion paper, Part 2 of this two-part paper.

1 Introduction

The climate forcing potential of atmospheric aerosols is well accepted by the global scientific community and policy makers (Boucher et al., 2013). This forcing can affect Earth's hydrological cycle (Ramanathan et al., 2001; Bollasina et al., 2011), increase the stability of the atmosphere (Ramanathan and Carmichael, 2008; Jacobson and Kaufman, 2006; Petäjä et al., 2016) and have a significant impact on the Indian summer monsoon (Lau and Kim, 2006). Along with these climatic impacts, aerosols are shown to have adverse effects on human health (Dockery et al., 1993; Seaton et al., 1995; Pope III et al., 2002). Accurate assessment of these impacts still remains a challenge, primarily due to the inadequate spatiotemporal coverage of the aerosol properties such as the aerosol optical depth (AOD) and single scattering albedo (SSA), and the large uncertainties prevailing in the available database. This is especially so over the Indian region, which is among regions having high aerosol loading that shows large heterogeneities in its spatial and temporal characteristics. These heterogeneities are primarily because of wide diversity in geography, anthropogenic activities and meteorological features at mesoscale and synoptic scale. As demonstrated by the past studies (Haywood and Shine, 1995, 1997; Heintzenberg and Helas, 1997; Russell et al., 2002; Takemura et al., 2002; Loeb and Su, 2010; Babu et al., 2016), a small change in the SSA can even alter the sign of aerosol radiative forcing (at the top of the atmosphere) from positive (warming) to negative (cooling), especially over highly reflecting surfaces. The large spatial heterogeneity in the surface reflectance of the land mass over this region and its seasonality makes the aerosol radiative forcing estimation all the more complex. Given this background, construction of gridded datasets of aerosol parameters, especially AOD and SSA, with reduced uncertainties and fairly homogeneous spatial and temporal distribution over the region, becomes imperative. One way to achieve this is data assimilation, a mathematical technique of generating a dataset with reduced uncertainties by systematically combining multiple datasets (which individually may have higher uncertainties) (Kalnay, 2003; Lewis et al., 2006).

Various space-borne sensors aboard remote sensing satellites (such as MODIS, MISR, OMI, etc.) provide the global datasets for spatial and temporal distributions of AOD and absorption AOD (AAOD) (Kaufman et al., 1997; Diner et al., 1998; Chu et al., 2002; Remer et al., 2005; Torres et al., 2007). Despite their wide spatial coverage, the satellite-retrieved data suffer from substantial biases and uncertainties due to cloud contamination, various assumptions made during the retrieval procedure, large spatial heterogeneity in the ground reflectance (over the heterogeneous landmass) and also due to very little information on variation during a day (due to snapshot nature of measurements). In addition, satellite retrievals suffer from issues regarding sensor calibration (Zhang and Reid, 2006; Jethva et al., 2014). Especially over land with heterogeneous surface reflectance, satellite-retrieved AODs depict higher uncertainties (Jethva et al., 2009). On the other hand, being direct measurements, AOD or black carbon (henceforth BC, which is the primary absorbing aerosol species) mass concentrations measured, respectively, using ground-based, periodically calibrated Sun photometers and aethalometers are quite accurate and have large temporal coverage in a day as well as over the years (Moorthy et al., 1989; Holben et al., 1998; Hansen and Novakov, 1990; Babu et al., 2004), along with smaller uncertainties than their satellite counterparts. However, their limited spatial representativeness (more like point measurements) calls for a dense network of observations for a reasonable spatial coverage; even then, remote and inaccessible areas remain undersampled. Moreover, practical constraints result in spatially non-uniform distribution of the ground-based stations.

These limitations of satellite-retrieved (SR) and ground-measured (GR) aerosol parameters restrict their applicability for climate impact assessment studies over heterogeneous regions, like the vast Indian region. However, the relative advantages of these two datasets could be effectively employed for improving regional radiative forcing estimation if these different independent datasets could be assimilated to generate a more accurate and spatiotemporally continuous gridded dataset following established statistical assimilation techniques.

There have been a few efforts in the past to combine AODs from various sources, regionally and globally. Collins et al. (2001) have assimilated AODs retrieved by the Advanced Very High Resolution Radiometer (AVHRR) with those simulated by the Multi-scale Atmospheric Transport and CHemistry (MATCH) model for generating forecasts of aerosols during INDian Ocean EXperiemt (INDOEX). Since then, a few more studies have focused on assimilating satellite-retrieved aerosol products with those simulated by regional/global chemistry transport models (Yu et al., 2003; Generoso et al., 2007; Niu et al., 2008; Zhang et al., 2008). Benedetti et al. (2009) have incorporated AOD assimilation as an integral part of the weather forecasting system at the European Centre for Medium-Range Weather Forecasts (ECMWF). However, in all these efforts, the focus was to use satellite products with chemistry transport models, and none of these studies have assimilated AOD observations from the ground-based network of Sun photometers with the corresponding satellite-retrieved parameters. Probably the first effort in this direction was by Chung et al. (2005), who have assimilated monthly mean AODs retrieved by MODerate Imaging Spectroradiometer (MODIS) with those simulated by a global chemistry transport model, and the resulting AODs are further integrated with monthly mean AOD measurements from Aerosol RObotic NETwork (AERONET) in order to generate the global merged AOD product. Over the Asian region, Adhikary et al. (2008) have assimilated monthly mean AERONET AODs with monthly averaged MODIS AODs and these combined AODs are further assimilated with those simulated by regional chemistry transport model. Nonetheless, over the Indian region (bounded between 0.5–34.5^∘ N and 65.5–96.5^∘ E; Fig. 1), both of these studies have employed ground-based measurements from just two AERONET stations (Kanpur: 26.51^∘ N, 80.23^∘ E; Hanimaadhoo: 6.74^∘ N, 73.17^∘ E). Due to this, the final assimilated AODs (for Chung et al., 2005 and Adhikary et al., 2008) over most parts of the Indian region are largely represented by satellite-retrieved AODs with their inherent large uncertainties as discussed earlier. More recently, Singh et al. (2017) have combined AODs simulated by ECMWF with those retrieved by MODIS and Multiangle Imaging SpectroRadiometer (MISR) as well as in situ measured AODs by a total of 35 AERONET stations spread over the Indian as well as Arabian regions. However, even in this case, employing about 17 AERONET stations over the Indian region, most of these stations were in the monsoon trough region and north-east India with no representation of other parts of the domain. The situation is still worse for SSA.

Thus, developing spatially and temporally continuous gridded products for AOD and SSA using long-term measurements from the dense network of ground-based aerosol observatories (covering most parts of the Indian region) and satellite-retrieved products still remained a dire necessity. This was recognized as one of the most important objectives of the South West Asian Aerosol Monsoon Interactions (SWAAMI) (https://gtr.ukri.org/projects?ref=NE%2FL013886%2F1, last access: 13 February 2019) (Morgan et al., 2016), a coordinated field campaign undertaken jointly by the Indian and UK scientists, and formed its important package.

Accordingly, we have used long-term (2001–2013) measurements of AOD at 550 nm from the two widely used space-borne sensors, MODIS and MISR, over the Indian region and the accurate, quality-checked AOD from a network of 44 ground-based Sun photometers – Aerosol Radiative Forcing over India NETwork (ARFINET) and AERONET – for the same period to generate a gridded dataset for AOD using a modified form of a well-established data assimilation technique. On the similar lines, we have also generated a spatially homogeneous gridded product for AAOD by combining the AAODs estimated using ground-based BC measurements and space-borne infrared radiance measurements (to delineate the dust contribution to AAOD) with AAODs (500 nm) derived from Ozone Monitoring Instrument (OMI) retrievals. These merged datasets for AAOD and AAOD are further employed to estimate columnar SSA at 1^∘ × 1^∘ over the domain.

In Part 1 of this two-part paper, we provide the details of datasets employed for merging and the assimilation methodologies used (for the merging process), in Sects. 2 and 3, respectively. The validation of merged AODs and AAODs against independent ground-based measurements, which did not take part in assimilation, is presented in Sect. 4.1. The merged datasets are then used to examine the spatial distribution of AOD and AAOD over the Indian region (Sect. 4.2) and to delineate their seasonality over the spatially homogeneous subregions of the study domain (Sect. 4.2 and 4.3). Further, using the validated, merged datasets, the gridded product for columnar SSA is derived and the regional-scale as well as subregional-scale SSA characteristics are presented in Sect. 4.3.

2 Database

2.1 Satellite-retrieved (SR) AOD

Monthly mean AOD (at 550 nm) products from MODIS aboard the Aqua and Terra satellites (Kaufman et al., 1997) (L3, Collection 6; https://modis.gsfc.nasa.gov/data/, last access: 13 February 2019), as well as from MISR aboard the Terra satellite (Diner et al., 1998) (L3; https://misr.jpl.nasa.gov/, last access: 13 February 2019), are used as the background data in this study. The MODIS AOD product having spatial resolution of 1^∘ × 1^∘ is constructed by merging AODs retrieved with enhanced Deep Blue algorithm (Hsu et al., 2013; Sayer et al., 2013) and Dark Target algorithm (Levy et al., 2013), in order to provide AODs over bright surfaces (deserts, arid regions, semiarid regions, etc.) as well as oceans. The MISR AOD product at 555 nm, with spatial resolution of 0.5^∘ × 0.5^∘, is regridded to 1^∘ × 1^∘ resolution for combining with MODIS AODs for minimizing data gaps in the background AODs due to non-availability of MODIS data. Being derived from finer-resolution measurements (1 km), the MODIS AOD product is preferred over MISR AOD while constructing background data. In the present study, AODs from MODIS Terra formed a first layer of background data with gaps in it being filled by AODs from MODIS Aqua, if present. Any data gaps existing further were filled with the regridded MISR AODs. Henceforth, the term “SR AOD” will refer to this integrated satellite-retrieved AOD (MODIS plus MISR).

2.2 Satellite-retrieved absorption AOD

This is obtained from OMI aboard the Aura satellite, which measures the upwelling radiations in the wavelength range of 270–500 nm, at the top of the atmosphere (Levelt et al., 2006). UV aerosol index, AOD and AAOD at 354 and 388 nm are then derived by incorporating the measured backscattered radiation into the inversion algorithm (OMAERUV), which makes use of precomputed reflectance by a set of aerosol models, as detailed by Levelt et al. (2006) and Torres et al. (2005). The AOD and AAODs are then extrapolated to 500 nm by considering the wavelength dependence of the respective retrievals, as specified in the corresponding aerosol models (Torres et al., 2005). In the present study, we have employed the monthly mean level-3 AAODs (500 nm) as the background data (for constructing the merged AAOD product), which is hereafter referred to as SR AAOD.

2.3 AOD from the ground-based Sun photometer network (GR AOD)

Ground-based measurements of AOD used in this study are obtained from the ARFINET observatories (Moorthy et al., 2013; Babu et al., 2013) established by Indian Space Research Organization (ISRO) as well as from the AERONET observatories established (over the Indian region) and operated by different institutions jointly with NASA (Holben et al., 1998). The locations of these 44 observatories, the AOD values from which are used in this study, are shown in Fig. 1.

Figure 1 Locations of the ground-based stations, AODs (a) and BC (surface-level) mass concentration measurements (b) which are used in the present study. In panel (a), blue and red dots represent, respectively, the ARFINET and the AERONET stations, the AODs from which are assimilated, while the pink dots represent the stations, the data from which are used for independent validation (and not used in the assimilation). In panel (b), blue and pink dots denote the ARFINET stations providing BC data for assimilation and validation, respectively.

The ARFINET observatories are set-up as a part of the Aerosol Radiative Forcing over India (ARFI) project of ISRO for investigating the spatial–temporal heterogeneities of aerosols, their spectral characteristics and size distributions, as well as to assess the impact of aerosols on regional radiative forcing. In order to obtain columnar AODs, each observatory in ARFINET is equipped with either a 10-channel multi-wavelength radiometer (MWR) (Moorthy et al., 1989, 2013) and/or Microtops Sun photometer (MSP) (Morys et al., 2001). The details of analysis of these data and intercomparison with other commercial and research-level Sun photometers are available in the literature (Shaw et al., 1973; Moorthy et al., 2013; Kompalli et al., 2010). The present study utilizes the monthly mean AOD data at 500 nm, measured at ARFINET stations which are detailed in Table S1 in the Supplement. The AERONET stations (Table S2 in the Supplement) are comprised of a network of automatic, Sun sky-scanning radiometers, set up and supervised by NASA (Holben et al., 1998, 2001). The present study employs level-2 monthly mean AODs at 500 nm provided by AERONET (https://aeronet.gsfc.nasa.gov/, last access: 13 February 2019).

2.4 AAOD from ground-based BC measurements

Unlike AOD, there are no direct ground-based measurements for absorption AOD. Hence, in order to construct a reliable dataset of AAOD over the Indian region, we have employed the regular BC mass concentration measurements (at surface level) performed at 34 ARFINET observatories (Fig. 1b). The list of these observatories along with their geographical coordinates and broad geographical features of the respective regions is provided in Table S3 (the Supplement). The reason behind using the BC measurements for the current purpose lies in the fact that black carbon is the primary light-absorbing aerosol species not only because of its ability to absorb the radiations over a wide wavelength range but also due to its longer atmospheric residence time (on the order of a few days to weeks) in the lower troposphere (Babu and Moorthy, 2002). In addition, BC can alter the properties of other aerosol species by mixing with them (Jacobson, 2001). The other strong absorbing aerosol species over this region is the mineral dust, which is perennially present, especially over the north-western arid regions and Indo-Gangetic Plain (IGP).

The continuous measurements for BC mass concentration are performed using the aethalometer (from Magee Scientific Inc., USA) (Hansen et al., 1984; Hansen and Novakov, 1990) at these 34 ARFINET stations (Fig. 1b). For maintaining consistency in measurements and ensuring data quality, these aethalometers are operated under a common protocol and are periodically intercompared.

The measured BC mass concentrations are used to estimate the AOD due to BC, making use of the Optical Properties of Aerosols and Clouds (OPAC) model (Hess et al., 1998). Presently, there is no empirical model available for the vertical distribution of BC over the Indian region. As such, in order to specify the vertical distribution of BC in OPAC, we have considered the representative vertical distribution of BC, based on commonly observed features of vertical heterogeneities of aerosols reported over the Indian region (Satheesh et al., 1999, 2008; Babu et al., 2011). Due to vertical mixing caused by the eddies within the daytime convective boundary layer, aerosols can be considered to be near-uniformly distributed within the planetary boundary layer (PBL). Aircraft and balloon measurements of BC over different regions of India and during different seasons (Satheesh et al., 2008; Suresh Babu et al., 2010; Babu et al., 2011) have shown that during daytime the near-uniform BC mass concentrations are observed until ≈2 km. Accordingly, we have considered the uniform mass concentration of BC within the PBL and an exponential decaying above it following the scale heights (seasonally varying) reported by Yu et al. (2010) using Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) observations over the Indian region.

For deriving the BC AODs from BC mass concentration using OPAC, the PBL heights (PBLHs) provided by the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) reanalysis dataset (Gelaro et al., 2017) are used. The MERRA-2 PBLH dataset has been validated by comparing with PBLH derived from radiosonde and GPS radio occultation measurements over the Indian region by Sathyanadh et al. (2017). Depending on the location, the slope of a linear fit between the MERRA-2 derived and measured PBLH varies between 0.75 to 0.93 (Sathyanadh et al., 2017). Nevertheless, the validation exercise performed by Sathyanadh et al. (2017) is representative of a quite limited period (May to September 2011). Therefore, we have validated MERRA-2 PBLHs with those estimated from radiosonde measurements (downloaded from http://weather.uwyo.edu/upperair/sounding.html, last access: 13 February 2019) performed at eight representative locations (Fig. 2), during the period of 2008 to 2018.

Figure 2 Locations of the ground stations, radiosonde measurements from which are used for the purpose of validating PBLH derived by MERRA-2. These subregional representative stations form a subset of ground-based observatories, AOD and BC mass concentration measurements from which are employed for construction of assimilated AOD and AAOD products.

The scatter plots between spatially collocated MERRA-2 PBLH and those derived from radiosonde measurements over the eight locations are presented in Fig. 3. It can be seen from Fig. 3 that MERRA-2 PBLHs are well correlated with those estimated from radiosonde measurements, although the correlation coefficient varies from 0.63 to 0.96, with respect to the location. The details about this validation exercise are provided in Appendix C.

Figure 3 Comparison of spatiotemporally collocated MERRA-2 PBLHs with those derived from radiosonde measurements performed at eight representative locations during 2008 to 2018. The correlation coefficient (R) (significant at 95 % confidence limit) and the equation of linear regression between the two PBLH estimates are provided in each of the figures.

After deriving BC AODs by incorporating BC mass concentration measurements and MERRA-2 PBLH in OPAC (Hess et al., 1998), the corresponding absorption AODs contributed by BC are then estimated by considering BC SSA as 0.22 (Hess et al., 1998).

Mineral dust aerosols form another important species contributing to aerosol absorption, not only the solar but also the outgoing terrestrial radiation (Satheesh et al., 2007; Deepshikha et al., 2006a, b). Mineral dust prevails over the central and northern parts of Indian region, locally produced as well as advected from west Asian and east African regions (Moorthy et al., 2005; Niranjan et al., 2007; Beegum et al., 2008). As a first step to estimate the dust absorption optical depth, we have computed the infrared difference dust index (IDDI) (Legrand et al., 2001), the reduction in infrared radiance sensed by the satellite-based instruments, due to atmospheric dust aerosols. Mathematically, IDDI is defined as shown in Eq. (1) following Legrand et al. (2001).

$$\begin{array}{}\text{(1)}& \mathrm{IDDI}=\mathrm{RD}\uparrow -\mathrm{RC}\uparrow \end{array}$$

Here, RD↑ and RC↑, respectively, denote the outgoing long-wave terrestrial radiations in dust-loaded and clear-sky (i.e. no aerosol and no clouds) conditions at the top of the atmosphere (TOA). Thus, IDDI is an indicator of amount of columnar dust loading in the atmosphere (Tanré and Legrand, 1991; Legrand et al., 2001). In this study, we estimated the IDDI from the brightness temperature corresponding to IR radiance (10.5 to 12.5 µm) measured by the Very High Resolution Radiometer (VHRR) aboard geostationary Indian satellites, Kalpana-1 and INSAT-3A, following Deepshikha et al. (2006a, b) and Srivastava et al. (2011). The daily estimates of IDDI are then used to infer the dust AOD (500 nm) (following Srivastava et al., 2011) from which dust absorption AODs are estimated by considering dust SSA (at 500 nm) as 0.91 (Moorthy et al., 2007). The monthly mean dust AAODs are then interpolated to the locations of ARFINET laboratories from which the BC measurements are employed. The final AAODs to be merged with the OMI AAODs are then constructed as shown in Eq. (2) and are henceforth referred to as GR AAODs.

$$\begin{array}{}\text{(2)}& \mathrm{GR}\mathrm{AAOD}=\mathrm{BC}\mathrm{AAOD}+\mathrm{dust}\mathrm{AAOD}\end{array}$$

3 Merging the different datasets: assimilation methodology

In generating merged datasets for AOD and AAOD by assimilating different datasets, we have adopted widely accepted statistical data assimilation techniques, along with a few physical constraints, as detailed below.

Several methods are available in the literature for combining scattered observations with gridded data, including variance-minimization-based methods (like 3D-Var) and heuristic methods like successive correction methods (SCMs) (Kalnay, 2003; Lewis et al., 2006). Methods based on variance minimization are mathematically sophisticated and perform the assimilation in such a way that the variance of the entire assimilated field is guaranteed to be lesser than those of the parent datasets. On the other hand, SCMs are empirical in nature and blend in the observations with background data locally within a specified radius of influence, although it does not assure the variance minimization.

As such, Mitra et al. (2003, 2009) have employed Cressman method (Cressman, 1959), a variant of an SCM, for forming daily analyses of rainfall over the Indian region by merging rain gauge measurements with satellite-retrieved precipitation. For AOD assimilation, the similar method has been employed by Chung et al. (2005), while optimal interpolation (OI) has been used by Collins et al. (2001), Adhikary et al. (2008) and Singh et al. (2017). Further, 3D-Var has been employed by Niu et al. (2008) and Zhang et al. (2008) for assimilation of dust aerosol properties and AODs, respectively. Nevertheless, analysed AODs produced by either variance minimization methods or SCM-based methods need not be always bounded by their parents.

In the present study, the datasets to be merged represent the same parameter although measured by different techniques. The ground-based measurements are comprised of data from a dense network of observatories representing all distinct environments over the region and has sufficient temporal coverage to smooth out any isolated events or episodes. Ground-based measurements also provide a fairly accurate spatiotemporal distribution, while the satellite products provide wide spatial coverage. In view of this, it is logical that merged datasets be bounded by the respective parent datasets. This also leads to a fairly smooth spatial variation that is needed for inputting the merged products to climate impact assessment models, without compromising on the accuracy of distinct spatial features.

3.1 Merging methodology for AOD

In order to achieve the variance minimization while ensuring the merged AODs to be bounded by the parent datasets (i.e. SR and GR AOD), we have employed SCM, with a variation, which we refer to as the weighted interpolation method (WIM). It expresses the merged AODs as the weighted average of SR and GR AODs, in such a way that the resulting merged AOD values are always bounded by the parents. While performing the weighted average, the weight given to GR AOD is inversely proportional to the distance between the location corresponding to a given grid point and a ground-based observatory, following the Cressman (1959). This inverse-distance-weighting method enables the merged AOD at a given location to be largely represented by the ground-based measurements from nearer stations with regard to farther ones. Depending on the weights given to GR AOD, WIM assigns weights for SR AOD such that the sum of weights for SR and GR AODs is always unity. This ensures that merged AODs are bounded by ground-based and satellite-retrieved AODs. The weighted average of SR and GR AODs is then performed in an iterative manner until the merged AODs interpolated at the locations of ground-based observatories match with respective GR AODs within their uncertainty limits. Mathematically, WIM is expressed as shown in Eq. (3).

$$\begin{array}{}\text{(3)}& {\mathit{X}}_{k+\mathrm{1}}=\mathbf{R}\left[{\mathit{X}}_k\right]+\mathbf{QW}\left[Z\right]\end{array}$$

Here, X_k+1 = vector (size n×1) of merged data at the k+1th iteration, where n represents total number of grid points in the domain. Mathematically, $n=n_x\times n_y$, where n_x and n_y denote number of nodes in longitude and latitude, respectively. Further, X_k represents vector (size n×1) of merged data at the kth iteration. During the first iteration of Eq. (3), X_k is equal to the vector of background data as provided by SR AOD. Z refers to the vector (size m×1) of GR AOD, where m represents number of ground-based measurements available at that instant in the whole spatial domain. H is an interpolation matrix (size m×n) which bi-linearly interpolates the gridded satellite data to the locations of ground-based observatories. The details about construction of H can be found in Kalnay (2003) and Lewis et al. (2006).

As the satellite-retrieved and ground-based AODs are not collocated, one needs to give appropriate weights to SR and GR AOD values during merging them. In Eq. (3), QW (size n×m) is the normalized weight matrix for GR AODs, and R (size n×n) is the weight matrix for SR AODs. The normalized weight matrix is constructed as a product of two matrices; Q, the normalization matrix (size n×n) and W, the weight matrix (size n×m) for GR AODs. The weights given to GR AODs, which form elements of matrix W, are computed using one of the widely accepted inverse-distance methods which is given by Cressman (1959) as shown in Eq. (4).

$$\begin{array}{}\text{(4)}& {\mathbf{W}}_{ij}=\left\{\begin{array}{ll}\frac{d^{\mathrm{2}}-r_{ij}^{\mathrm{2}}}{d^{\mathrm{2}}+r_{ij}^{\mathrm{2}}},& \mathrm{if}{r}_{ij}\le d\\ \mathrm{0},& \mathrm{otherwise}\end{array}\right.\end{array}$$

Here, W_ij denotes the weight given to GR AOD from jth ground-based observation location during merging with SR AOD at the ith grid point. This weighting strategy (Eq. 4) ensures that the contribution of ground-based measurements to merged AOD is higher (lesser) if the distance between a ground station and a grid point (referred to as r_ij in Eq. 4) is lower (higher). The weight matrix W also makes sure that merged AODs are not contributed by the GR AODs from the ground-based stations lying outside the radius of influence which is denoted by d in Eq. (4). It is to be noted here that the radius of influence corresponds to the region surrounding a given location, within which an AOD from that location can be considered to be largely representative. The details about the choice of radius of influence for the present study are described in the next section describing modified weight matrix formulation.

The weight matrix (W) thus computed (Eq. 4) needs to be normalized to make sure that the sum of the weights given to GR AODs from all the stations is less than unity. This normalizing process thus constrains the merged AODs by the available ground-based measurements. In order to perform this normalization, the diagonal matrix Q (which multiplies to W as shown in Eq. 3) is constructed as follows (Cressman, 1959; Kalnay, 2003; Lewis et al., 2006).

$$\begin{array}{}\text{(5)}& {\mathbf{Q}}_{ii}^{-\mathrm{1}}=\left[\sum _{j=\mathrm{1}}^m{\mathbf{W}}_{ij}+{\displaystyle \frac{{\mathit{\sigma }}_o^{\mathrm{2}}}{{\mathit{\sigma }}_B^{\mathrm{2}}}}\right]\end{array}$$

The first term on the right-hand side of Eq. (5) is the summation over the weights given to all GR AODs from all the stations within the radius of influence from the ith grid point, and the second term is the ratio of error variances in GR AOD (${\mathit{\sigma }}_o^{\mathrm{2}}$) to that in SR AODs (${\mathit{\sigma }}_B^{\mathrm{2}}$). This normalization strategy also ensures that the weights for a parent dataset reduce with the increase in its uncertainty. As a major portion of background data is formed by MODIS AODs, σ_B is represented as rms (root mean square) uncertainty in MODIS AODs, which is given as 0.03+0.2τ_sat by Sayer et al. (2013), where τ_sat is the MODIS AOD. The σ_o term is formed by uncertainty of GR AOD measurements made at ARFINET and AERONET observatories. As the uncertainties in GR AODs at different wavelengths are in the range of 0.01 to 0.03 (Holben et al., 1998; Babu et al., 2013), the maximum uncertainty (i.e. 0.03) is considered as σ_o .

After computing the normalized weights for GR AOD, we have calculated the weights to the SR AODs such that the sum of the weights for SR and GR AODs is unity. In other words, the weights for SR AODs are computed such that the merged AODs are guaranteed to be a convex combination of the parent datasets. Mathematically,

$$\begin{array}{}\text{(6)}& {\mathbf{R}}_{ii}={\mathbf{I}}_{ii}-\sum _{j=\mathrm{1}}^m\left[{\mathbf{Q}}_{ii}{\mathbf{W}}_{ij}\right].\end{array}$$

Here, R is the diagonal matrix with its ith diagonal element referring to the weight given to SR AOD at grid point at the ith grid point. I is the identity matrix of size n×n.

It is also to be mentioned here that WIM (Eq. 3) makes sure that uncertainties in merged (MG) AODs are either less than uncertainties in SR and GR AODs or at least less than the largest of the two. The theoretical proof for this is given in Appendix A.

Modified weight matrix formulation

The weight matrix formulation (Eq. 2) involves the distance (r_ij) between a grid point and a network observatory as well as the radius of influence (d) from that grid point. However, both of these quantities (r_ij and d) can be computed/estimated based on either only horizontal or horizontal and vertical coordinates of the corresponding grid point and ground station, depending upon the nature of problem. In the current problem of AOD merging, apart from horizontal distance, it is essential to consider the altitude difference between a grid point and a ground station. This is because AODs measured at an aerosol observatory located over a sharp peak situated over a large plain terrain may not be representative of AOD corresponding to adjoining grid points over the plains due to sharper variations in aerosol concentrations in the vertical than horizontal direction. Hence, in order to have realistic merging, especially in cases involving merging of AODs from two locations within the horizontal radius of influence yet differing in altitudes significantly, it is necessary to take into account the vertical distribution of aerosols and an associated length scale.

Due to the dynamics of the daytime convective boundary layer and the associated updrafts, aerosols can be considered to be near-uniformly distributed within the PBL; but above this, vertical heterogeneities are possible (for example, Satheesh et al., 2008; Suresh Babu et al., 2010; Babu et al., 2011). As such, the planetary boundary layer height is considered as the region within which the aerosol distribution is near homogeneous in the vertical. However, the vertical gradients in aerosols could be much sharper than the horizontal variations, especially above the PBL. Concentration of aerosols may significantly differ above the top of the PBL, which acts as a virtual lid (although leaky) shielding the free troposphere from surface-based emissions, significantly. A typical example of such a case is the Nainital station located over the mountain peak of nearly 2 km elevation above mean sea level, at a radial distance of < 50 km from the Indo-Gangetic Plain on the south, east and west of it. Similar is the case with a few other stations such as Shillong, Ooty, etc.

As such, we expressed the weight matrix (W) in Eq. (3) as the product of two matrices (scalar product) of the same order, W1 and W2, which take into account the horizontal and vertical variations of aerosols, respectively. The details are as given below.

For the horizontal component (W1), this weight matrix (size n×m) is defined in terms of horizontal radius of influence (d) and the horizontal separation between a grid point and a ground network observatory (r_ij), as given below in Eq. (7).

$$\begin{array}{}\text{(7)}& \mathbf{W}{\mathrm{1}}_{ij}=\left\{\begin{array}{ll}\frac{d^{\mathrm{2}}-r_{ij}^{\mathrm{2}}}{d^{\mathrm{2}}+r_{ij}^{\mathrm{2}}},& \mathrm{if}{r}_{ij}\le d\\ \mathrm{0},& \mathrm{otherwise}\end{array}\right.\end{array}$$

In the present study, the radius of influence is considered to be 250 km (for the first iteration of Eq. 3) following Winker et al. (1996), who have suggested the global horizontal correlation length scale for aerosol to be −200 km, using observations from the Lidar In-space Technology Experiment (LITE). In the current work, the radius of influence is reduced by 50 km during each successive iteration of Eq. (3), to make sure that GR AODs from the location nearest to the given grid point are merged with the background data to the maximum possible extent while iterations converge.

The vertical component of the weight matrix (W2) (size n×m) is configured in terms of height of influence (H) and the altitude difference between a grid point and a network observatory (h_ij) as given below.

$$\begin{array}{}\text{(8)}& \mathbf{W}{\mathrm{2}}_{ij}=\left\{\begin{array}{ll}\frac{H^{\mathrm{2}}-h_{ij}^{\mathrm{2}}}{H^{\mathrm{2}}+h_{ij}^{\mathrm{2}}},& \mathrm{if}{h}_{ij}\le \mathbf{H}\\ \mathrm{0},& \mathrm{otherwise}\end{array}\right.\end{array}$$

Here, h_ij represents the difference between altitudes of the ith grid point and jth observation location. H is the height of influence defined as PBLH + τ , where τ is the height of layer measured above PBL and in which the aerosol concentrations are considered to be decreasing rapidly from the near-constant value within the PBLH which is specified using the MERRA-2 reanalysis dataset. The value of τ has been taken from the variance of PBLH given by MERRA-2. For this, we computed the covariance matrix (size n×n) from monthly mean PBLH over the Indian region. The diagonal elements of this matrix provide variances in the PBLH data at each grid point. After computing standard deviations (σ) from variance values, τ values are taken as 2σ. Based on the relation between H, h_ij and PBLH, any of the three following cases can arise, with each defining the distinctive way in which W1 and W2 contribute to the resultant weight matrix (i.e. W).

1. If h_ij≤PBLH

   In this case, the jth ground station is located within the PBL of the ith grid. In accordance with the consideration of the well-mixed boundary layer, W2 would have the highest weight (=1) and the resultant weight matrix would be solely determined by its horizontal component (Eq. 7). This leads to an element of the resultant weight matrix being expressed as shown in Eq. (9).

   $$\begin{array}{}\text{(9)}& {\mathbf{W}}_{ij}=\mathbf{W}{\mathrm{1}}_{ij}\end{array}$$

2. If $H\ge h_{ij}>\mathrm{PBLH}$

   In this case, the jth ground station is at an altitude just above the PBL of the ith grid but not high enough to be considered uninfluenced by variations at the ith grid point; rather, its influence would be rapidly decreasing. In this case, the resultant weight matrix W is expressed as the scalar product of W1 and W2, so that

   $$\begin{array}{}\text{(10)}& {\mathbf{W}}_{ij}=\mathbf{W}{\mathrm{1}}_{ij}\cdot \mathbf{W}{\mathrm{2}}_{ij}.\end{array}$$

3. If h_ij>PBLH

   In this case, altitude of the jth ground station is high enough above the PBL at the ith grid point such that AOD measured at such a ground station has hardly any relevance to the AOD at the ith grid point. As such, the grid point and ground station are considered to be independent, and W2 is set to zero; hence, the resultant weight matrix, which is the product of W1 and W2, becomes zero (Eq. 11).

   $$\begin{array}{}\text{(11)}& {\mathbf{W}}_{ij}=\mathrm{0}\end{array}$$

After constructing W following the above considerations, normalizing matrix Q is computed by substituting corresponding W into Eq. (5). This is followed by computation of the weight matrix for background data, R, as given in Eq. (6). Following the construction of W, Q and R matrices, Eq. (3) is solved iteratively, with SR AOD being the background data for the first iteration. The solution of the first iteration forms the background data for the second iteration and the procedure is repeated until the norm of the vector of absolute difference between GR AOD and MG AOD interpolated to locations of ground observatories (mathematically, norm(Z−HX_k)) reaches a preset limiting value of 0.02, which is the mean uncertainty in ground-measured AODs (Holben et al., 1998; Babu et al., 2013). In the present study, this condition is satisfied within ≈10 iterations. Nevertheless, in some of the cases, norm(Z−HX_k) gets levelled off before reducing to the limiting error value. This occurs due to AOD measurements at some of the ground stations being unassociated with AODs corresponding to grid points surrounding them. In such cases, iterations of Eq. (3) are performed until the absolute difference between errors (i.e. norm(Z−HX_k) during successive iterations is less than 10⁻³.

3.2 Merging datasets for AAODs

For AAOD merging, the method slightly differed from the above, as the ground-based observations are only of the BC mass concentration measurements which are representative only of surface-level black carbon, unlike AOD which has been columnar for both space-based and ground-based measurements. As detailed in Sect. 2.4, the columnar absorption optical depth for BC is estimated by incorporating the BC number concentration (corresponding to BC mass concentration measurements) into OPAC (Hess et al., 1998) in which the vertical distribution of BC is specified using commonly observed characteristics of vertical heterogeneities of aerosols reported over the Indian region (Satheesh et al., 1999, 2008; Babu et al., 2011). In addition to BC, the contribution of dust is also taken into account to construct the columnar AAODs to be merged with OMI AAODs (Sect. 2.4).

Unlike the ground-based AODs demonstrating much stronger correlation (R=0.77; Fig. 4a) with satellite AODs, the above-mentioned GR AAODs are relatively weakly correlated with OMI AAODs (R=0.35; Fig. 4c). Due to these differences, employing WIM for AAOD assimilation is observed to generate non-smooth and highly discontinuous merging patterns. As such, we have employed one of the widely used data assimilations methods, 3D-Var (Niu et al., 2008; Zhang et al., 2008), which is based on the principle of least-square error minimization. In 3D-Var, the merged AAODs are estimated as a solution of the minimizer of the following objective function (referred to as J) which expresses the weighted sum of the departures in merged AAODs from GR and OMI AAODs, as shown in Eq. (12).

$$\begin{array}{}\text{(12)}& \begin{array}{rl}J\left(X\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\right(\mathit{X}-{\mathit{X}}_b{)}^T{\mathbf{B}}^{-\mathrm{1}}(\mathit{X}-{\mathit{X}}_b)\\ & +(\mathit{Z}-\mathbf{H}\mathit{X}{)}^T{\mathbf{O}}^{-\mathrm{1}}(\mathit{Z}-\mathbf{H}\mathit{X}\left)\right]\end{array}\end{array}$$

Here, X and X_b refer to MG AAOD and OMI AAOD vectors, respectively, of size n×1, where n is the total number of grid points in the spatial domain (n=1120 for the current study). Further, in Eq. (12), Z denotes the vector of GR AAODs, which is of size m×1, where m is the number of ground-based observatories, data from which are available during the respective month. The map between the grid space and the observation space is provided by the interpolation matrix referred to as H (size m×n) (Eq. 12) (Lewis et al., 2006; Kalnay, 2003). Finally, B (size n×n) and O (size m×m) represent the error covariance matrices for OMI AAODs and GR AAODs, respectively. The minimizer to the above-mentioned objective function is estimated by solving Eq. (13).

$$\begin{array}{}\text{(13)}& \left[{\mathbf{B}}^{-\mathrm{1}}+{\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}\mathbf{H}\right]\mathit{X}=\left[{\mathbf{B}}^{-\mathrm{1}}{\mathit{X}}_b+{\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}\mathit{Z}\right]\end{array}$$

Further details about 3D-Var can be found in Kalnay (2003) and Lewis et al. (2006).

Constructing error covariance matrices (B and O) is a fundamental element of 3D-Var data assimilation. This is mainly because the underlying correlation structure and the actual variance values not only dictate the pattern in which observations get merged with the background data but also decide the weights given to each of the parent datasets during the merging process. In the present study, the observation error covariance matrix (O) is considered to be diagonal, implying that errors in GR AAODs from different ground-based stations are uncorrelated, which is generally true and is followed earlier also (Niu et al., 2008; Zhang et al., 2008; Singh et al., 2017). As the diagonal terms of the covariance matrix refer to variance of the corresponding data, the diagonal terms of O are formed by taking the square of uncertainties in the GR AAODs which are estimated as explained below.

It can be understood that the uncertainties in BC AAOD arise largely from the uncertainties in BC mass concentration measurements, the assumed vertical distribution of BC as well as uncertainties associated with the OPAC model. So, in order to estimate the uncertainties in BC AAODs, we perturbed BC mass concentration measurements, MERRA-2 PBLH and scale height within their respective uncertainty limits to compute the multiple realizations for a set of BC AAODs. For this exercise, the uncertainties in BC measurements are considered to be 2 % to 5 % (Hansen and Novakov, 1990; Babu et al., 2004; Dumka et al., 2010) and those in MERRA-2 PBLH are estimated to be 5 % to 20 %, while the uncertainties in scale height for vertical distribution of aerosols (derived from CALIPSO measurements) are considered to be ≈100 m (Kim et al., 2008). The standard deviation of the multiple realizations for a given BC AAOD is adopted as the uncertainty in the corresponding BC AAOD. This analysis showed that the uncertainties in BC AAODs vary from around 11 % to 20 % with its mean, i.e. 15 % being considered as the uncertainty in BC AAOD. Similarly, the uncertainties in dust AAOD, which are largely emanating from the uncertainties in vertical heterogeneity of dust and its optical properties are estimated to be around 25 % of dust AAOD. The diagonal terms of error covariance matrix for observations (O) are thus constructed as shown in Eq. (14).

$$\begin{array}{}\text{(14)}& O_{ii}=(\mathrm{0.15}\cdot \mathrm{BC}{\mathrm{AAOD}}_i{)}^{\mathrm{2}}+(\mathrm{0.25}\cdot \mathrm{dust}{\mathrm{AAOD}}_i{)}^{\mathrm{2}}\end{array}$$

Here, i represents the index varying from 1 to the number of ground-based stations.

As the observation error covariance matrix is diagonal, the patterns of merging GR AAODs with OMI AAODs are fully dependent on the background error covariance matrix (B). In view of this, we have estimated the background error covariance matrix from historical time series (2005 to 2016) of OMI AAOD at 500 nm. This covariance matrix provides the spatial structure of the correlation between OMI AAODs at n grid points as well as the variances of the background data. The details about the construction of seasonally varying B from time series of OMI AAOD are provided in Sect. S2 in the Supplement.

It is to be noted that 3D-Var assures the variance of merged estimates to be lesser than those of both parent datasets (Lewis et al., 2006; Kalnay, 2003). However, the theoretical proof for variance in analysed (assimilated) estimates (constructed by 3D-Var) being smaller than those in parent datasets is provided in Appendix B.

This translates to the uncertainties in the merged AAODs being guaranteed to be smaller than those in OMI AAODs and GR AAODs.

4 Results and discussion

Following the above methods, we constructed spatially and temporally homogeneous gridded data products of AOD and AAOD over the study domain, for each month of the year. Before examining the products for their basic features, it is essential to validate them with independent measurements.

4.1 Validation of merged products

For the validation purpose, we have evaluated the performance of the merged datasets against independent ground-based measurements from subregional representative locations, the data from which did not enter the assimilation process. The merged AOD (AAOD) product constructed by assimilating long-term, ground-based AODs (AAODs) from the network of 36 (26) observatories are validated against the AODs (AAODs) from eight (eight) independent representative observatories, which are shown by pink dots in Fig. 1. As the grid nodes and locations of ground-based observatories are not collocated, we have interpolated the merged AODs and AAODs from the grid nodes contained by the 3^∘ × 3^∘ box surrounding the locations of respective ground locations used for validation.

The comparisons of collocated merged AODs and AAODs with the respective, independent ground-based estimates are shown by scatter plots in Figs. 4b and 5b. To assess the quality improvement due to present assimilation, we have shown the scatter plots of the satellite-retrieved AOD and AAOD (interpolated to ground station locations) against those from the corresponding independent measurements in Figs. 4a and 5a, respectively.

The regression lines and the ideal 1:1 lines are also drawn in the respective panels and the corresponding statistics (regression coefficients and correlation coefficient) are also provided in Figs. 4 and 5.

Figure 4 Validation of merged AOD product: (a) satellite AOD vs. GR AOD; (b) merged AOD vs. GR AOD. Red lines are regression fitted to the point, while the dotted black lines represent the ideal 1:1 case.

Figure 5 Validation of merged AAOD product: (a) satellite AAOD vs. GR AAOD; (b) merged AAOD vs. GR AAOD. Red lines are regression fitted to the point, while the dotted black lines represent the ideal 1:1 case.

The merged products are demonstrating improved agreement (stronger correlation) with independent ground-based datasets than that shown by respective satellite products (Figs. 4, 5). This highlights the significant advantage of the assimilation and is all the more important for AAOD (Fig. 5), the most important parameter for the accurate estimation of atmospheric forcing. With this confidence established through statistical means, we proceeded then to include these independent stations also into the group of ground locations used for merging, and the whole assimilation process (Sect. “Modified weight matrix formulation”) is repeated to generate the final gridded, merged AOD and AAOD datasets. These final merged AODs and AAODs constructed by assimilating ground-based data from 44 stations for AOD and 34 stations for AAOD are henceforth referred to as MG AOD and MG AAOD, respectively.

4.2 Spatiotemporal characteristics of merged products

Having generated the harmonized gridded datasets of AOD and AAOD, we examined the spatiotemporal features for their fidelity in reproducing the already-reported characteristics over this region from several subregional studies. In Figs. 6 and 7, we present the spatial variation of MG AOD, respectively, for January 2009 (representative of winter, the season with lowest vertical mixing) and May 2009 (representative of the pre-monsoon summer season, when the convective mixing is very strong). The corresponding features for AAOD are shown in Figs. 8 and 9. In all four figures, locations of ground stations, data from which are assimilated, are indicated by circles. In Figs. 6 and 7 (Figs. 8 and 9), the panels from left to right indicate SR AOD (SR AAOD), MG AOD (MG AAOD) and dAOD (dAAOD) which is the difference between merged and satellite-retrieved AOD (SR AAOD).

Mathematically,

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathrm{dAOD}=\mathrm{MG}\mathrm{AOD}-\mathrm{SR}\mathrm{AOD}\text{(16)}& {\displaystyle }\mathrm{dAAOD}=\mathrm{MG}\mathrm{AAOD}-\mathrm{SR}\mathrm{AAOD}.\end{array}$$

Figure 6 Spatial variation of monthly mean SR AOD (a), MG AOD (b) and dAOD (c) for January 2009.

Figure 7 Spatial variation of monthly mean SR AOD (a), MG AOD (b) and dAOD (c) for May 2009.

Figure 8 Spatial variation of monthly mean SR AAOD (a), MG AAOD (b) and dAAOD (c) for January 2009.

Figure 9 Spatial variation of monthly mean SR AAOD (a), MG AAOD (b) and dAAOD (c) for May 2009.

Figures 6 to 9 clearly show that the broad spatial features are consistent between the merged and respective satellite products. For instance, higher AODs and AAODs are exhibited over the IGP than those over central and peninsular India by merged products, which is in line with their respective gridded parents and also with reports from several past studies (Babu et al., 2013).

However, in most of the cases, the merged AODs (Figs. 6b and 7b) show higher values than the respective satellite AODs (Figs. 6a and 7a) as the ground station points are approached. This is in line with the general observation about satellite-retrieved AODs being underestimated over this region (Jethva et al., 2005, 2007; Tripathi et al., 2005). However, sufficiently further away, where no ground-based measurements are available, the merged products tend to be close to the corresponding satellite data. This brings in the need for improving the density of the ground network to further improve the accuracy of regional AOD for providing even better inputs to climate models.

Further, we estimated the variance in merged AOD and AAODs and compared it with that in respective satellite products. As assured by the assimilation methodologies employed, the uncertainties (square root of variance) in merged AOD and AAODs are observed to be substantially lower than those in the corresponding satellite data. For the above-shown representative cases, the uncertainties in merged AODs are observed to be even as small as ≈13 % of those in SR AOD. The uncertainties in merged AAODs are estimated to be as small as ≈82 % and ≈56 % of those in corresponding satellite products during January and May 2009, respectively.

4.3 SSA estimation

The merged gridded datasets of AOD and AAOD over the domain enable estimation of SSA, the critically important aerosol parameter for radiative forcing estimation. The importance of the accurate estimation of SSA for climate impact assessment of aerosols has been underlined by numerous studies in the past (Haywood and Shine, 1995, 1997; Heintzenberg and Helas, 1997; Russell et al., 2002; Takemura et al., 2002; Loeb and Su, 2010; Babu et al., 2016). Takemura et al. (2002) have shown that the small changes in SSA can even alter the sign of aerosol radiative forcing (ARF) at TOA, while Loeb and Su (2010) have demonstrated that uncertainties in SSA could even be the largest contributor to the uncertainties in total direct ARF in clear-sky as well as all-sky conditions.

Even though the OMI SSA (Torres et al., 2007) provides a wide spatial coverage, OMI retrievals suffer from uncertainties emanating largely from subpixel cloud contamination as well as assumptions regarding height of an aerosol layer and surface albedo (Satheesh et al., 2009; Jethva et al., 2014; Torres et al., 2007). On the other hand, the highly accurate SSAs derived from airborne measurements of scattering and absorption coefficients (Babu et al., 2016) are location- and season-specific and thus lack the spatiotemporal coverage necessary for the regional ARF estimation. In this context, the merged and validated gridded AOD and AAOD products, generated above, assume importance.

The gridded data for columnar SSA are derived from the merged AODs and AAODs using Eq. (17) as follows.

$$\begin{array}{}\text{(17)}& \mathrm{SSA}={\displaystyle \frac{\mathrm{MG}\mathrm{AOD}-\mathrm{MG}\mathrm{AAOD}}{\mathrm{MG}\mathrm{AOD}}}\end{array}$$

The spatial variation of the above-estimated SSA is presented for the representative months of January and May 2009 in Fig. 10.

Figure 10 Spatial variation of monthly mean columnar SSA estimated using merged AOD and AAODs for January 2009 (a) and May 2009 (b).

It can be seen from Fig. 10 that SSA over the Indo-Gangetic Plain and northern as well as north-western India is lower than that over southern India during both representative months. This is in line with the regional distribution of SSA reported by Narasimhan and Satheesh (2013) using the gridded SSA retrieved using the joint OMI-MODIS algorithm (Satheesh et al., 2009). The lower SSA values over IGP and north-western India indicate a higher load of absorbing aerosols, mainly BC (over IGP and northern India) and mineral dust (over the north-western Indian region consisting of the Thar Desert). The increased presence of BC over IGP and northern India can be largely attributed to emission from thermal power plants, the increasing number of motorized vehicles as well as biomass burning. Further inspection of Fig. 10 reveals that the eastern coast of India is demonstrating lower SSA with regard to the western coast, especially during the pre-monsoonal month of May 2009 (Fig. 10b). In addition, consistently lower SSA can also be seen over the parts of Myanmar and the surrounding regions during both representative months (Fig. 10).

On the background of sensitivity of aerosol radiative effect to the changes in SSA (Haywood and Shine, 1995, 1997; Heintzenberg and Helas, 1997; Russell et al., 2002; Takemura et al., 2002; Loeb and Su, 2010; Babu et al., 2016), it would be imperative to assess the uncertainty in the above-estimated SSA (Eq. 17). For this purpose, we have perturbed MG AOD and AAOD within their respective uncertainty limits to derive multiple realizations for a given SSA (using Eq. 17) and the standard deviation across these multiple realizations is adopted as an uncertainty in the respective SSA. For the above-shown representative cases (Fig. 10), the root mean square (rms) uncertainty in SSA is 0.03 and 0.02 for January and May 2009, respectively, which is lower than that in OMI SSA (0.05 to 0.1) (Torres et al., 2002) and comparable to that in AERONET SSA (≈0.03 for AOD_440 nm>0.2 and solar zenith angle larger than 50^∘) (Dubovik et al., 2000).

4.3.1 Seasonality in SSA

In view of the known seasonality in aerosol types arising from the seasonal nature of aerosol sources, transport pathways and the mesoscale and synoptic meteorology, it would be important to examine the seasonality of SSA over the study domain in light of already-published data. For this, we have considered four representative and fairly homogeneous subregions of the Indian domain to assess the seasonality, as shown in Table 1 (and depicted in Fig. S1 in the Supplement). In line with the seasonal variation in synoptic meteorology influencing the aerosol field over the Indian region, we have considered three seasons: pre-monsoon, which is comprised of March–April–May and referred to as PrM (characterized by strong heating, deeper planetary boundary layer and prevailing westerlies over the region); winter season, which is comprised of December–January–February (characterized by relatively lesser solar heating, shallower planetary boundary layers and easterly winds); and post-monsoon (referred to as PoM), which is comprised of October–November, which marks the transition from summer monsoon to winter season.

Table 1 Details of subregions considered.

The SSA values derived from the merged datasets following Eq. (17) are averaged over the subregions (Table 1) and seasons are presented in the right-most panels of Fig. 11. The corresponding seasonal mean subregion-averaged values of merged AODs and AAODs appear, respectively, in the left and middle panels of Fig. 11. The error bars represent the standard deviation and hence the spread of the respective quantities in the spatiotemporal domain.

Figure 11 Climatological seasonal cycle of AOD (first column), AAOD (second column) and derived SSA (third column) averaged over IGP (first row), NE (second row), PI (third row) and WAR (fourth row); see Table 1 for acronym definitions.

The figures clearly demonstrate that amongst all the subregions, the highest seasonality in SSA occurs over the IGP (Fig. 11c), while the seasonality is lowest in the NE subregion (Fig. 11f). Over most of the regions, SSA is lowest in the pre-monsoon season and highest in the winter (except for subregion PI; Fig. 11i). The lower SSA over IGP which translates to increased aerosol absorption (Fig. 11c) during the pre-monsoon could be largely because of transport of mineral dust from the Thar Desert and west Asian aid regions to IGP as has been reported earlier (Moorthy et al., 2005; Niranjan et al., 2007; Beegum et al., 2008; Jethva et al., 2005). These mix with the local emissions and get distributed vertically deep in the atmosphere (due to vigorous convective motions in the pre-monsoon season) to the regions above low-level clouds (Satheesh et al., 2008), leading to further enhancements of atmospheric absorption (Chand et al., 2009). It is to be noted that the seasonality of SSA over IGP, as shown in Fig. 11c, is in line with Babu et al. (2016) and Vaishya et al. (2018), who have reported lesser columnar SSA during spring with regard to winter over various locations (Lucknow, Ranchi, Patna and Dehradun) in IGP, based on airborne measurements of scattering and absorption coefficients.

Contrary to all the subregions, the seasonal cycle of SSA over PI shows the maxima (Fig. 11i) during PoM. However, SSA over WAR (Fig. 11l) demonstrates substantial seasonality of a similar kind to that over IGP (Fig. 11c), despite having considerably different seasonal variation in AOD (Fig. 11j) than that over IGP (Fig. 11a).

5 Summary

Gridded datasets of monthly mean aerosol optical depth and absorption aerosol optical depth have been generated for the first time over the Indian region, as a part of the SWAAMI project, by merging long-term measurements from the dense network of ground-based stations with corresponding satellite data. The merging of datasets is performed employing well-established data assimilation methods modified following a weighted interpolation scheme to account for the vertical distribution of aerosols. The gridded data demonstrated improved accuracy and conformity with independent ground-based measurements over different subregions than the corresponding satellite datasets. The merged products also demonstrate substantially less uncertainties than those in respective satellite products, as ensured by the assimilation methodologies employed. These benefits of merged products emphasize their superiority for inputting into regional climate models. The merged AODs and AAODs reproduced the widely reported spatiotemporal features of aerosols over this region despite being significantly different (in terms of AOD and AAOD values) from their gridded parent. The columnar SSA values have been derived from the harmonized products, and their spatiotemporal variation across the domain is examined at regional and subregional scales. The application of these quality-enhanced, merged datasets for regional radiative forcing estimation would be discussed in Part 2 of this two-part paper.

Appendix A: Variance in merged AODs constructed by WIM

As explained in Sects. 1 and 2, satellite-retrieved AOD has higher uncertainties than ground-based AOD measurements, due to several reasons. As the merged AOD product is developed by systematically combining SR and GE AOD by weighted interpolation method, it would be interesting to analyse and compare uncertainty in MG AOD with respect to that in its parent datasets.

For simplicity, we assume a case in which SR AOD at a given grid point (X₁) is being merged with a GE AOD (X₂) from ground-based aerosol observatory lying within the radius of influence from the grid point. So, following the basic equation for WIM, we can write

$$\begin{array}{}\text{(A1)}& \stackrel{\mathrm{̃}}{\mathit{X}}=A{\mathit{X}}_{\mathrm{1}}+B{\mathit{X}}_{\mathrm{2}}.\end{array}$$

Here, $\stackrel{\mathrm{̃}}{\mathit{X}}$ is MG AOD at the given grid point, and A and B are weights (real, positive-valued scalars of size 1×1) for corresponding SR and GE AOD, respectively. As the weighted interpolation method expresses MG AOD as a convex combination of SR and GE AOD, we can write

$$\begin{array}{}\text{(A2)}& A+B=\mathrm{1},\end{array}$$

which means

$$\begin{array}{}\text{(A3)}& A^{\mathrm{2}}+B^{\mathrm{2}}=\mathrm{1}-\mathrm{2}AB.\end{array}$$

Taking variance on both sides of Eq. (A1), we can write

$$\begin{array}{}\text{(A4)}& \mathrm{var}\left(\stackrel{\mathrm{̃}}{\mathit{X}}\right)=A^{\mathrm{2}}\mathrm{var}\left({\mathit{X}}_{\mathrm{1}}\right)+B^{\mathrm{2}}\mathrm{var}\left({\mathit{X}}_{\mathrm{2}}\right).\end{array}$$

Equation (A4) expresses variance in MG AOD as a linear combination of variance in SR and GE AOD with A² and B² being respective weights. As the sum of A² and B² can never be greater than unity, as can be seen from Eq. (A3), the following inferences can be drawn regarding variance in MG AOD.

1. Variance of MG AOD can never be greater than variance of both of its parents.

2. Variance of MG AODs will always be lesser than that of SR AODs if GE AODs are available. Although it is dependent on values of A and B, variance of MG AOD may or may not be lesser than that of GE AOD.

3. If GE AOD is unavailable at a given location (i.e. B=0), then MG AOD and its variance are exactly equal to SR AOD and its variance, respectively.

These observations can be easily verified for the general case in which observations from multiple ground stations are being assimilated with the background data at a given grid point (i.e. X₂ is a vector).

Appendix B: Variance in analysed estimate by 3D-Var

The 3D-Var method constructs an analysis estimate such that the squared, weighted departures in both parents from the analysis estimate are minimized. Here, the weights for departure in each of the parents are expressed as the inverse of the error covariance matrix for the corresponding parent datasets. If both parent datasets are providing unbiased estimates, the analysis estimate constructed by 3D-Var guarantees to have minimum variance. In this section, we prove that this minimum variance for analysis estimate is guaranteed to be smaller than variances in both parent datasets.

Let $\stackrel{\mathrm{̃}}{\mathit{X}}$ be the unknown random variable denoting the analysed (i.e. assimilated) estimate of AAOD with mean M and variance σ². Let X₁ and X₂ be two random variables denoting satellite-retrieved AAOD (referred to as SR AAOD) and ground-based AAOD (referred to as GR AAOD) which are the two available unbiased estimates of AAOD with mean M and variances ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$, respectively. As both SR and GR AAODs are completely independently achieved estimates, we can consider X₁ and X₂ to be uncorrelated.

The goal is to derive the AAOD estimate, $\stackrel{\mathrm{̃}}{\mathit{X}}$, from the linear combination of X₁ and X₂ such that

1. $\stackrel{\mathrm{̃}}{\mathit{X}}$ is a linear, unbiased estimate of AAOD, i.e. $E\left[\stackrel{\mathrm{̃}}{\mathit{X}}\right]=M$ and

2. variance of $\stackrel{\mathrm{̃}}{\mathit{X}}$ is minimum.

Let

$$\begin{array}{}\text{(B1)}& \stackrel{\mathrm{̃}}{\mathit{X}}=a_{\mathrm{1}}{\mathit{X}}_{\mathrm{1}}+a_{\mathrm{2}}{\mathit{X}}_{\mathrm{2}}\end{array}$$

be the assimilated estimate for AAOD where a₁ and a₂ are to be determined such that the above conditions are satisfied. Taking expectations on both sides of Eq. (B1), we get

$$\begin{array}{}\text{(B2)}& {\displaystyle }E\left[\stackrel{\mathrm{̃}}{\mathit{X}}\right]=a_{\mathrm{1}}E\left[{\mathit{X}}_{\mathrm{1}}\right]+a_{\mathrm{2}}E\left[{\mathit{X}}_{\mathrm{2}}\right]\text{(B3)}& {\displaystyle }M=M(a_{\mathrm{1}}+a_{\mathrm{2}});\mathrm{thus},\mathrm{we}\mathrm{get}\text{(B4)}& {\displaystyle }\mathrm{1}=a_{\mathrm{1}}+a_{\mathrm{2}}.\end{array}$$

Taking variance on both sides of Eq. (B1),

$$\begin{array}{}\text{(B5)}& \mathrm{var}\left(\stackrel{\mathrm{̃}}{\mathit{X}}\right)=a_{\mathrm{1}}^{\mathrm{2}}\mathrm{var}\left({\mathit{X}}_{\mathrm{1}}\right)+a_{\mathrm{2}}^{\mathrm{2}}\mathrm{var}\left({\mathit{X}}_{\mathrm{2}}\right).\end{array}$$

We need to find a₁ and a₂ such that $\mathrm{var}\left(\stackrel{\mathrm{̃}}{\mathit{X}}\right)$ is minimized. Therefore, we differentiate Eq. (B5) with respect to a₁ and equate it to zero in order to get expressions for a₁ and a₂ as

$$\begin{array}{}\text{(B6)}& {\displaystyle }{a}_{\mathrm{1}}={\displaystyle \frac{{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}}{{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}+{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}}}\text{(B7)}& {\displaystyle }{a}_{\mathrm{2}}=\mathrm{1}-a_{\mathrm{1}}={\displaystyle \frac{{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}}{{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}+{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}}}.\end{array}$$

Substituting the above expressions for a₁ and a₂ (Eqs. B6 and B7, respectively) in Eq. (B5), we get an expression of variance in $\stackrel{\mathrm{̃}}{\mathit{X}}$ as

$$\begin{array}{}\text{(B8)}& \mathrm{var}\left(\stackrel{\mathrm{̃}}{\mathit{X}}\right)={\left({\displaystyle \frac{{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}}{{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}+{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}}}\right)}^{\mathrm{2}}{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}+{\left({\displaystyle \frac{{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}}{{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}+{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}}}\right)}^{\mathrm{2}}{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}.\end{array}$$

Equation (B8) can be further rearranged as

$$\begin{array}{}\text{(B9)}& \mathrm{var}\left(\stackrel{\mathrm{̃}}{\mathit{X}}\right)=\left({\displaystyle \frac{{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}}{{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}+{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}}}\right).\end{array}$$

It can be verified that

$$\begin{array}{}\text{(B10)}& \left({\displaystyle \frac{{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}}{{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}+{{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}}}}\right)\le min\left({{\mathit{\sigma }}_{\mathrm{1}}}^{\mathrm{2}},{{\mathit{\sigma }}_{\mathrm{2}}}^{\mathrm{2}}\right).\end{array}$$

Therefore, the variance of the linear, unbiased and minimum variance estimator (which is the case for 3D-Var for the present problem) is guaranteed to be smaller that those of parent datasets. This proves that uncertainties (square root of variance) in merged AAODs constructed using 3D-Var are guaranteed to be less than those in SR and GR AAODs.

Appendix C: Validation of MERRA-2 PBLH

We have validated the MERRA-2 PBLHs for the duration of 11 years (2008 to 2018) with those estimated using radiosonde measurements (downloaded from http://weather.uwyo.edu/upperair/sounding.html, last access: 27 June 2019) over the Indian region. However, due to unavailability of continuous radiosonde measurements over many of the locations of ground-based ARFINET and/or AERONET stations, we have considered radiosonde measurements from eight subregional representative locations (Fig. 1), the AOD and BC measurements which are employed for constructing assimilated products. The details regarding lat–long coordinates and broad geographical features for these stations can be found in Tables S1 and S3 in the Supplement, along with other ARFINET and AERONET stations, data from which are used for the assimilation study.

The radiosonde measurements at these stations (Fig. 2) are usually performed twice a day, at 00:00 and 12:00 GMT, and provide vertical distribution temperature, pressure and relative humidity. Further, these fundamental thermodynamic fields are used to derive the vertical profiles for virtual potential temperature (θ_v), which are also provided in the respective data files.

In order to estimate PBLH from the radiosonde data, we have computed the gradient in the virtual potential temperature (Δθ_v) at each given altitude. The height (above surface) at which the Δθ_v exceeds 3 K km⁻¹ is considered as PBLH (Kompalli et al., 2014; Nair et al., 2011) at that location. The planetary boundary layer is likely to be deeper during daytime with regard to nighttime, due to stronger solar heating during the day. Due to this, shallower PBL occurring in the early morning (00:00 GMT) may not be always captured with the provided radiosonde profiles. In view of this, we have employed PBLH estimated using radiosonde measurements during daytime (12:00 GMT) only for the present validation purposes.

The hourly averaged PBLHs (12:00 GMT) given by MERRA-2 for that particular day are bi-linearly interpolated to the locations of stations shown in Fig. 2, in order to get spatiotemporally collocated estimate of MERRA-2 PBLH. The scatter plots between the collocated PBLHs and those estimated from radiosonde measurements for eight locations, during 2008 to 2018, are presented in Fig. 3.

It can be seen from Fig. 3 that PBLHs provided by the MERRA-2 dataset are well correlated with those estimated using radiosonde data, although the correlation coefficient varies from 0.63 to 0.96 with respect to the location. The equations for linear regression between the two PBLH estimates suggest that PBLHs given by MERRA-2 are underestimated over the majority of the stations (Fig. 3a to e), which is in line with the general observation made by the reviewer. Nonetheless, substantially overestimated PBLH values by MERRA-2 are apparent for some of the stations (Fig. 3f to h).