Retrieval of black carbon aerosol surface concentration using integrated MODIS and AERONET data

Black Carbon (BC) is a carbonaceous aerosol that strongly absorbs solar radiation. The high emissions of these highly absorbent particles exacerbate regional air quality and pose significant threats to global climate, both in the short and long term. Therefore, accurately quantifying the spatial distribution of BC is crucial for improving regional air quality and mitigating the climate change impacts driven by human activities. In this study, we developed a novel algorithm for retrieving BC surface concentration jointly using MODIS and AERONET data. Firstly, the algorithm employed the K-means clustering method to determine seasonal background aerosols model based on AERONET V3 daily products. Then, the Maxwell–Garnett effective medium approximation model was utilized to calculate the complex refractive index of the internally mixed aerosols. Subsequently, the lookup tables were established using the 6SV2.1 radiative transfer code to estimate optimal BC fraction and column concentration. Next, the column concentration data were converted to surface concentration using a conversion coefficient derived from MERRA-2. Finally, the retrieved MODIS BC surface concentration was validated with in-situ Aethalometer measurements. The validation showed a correlation coefficient (R) of 0.727, a root mean square error (RMSE) of 0.353, a mean absolute error (MAE) of 0.211, and a linear fit function of $y=\mathrm{0.718}x+\mathrm{0.015}$. These statistical parameters outperform those obtained from MERRA-2 BC data (R = 0.655, RMSE = 0.487, MAE = 0.381, and $y=\mathrm{0.686}x+\mathrm{0.400}$), demonstrating the superior performance of the proposed algorithm in this study area.

1 Introduction

Black carbon aerosol (BC) particles have important impacts on the global climate, air pollution, and human health (Ramanathan and Carmichael, 2008). According to the latest report of the Intergovernmental Panel on Climate Change (IPCC) (Everett et al., 2022), the overall climate effect of aerosols is cooling, but the BC particles therein have a significant warming effect, and play an extremely important role in promoting glacier melting in the Arctic region (Flanner, 2013). In addition, BC emitted from human activities also significantly affects air quality (Cao et al., 2007). It is one of the main components of fine particulate matters (PM_2.5) (Cai et al., 2020), which has a great impact on human health and is considered an important factor leading to obesity (Guo et al., 2022). Therefore, it is of great research significance to obtain accurate and reliable BC spatiotemporal distribution, especially BC surface concentration data that is extremely relevant to human activities.

Currently, the most common method to obtain the chemical composition of BC relying on in-situ measurements, such as Aethalometers (AE31/33) (Pavese et al., 2020), aerosol mass spectrometer (Wang et al., 2018), scanning electron microscopy (Brodowski et al., 2005), etc. These methods are accurate and reliable but can only obtain BC concentrations in specific areas. In addition, using chemical transport models to simulate BC concentrations is also an important method (Xu et al., 2020). Still, the current assimilation products have low spatial resolution and the results are highly uncertain due to the deviation of model input data from actual atmospheric conditions (Sato et al., 2003).

Over the past 20 years, obtaining aerosol physical and chemical properties based on ground-based remote sensing and satellite remote sensing methods has been a hot topic in atmospheric science research (Remer et al., 2024). In ground-based remote sensing, some scholars used the Aerosol Robotic Observation Network (AERONET) to invert the BC-dominated absorption component based on the observation results of the complex refractive index (RI), combined with the three aerosol characteristics of BC, ammonium sulfate, and aerosol water (Sato et al., 2003; Schuster et al., 2005). Later, the single scattering albedo (SSA) was introduced based on the RI, and a five-component inversion model was established (Wang et al., 2013). On this basis, Xie et al. (2017) introduced the log-normal volume size distribution, and then using fine mode volume concentration and coarse mode volume concentration, combined with aerosol optical depth (AOD) correction, successfully separated the proportion of coarse and fine particles, and simulated the concentration and volume distribution of seven different aerosol types. Zhang et al. (2024) simultaneously obtained the long-term series BC column concentration of global AERONET stations based on the Generalized Retrieval of Aerosol and Surface Properties (GRASP)/Component algorithm (Dubovik et al., 2011). However, ground-based remote sensing cannot accurately describe the spatial variation of BC. Therefore, multi-temporal and wide-area observations using satellite remote sensing have the potential to monitor the large-scale spatiotemporal variation of BC. Some studies have preliminarily proposed BC concentration inversion algorithms based on satellite data, such as applying the GRASP algorithm to polarized satellite data such as Polarization and Directionality of the Earth's Reflectances (POLDER) (Bao et al., 2019; Li et al., 2019, 2020) and Directional Polarimetric Camera (DPC) instruments (Li et al., 2022). Based on Moderate Resolution Imaging Spectroradiometer (MODIS) data, the BC column concentration was estimated in China using the lookup table method (LUT) (Bao et al., 2020). Based on the geostationary satellite Himawari-8 data, the hourly BC surface concentration in North China was estimated using the critical reflectance method (Bao et al., 2023). Choi et al. (2024) used the inferred BC volume fraction and particle mass concentration using the Earth Polychromatic Imaging Camera (EPIC) to infer the volume fraction and particle mass concentration of BC based on the spectral absorption provided by the Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm. In this paper, we utilized MODIS data, long-term AERONET aerosol optical property observation data, considering the seasonal differences of background aerosols (BAs) across the study area. The K-means method was applied to categorize the optical properties of BAs for different seasons, and multiple LUTs were established with varying BC fraction. This enables the development of a novel BC surface concentration inversion algorithm tailored to the region. Given that MODIS has provided a substantial volume of long-term data, this new algorithm holds significant potential for investigating long-term spatiotemporal changes in BC concentrations. Section 2 describes data source; Sect. 3 illustrates the methodology, including forward model, calculation strategy for physical properties of internal mixed aerosols, inverse method, and sensitivity studies; Sect. 4 contains the retrieved MODIS BC surface concentration results, validation, and uncertainty analysis. Section 5 represents conclusion.

2 Data

2.1 MODIS data

MODIS has been recording data on the Aqua and Terra satellites launched by NASA and has been providing a large amount of observations since 1999 (Remer et al., 2005). This study used MODIS data from November 2023 to June 2024 in the study area of 5–20° E, 30–50° N, including three types of datasets: MO/YD02 (L1B data), MO/YD03 (Geolocation data), and MO/YD04 (AOD data based on Dark Target algorithm (DT), 0.55 µm). The DT used the linear relationship between the surface reflectance of 0.47, 0.66 and 2.12 µm to retrieve AOD. This product has been widely used in atmospheric remote sensing and climate change research due to its reliable accuracy and long time series. The DT used the linear relationship between the surface reflectance of 0.47, 0.66 and 2.12 µm to retrieve AOD (Levy et al., 2013). This product has been widely used in atmospheric remote sensing and climate change research due to its reliable accuracy and long time series. These datasets can be obtained from this website https://ladsweb.modaps.eosdis.nasa.gov/ (last access: 5 July 2025).

2.2 AERONET data

AERONET is the world's most widely used ground-based aerosol physical characteristics observation network, providing long-term aerosol optical and physical property observation data from thousands of stations for nearly 30 years (Dubovik et al., 2000). This study used the AERONET V3 daily dataset of 32 stations in the study area to obtain BAs characteristic data, which was used as aerosol model input data in the atmospheric radiation transfer model. The locations of these sites are shown as red dots in Fig. 1, and the detailed site information is shown in Table 1. This dataset can be downloaded at https://aeronet.gsfc.nasa.gov/ (last access: 5 July 2025).

Figure 1 Location distribution of AERONET and AE33 ground observation stations. The small red dots, large green dots, and text in the figure represent AERONET sites, AE33 sites, and AE33 site names, respectively.

Table 1 AE33 and AERONET ground observation stations parameters.

2.3 AE33 data

AE33 aethalometer is based on the principle of light absorption and quantifies BC surface concentration by measuring the light absorption characteristics of aerosol samples at multiple wavelengths (Yus-Díez et al., 2021). The instrument typically conducts real-time continuous light absorption measurements at seven wavelengths, ranging from ultraviolet to near-infrared, allowing it to distinguish between different sources of BC and aerosol components, thus improving data accuracy (Rajesh and Ramachandran, 2018). In this study, we used the BC surface concentration data from 6 sites equipped with AE33 and located in the study area, with the measurement wavelength at 637 nm. The locations of the 6 AE33 sites are shown in Fig. 1, and the detailed site information is shown in Table 1. The AE33 BC surface concentration data can be obtained from this website https://ebas.nilu.no/data-access/ (last access: 5 July 2025).

2.4 MERRA-2 data

MERRA-2 is a global atmospheric reanalysis dataset developed by NASA. It is specifically designed to provide high-quality historical datasets for the study of atmospheric and climate processes (Gelaro et al., 2017). In this paper, we used water vapor and ozone data to correct the absorption of MODIS L1B data, and BC column concentration and surface concentration data were used to obtain a priori ratio and for comparison. The MERRA-2 datasets can be downloaded from https://search.earthdata.nasa.gov/ (last access: 5 July 2025).

3 Methodology

3.1 Forward model

In this study, we used the 6SV2.1 model (Vermote et al., 2016), which is widely used in current aerosol remote sensing inversion. When the zenith angle does not exceed 75°, the estimation error of top of atmosphere reflectance (TOA) does not exceed 0.4 % (Jiang et al., 2022). The 6SV2.1 model equation is as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathit{\rho }}_{\mathrm{TOA}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi },\mathit{\tau },f_{\mathrm{BC}},\mathrm{BAs}\right)=\\ & T_{\mathrm{g}}\left[\begin{array}{l}{\mathit{\rho }}_{\mathrm{atm}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi },\mathit{\tau },f_{\mathrm{BC}},\mathrm{BAs}\right)\\ +{\displaystyle \frac{T_{\mathrm{1}}\left({\mathit{\theta }}_{\mathrm{s}},\mathit{\tau },f_{\mathrm{BC}},\mathrm{BAs}\right)T_{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{v}},\mathit{\tau },f_{\mathrm{BC}},\mathrm{BAs}\right){\mathit{\rho }}_{\mathrm{s}}}{\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{s}}S\left(\mathit{\tau },f_{\mathrm{BC}},\mathrm{BAs}\right)\right)}}\end{array}\right]\end{array}\end{array}$$

In the Eq. (1), θ_s, θ_v, φ, τ and f_BC denote solar zenith angle, satellite zenith angle, relative azimuth angle, AOD, and BC volume fraction, respectively. ρ_TOA, ρ_atm, ρ_s, S, T_g, T₁, and T₂ represent TOA, atmospheric path reflectance, surface reflectance, atmospheric spherical albedo, gaseous transmission, downward atmospheric transmission, and upward atmospheric transmission, respectively.

3.2 Estimation of optical properties of mixture aerosol

Maxwell−Garnett effective medium approximation model (MG-MEA) is used to estimate the Maxwell-Garnett dielectric function of aerosol mixtures (Schuster et al., 2005). In this paper, the schematic diagram of the mixture including BC surrounded by BAs is shown in Fig. 2. The MG-EMA equation is as follows:

$$\begin{array}{}\text{(2)}& {\mathit{\epsilon }}_{\mathrm{MG}}=\left[\mathrm{1}+{\displaystyle \frac{\mathrm{3}{f}_{\mathrm{BC}}\left(\frac{\left({\mathit{\epsilon }}_{\mathrm{BC}}-{\mathit{\epsilon }}_{\mathrm{BAs}}\right)}{\left({\mathit{\epsilon }}_{\mathrm{BC}}+\mathrm{2}{\mathit{\epsilon }}_{\mathrm{BAs}}\right)}\right)}{\left(\mathrm{1}-f_{\mathrm{BC}}\frac{\left({\mathit{\epsilon }}_{\mathrm{BC}}-{\mathit{\epsilon }}_{\mathrm{BAs}}\right)}{\left({\mathit{\epsilon }}_{\mathrm{BC}}+\mathrm{2}{\mathit{\epsilon }}_{\mathrm{BAs}}\right)}\right)}}\right]\end{array}$$

where ε_MG, ε_BAs, f_BC and ε_BC indicate mixture aerosol dielectric function, BAs complex dielectric function, BC volume fraction and BC complex dielectric function, respectively. For ε_BAs and ε_BC values can be obtained through RI:

$$\begin{array}{}\text{(3)}& {\mathit{\epsilon }}_j={\mathrm{RI}}_j^{\mathrm{2}}\end{array}$$

where j = BC, BAs, represents different component.

Figure 2 Example diagram of internal mixing of aerosols.

After obtaining ε_MG, the mixture aerosol RI can be calculated by Eqs. (4) and (5):

$$\begin{array}{}\text{(4)}& {\displaystyle }n=\sqrt{{\displaystyle \frac{\left(\sqrt{{\mathit{\epsilon }}_{\mathrm{r}}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{i}}^{\mathrm{2}}}+{\mathit{\epsilon }}_{\mathrm{r}}\right)}{\mathrm{2}}}}\text{(5)}& {\displaystyle }k=\sqrt{{\displaystyle \frac{\left(\sqrt{{\mathit{\epsilon }}_{\mathrm{r}}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{i}}^{\mathrm{2}}}-{\mathit{\epsilon }}_{\mathrm{r}}\right)}{\mathrm{2}}}}\end{array}$$

where ε_r and ε_i are real part and imaginary part of ε_MG, n and k are real part and imaginary part of the mixture aerosol RI.

Figure 3 BAs particle optical properties in different seasons. Panels (a) and (b) denote particle volume size distribution and SSA, respectively. DJF represents winter (December–January–February); MAM represents spring (March–April–May); JJA represents summer (June–July–August); SON represents autumn (September–October–November).

Table 2 RI clustering results of BAs in different seasons.

Table 3 Particle volume size distribution parameters clustering results of BAs in different seasons. Vol-m, VMR-m, and SD-m (m = F, C; F = Fine mode, C = Coarse mode) represent particle volume concentration, volume median radius, and standard deviation, respectively.

We used the ${\mathrm{RI}}_{\mathrm{BC}}=\mathrm{1.95}-\mathrm{0.79}i$ (Bond and Bergstrom, 2006). As regards RI_BAs, we obtained them from AERONET based on K-means cluster method (Russell et al., 2014). In AERONET V3 daily product, we used the following criteria to remove strong absorbing fine aerosol particles data before clustering: (1) at a wavelength of 440 nm, many aerosol particles exhibit strong absorption, making it challenging to separate strongly absorbing BC particles. However, BC still exhibits strong absorption in the 675–1020 nm range, and this characteristic can be leveraged to effectively isolate strongly absorbing fine particles (Bond et al., 2013), so we choose to remove the data with SSA (675–1020 nm) < 0.85 and Fine mode fraction (FMF) > 0.4; (2) in some biomass combustion and industrial cases, SSA values range from 0.85 to 0.95 (Dubovik et al., 2002). In order to reduce the impact of BC aerosols on classification, fine particles whose SSA decreases with wavelength climbing (Ångström Exponent, AE > 1.5) are also removed. Figure 3 shows the particle volume size distribution and SSA of BAs at different times. It can be observed that during spring and summer, the volume concentration of coarse-mode particles is higher in BAs, which is associated with the frequent occurrence of dust aerosols from North Africa during from March to June every year (Meloni et al., 2008). Moreover, the changes in SSA across different seasons are quite pronounced, with the absorption of fine aerosol particles being higher in winter. Using data from all seasons for clustering could introduce significant errors in the estimation of BAs. Therefore, this study clustered the AERONET data by season to obtain accurate seasonal variations in the physical properties of BAs. Tables 2 and 3 show RI clustering results and particle volume size distribution parameters clustering results of BAs in different seasons.

In the 6SV2.1 model, we need to input aerosol mixture RI and particle volume size distribution. The particle volume size distribution equation is as follows (Dubovik and King, 2000):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}\ln r}}=\sum _{i=\mathrm{1}}^n\left({\displaystyle \frac{C_i}{\sqrt{\mathrm{2}\mathit{\pi }}\ln {\mathit{\sigma }}_i}}\right)\exp \left[-\mathrm{0.5}{\left({\displaystyle \frac{\left(\ln r-\ln r_{\mathrm{m},i}\right)}{\ln {\mathit{\sigma }}_i}}\right)}^{\mathrm{2}}\right]\end{array}$$

In Eq. (6), i represents components, including BC, fine BAs, and Coarse BAs; r represents particle radius; C_i, σ_i, and r_m,j represent particles volume concentration, standard deviation, and volume median radius of different components, respectively. Particles volume size distribution parameters for BAs have been shown in Table 3. For BC, $r_{\mathrm{m},\mathrm{BC}}=\mathrm{0.095}$ µm and σ_BC=1.80 µm (Ganguly et al., 2009). Because the C_total is normalization parameter, C_i is equal to the volume fraction of each component.

Figure 4 BC surface concentration retrieval algorithm workflow chart.

3.3 Inverse method

After geometric correction of the MODIS L1B data, the water vapor and ozone data of MERRA-2 were used to correct the gas absorption of the band (Xie et al., 2020). Considering that at 0.47 µm, there is strong absorption by other particles (like BrC et al.) besides BC absorption (Chung et al., 2012), this study used the band of 0.66 and 2.12 µm for inversion. The cloud mask algorithm used multiple wavelengths from visible to near-infrared for cloud identification (Xue et al., 2014). Then we used the 6SV2.1 to build LUT based on mixture aerosol optical properties. Since the fraction of BC in mixture aerosol particles generally does not exceed 6 % (Bao et al., 2020), when generating the LUT, the value of f_BC is from 0 to 0.06, and the step size is 0.01. In the retrieval process, we used DT AOD to input for finding optimal results, the cost function is as follows:

$$\begin{array}{}\text{(7)}& \mathit{\chi }=min\sum _{j=\mathrm{1}}^{\mathrm{2}}{\left({\mathit{\rho }}_{\mathrm{s},j}^{\mathrm{cal}}-{\mathit{\rho }}_{\mathrm{s},j}^{\mathrm{DT}}\right)}^{\mathrm{2}}\end{array}$$

where j represents 2 band (0.66 and 2.12 µm), ${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{cal}}$ and ${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{DT}}$ are calculated surface reflectance and surface reflectance generated by DT algorithm, respectively.

Figure 5 Sensitivity analysis of BC inversion based on 6SV2.1 model. Panels (a)–(d) represent the ρ_s at 0.66 µm variation step sizes, which are 0.02, 0.10, 0.20, and 0.30, respectively. Panel (e) represents SSA and Standard deviation of TOA changes in different f_BC.

Once the f_BC is obtained, the BC column concentration can be calculated according to the following equation:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathrm{BC}}_{\mathrm{column}}=f_{\mathrm{BC}}\times C_{\mathrm{total}}\times {\mathit{\rho }}_{\mathrm{BC}}\text{(9)}& {\displaystyle }{C}_{\mathrm{total}}={\displaystyle \frac{C_{\mathrm{BAs}}\left(\mathit{\tau }\right)}{\left(\mathrm{1}-f_{\mathrm{BC}}\right)}}\end{array}$$

In Eq. (8), ρ_BC=1.0 g cm⁻³ (Ganguly et al., 2009), represents mass density of BC particles. C_BAs is the integral of the volume size distribution of BAs obtained from AERONET clustering.

Since AE33 measures BC surface concentration (BC_surface), the inverted BC_column needs to be converted. Previous studies assumed that BC was continuous uniform distribution below the atmospheric boundary layer, and directly divided BC_column by the boundary layer height to obtain BC_surface (Li et al., 2020; Bao et al., 2019). However, it is well known that the vertical distribution of BC is not uniform (Yuan et al., 2022), so this study used MERRA-2 data to obtain the ratio K at each pixel to improve the accuracy of the conversion. The conversion equation is as follows:

$$\begin{array}{}\text{(10)}& {\mathrm{BC}}_{\mathrm{surface}}=K\times {\mathrm{BC}}_{\mathrm{column}}\end{array}$$

The overall inversion process is shown in Fig. 4.

3.4 Model sensitivity analysis

Based on the aerosol model physical properties obtained above, we took the aerosol model in the DJF period as an example to conduct sensitivity analysis on 6SV2.1. The solar zenith angle, satellite zenith angle and relative azimuth angle are 30, 30 and 12° respectively, and the surface reflectance of 0.66 µm change steps are 0.02, 0.10, 0.20, and 0.30. The relevant results are shown in Fig. 5.

Figure 6 Maps of monthly BC surface concentration distribution. Panels (a)–(h) respectively represent from November 2023 to June 2024.

As shown in Fig. 5a–d, the sensitivity analysis results indicate that as AOD increases, the estimated TOA standard deviation under different surface conditions gradually increases, suggesting that the theoretical inversion accuracy is higher under high aerosol loading conditions. However, when ρ_s=0.10 and the aerosol loading is high (AOD > 1.0), if the BC fraction is high (f_BC≥0.04), the TOA will basically not change with the increase of AOD, which will lead to an unsatisfactory inversion effect under such conditions. In Fig. 5e, as the BC fraction increases, the SSA, which is independent of aerosol load, decreases notably, suggesting that BC content has a substantial impact on the overall aerosol absorption properties. Additionally, under low f_BC conditions, the standard deviation of the estimated TOA for dark surfaces is higher, while under high f_BC conditions, the standard deviation of the estimated TOA for bright surfaces is higher. This indicates that bright surfaces are more sensitive to absorbing aerosols and are more conducive to estimating strongly absorbing BC particles.

4 Results and discussion

4.1 Inverse results

Figure 6 shows the monthly variations in MODIS BC surface concentration from November 2023 to June 2024. It is evident that in the regions surrounding northern Italy, the BC concentration exhibits a pattern of first increasing and then decreasing from November 2023 to March 2024, with emission levels significantly higher than in other areas. This trend is likely related to the region's high population density, developed industry, and low temperatures, which hinder the timely dispersion of emitted BC. Additionally, the presence of the northern and western Alps as well as the southern Apennine Mountains determines weak wind conditions and frequent temperature retrogrades, which hinder atmospheric diffusion and trap pollution on the ground (Renna et al., 2024). From April to June, the overall BC concentration in the study area remains at a relatively low level.

4.2 Validation

To specifically compare the differences between MODIS BC and AE33 BC, this study employed a spatio-temporal matching method. For MODIS, the average value of valid pixel data within a 50 km radius centered on the location of the ground station was used. For AE33, the average value was calculated from the data within 1 h before and after the satellite's transit time (Remer et al., 2005). Figure 7 presents a comparison of the trend changes in the observed valid values of MODIS BC and AE33 BC from November 2023 to June 2024 at three AE33 stations (ISAC Bologna II, Marseille Longchamp, and Milano Pascal) located in high BC emission areas. It is evident that the fluctuation trends of MODIS BC and AE33 BC are generally consistent, although MODIS BC tends to be lower than AE33 BC most of the time. Therefore, the inversion results based on this algorithm can accurately capture the spatiotemporal variations of BC in high-emission areas.

Figure 7 Comparison of consistency changes between MODIS BC and AE33 BC. Panels (a), (b), and (c) represent ISAC Bologna II, Marseille Longchamp, and Milano Pascal, respectively.

Figure 8 Scatter plot of MODIS BC and AE33 versus MERRA-2 BC and AE33. Panels (a) and (b) represent MODIS BC and MERRA-2 BC, respectively. The black dashed line and the green solid line represent the 1 : 1 line and the linear regression function respectively.

Figure 8a presents a scatter plot of all valid values between MODIS BC and the six AE33 BC sites located in the study area from November 2023 to June 2024, used to quantitatively evaluate the accuracy of the algorithm. Additionally, Fig. 8b shows the BC surface concentration verification accuracy of MERRA-2 for comparison. Statistical parameters include the total number of matching points (N), correlation coefficient (R), root mean square error (RMSE), mean absolute error (MAE), and linear fitting function (Jiang et al., 2024). The corresponding accuracy validation results are summarized in Table 4. The correlation coefficient (R) for MODIS BC is 0.727, while for MERRA-2 BC it is 0.655, indicating that our algorithm performs better in terms of relevance. The RMSE for MODIS BC is 0.353, compared to 0.487 for MERRA-2 BC, and the MAE for MODIS BC is 0.211, whereas for MERRA-2 BC it is 0.381. These results suggest that the numerical difference between MODIS BC and AE33 is smaller, indicating better accuracy for MODIS BC. From the perspective of linear regression and scatter point density, MERRA-2 tends to overestimate the surface concentration of BC, while MODIS BC shows the opposite trend. This may be related to the MG-EMA model only considers BC internal mixing state, but there may still be a small amount of fresh and exposed BC externally mixed in the atmosphere (China et al., 2013), which may result in an underestimate of BC.

Table 4 Summary of statistical parameter results from Fig. 8.

4.3 Uncertainty analysis

Figure 9 analyzes the influence of AOD and surface reflectance on the Bias (MODIS–AE33) of the retrieved BC surface concentration. The results reveal that the BC bias exhibits different patterns of change. As shown in Fig. 9a, when AOD is low, the uncertainty of the bias is high, but as AOD increases, the bias decreases. When AOD > 0.75, the overall bias approaches zero, and the uncertainty becomes very low. This trend aligns with the conclusion in Sect. 3.4, where the retrieval accuracy improves under high AOD conditions due to the stronger aerosol signal. In Fig. 9b, the uncertainty of the bias gradually increases as surface reflectance rises. When ρ_s > 0.08, the uncertainty of the bias increases significantly, suggesting that the algorithm's applicability in relatively high surface reflectance areas still needs improvement. However, it is worth noting that due to the lack of data from AE33 stations in high brightness areas (ρ_s≥0.2), the surface reflectance of the AE33 stations in this study is below 0.12, and when ρ_s > 0.1, the data also is insufficient. Therefore, uncertainty analysis confirms that this retrieval algorithm has better performance under high AOD conditions. However, due to the lack of ground-based AE33 observation data in high-brightness surface areas, the accuracy under this surface condition still lacks effective validation.

Figure 9 The box plot of BC surface concentration Bias (MODIS–AE33) independence analysis on the (a) AOD and (b) surface reflectance. The red box represents the interquartile range (IQR, 25th–75th percentiles), and the black whiskers extend to the most extreme data points within 1.5 × IQR from the quartiles.

5 Conclusion

In this paper, we proposed a new algorithm for retrieving BC surface concentrations based on joint MODIS and AERONET data. First, the algorithm derived the optical properties of internally mixed BAs for each season from the AERONET V3 daily product. It used the particle volume size distribution and MG-MEA approximation equation to calculate the RI and volume concentration of the internal mixed aerosol, which were the 6SV2.1 input data of the aerosol model. Then, the sensitivity analysis was conducted for AOD, surface reflectance, and BC fraction. During the inversion process, multiple LUTs were established based on different BC fraction, and DT AOD values were inputted to iteratively find the optimal f_BC. Next, the BC column concentration and surface concentration values generated by MERRA-2 were used to convert the inverted BC column concentration to surface concentration. Finally, the retrieved BC surface concentrations were validated against AE33 observation data. The BC surface concentrations obtained by this algorithm show relatively high reliability and accuracy (R = 0.727, RMSE = 0.353, MAE = 0.211), though there is a slight overall underestimation compared to high-precision ground-based in-situ measurements. This might be due to a small number of BC particles being exposed on the outside of the shell, which led to the failure to estimate the relevant aspect. Additionally, uncertainty analysis of the inversion results indicates that the algorithm is more suitable for high AOD conditions. However, since there is no AE33 site data in the bright surface area, the performance of the inversion results on the bright surface still needs further verification. Therefore, future work will focus on improving the algorithm's performance low aerosol loading conditions and evaluating inversion results accuracy in bright surface.