Constraining the vertical profiles of aerosol component and heating rate over east Asia through assimilating CALIOP observations

The vertical distribution of distinct aerosol components fundamentally governs atmospheric shortwave heating rates and radiative effects. However, traditional data assimilation (DA) methods typically rely on total aerosol optical thickness (AOT) or extinction, which fails to constrain the specific aerosol composition and often leads to the misallocation of aerosol between scattering and absorbing species. To address this limitation, we develop a component-resolved data assimilation system based on the local ensemble transform Kalman filter (LETKF) with spatial and observational constraints within the WRF-Chem model. This novel system assimilates CALIOP-MODIS synergistic retrievals of species-specific extinction profiles (dust, sea salt, black carbon, and water-soluble aerosols) to explicitly optimize the three-dimensional distributions of individual components over East Asia. Results demonstrate that this DA system successfully reconstructs the complex vertical layering of multi-component aerosols. Notably, it effectively corrects severe underestimations of elevated black carbon (BC) plumes, capturing persistent free-tropospheric BC layers over South Asia that traditional models typically miss. Independent validations against ground-based AERONET and AD-Net lidar observations confirm significant improvements not only in total AOT and extinction but also in the single scattering albedo (SSA). By independently adjusting the mass of scattering and absorbing species, the system corrects local biases in aerosol optical properties. Consequently, the optimized component-specific profiles profoundly affect the atmospheric shortwave radiative heating. The elevated BC plumes induce a pronounced mid-tropospheric warming accompanied by a reduction in lower-tropospheric heating due to the attenuation of downward solar radiation. This study highlights the importance of component-specific vertical constraints for accurately assessing aerosol-induced atmospheric heating and its vertical structure.

1 Introduction

Aerosols significantly affect the Earth's energy budget by scattering and absorbing solar shortwave radiation, thereby modulating atmospheric heating rates and radiative effects. Recent research identifies the vertical distribution of aerosols as a dominant driver of uncertainties in these radiative processes (Li et al., 2022). The vertical profile governs the altitude at which solar energy is absorbed and scattered, capable of altering the atmospheric thermal structure (Samset and Myhre, 2011; Samset et al., 2013). Specifically, this profound sensitivity is dictated by the vertically resolved profiles of distinct aerosol components, particularly absorbing species. Over dark surfaces, absorbing aerosols located in the upper atmosphere absorb more incoming solar shortwave radiation than those in the lower atmosphere, leading to stronger warming at higher altitudes. In contrast, over bright surfaces, absorbing aerosols primarily interact with the reflected upward radiation, producing a more pronounced warming effect in the lower atmosphere (Haywood and Ramaswamy, 1998; Mishra et al., 2015). Consequently, the radiative heating induced by black carbon in the stratosphere can exceed its near-surface effect by over an order of magnitude (Samset and Myhre, 2011). Furthermore, atmospheric heating exhibits strong non-linear responses to the relative partitioning of these species; a mere ±3 % variation in the scattering-to-extinction ratio can perturb the radiative effect by 50 % (Loeb and Su, 2010), and multi-component mixing can alter regional radiative effect estimates by up to 100 % (Mishra et al., 2015; Zhang et al., 2022). Therefore, neglecting the distinct vertical mechanisms of individual components introduces substantial uncertainties into atmospheric heating assessments. Despite its critical importance, the severe lack of three-dimensional observational data for aerosol speciation currently constrains our ability to quantitatively assess how component-specific vertical structures affect atmospheric heating.

Current approaches to characterizing aerosol vertical structure primarily rely on active remote sensing and numerical modelling, yet both exhibit significant limitations. Ground-based lidar networks, such as Asian Dust and aerosol lidar observation Network (AD-Net) and the Micro-Pulse Lidar Network (MPLNET), can provide valuable insights into regional vertical distributions (Welton et al., 2002; Shimizu et al., 2004; Wang et al., 2022a). However, the sparse distribution of these stations restricts their horizontal spatial coverage. Satellite-borne lidars, such as Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) and the Cloud-Aerosol Transport System (CATS), face a similar spatial limitation despite providing global high-resolution vertical profiles (Campbell et al., 2010; Proestakis et al., 2019). Due to their nadir-viewing geometry and narrow footprints, these instruments leave unobserved gaps between adjacent ground tracks. Coupled with their temporal discontinuity, these spatial gaps significantly restrict their spatiotemporal representativeness. Although numerical models such as Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) offer continuous spatiotemporal coverage, their simulations still frequently deviate from observations due to uncertainties in emission inventories, coarse grid resolutions, and imperfect sub-grid physical parameterizations (Hu et al., 2010; Koffi et al., 2016; Zeng et al., 2020). Consequently, the lack of continuous and three-dimensional observations, combined with inherent model biases, severely restricts the accurate quantification of aerosol vertical structures and their radiative effects.

Data Assimilation (DA) offers a robust strategy to bridge the gap between observations and models. However, traditional DA studies have predominantly focused on assimilating column-integrated Aerosol Optical Thickness (AOT) (Rubin et al., 2017; Di Tomaso et al., 2017; Yumimoto et al., 2018; Chen et al., 2022; Huang et al., 2023) or surface mass concentrations (Peng et al., 2018; Jin et al., 2022). Since these observations lack vertical information, the assimilation system is forced to distribute increments based on the model's background profile, leaving errors in the vertical structure largely uncorrected (Cheng et al., 2019a; Ye et al., 2021; Gao et al., 2025; Pang et al., 2025). Recent advances in vertical DA using lidar extinction profiles have proven effective in improving three-dimensional aerosol structures (Zhang et al., 2011; Wang et al., 2022a; Escribano et al., 2022) and correcting biases in the boundary layer (El Amraoui et al., 2020). Despite these improvements, a critical gap persists in vertical DA methodologies. Conventional frameworks typically assimilate total aerosol extinction by distributing increments proportionally across different aerosol species, thereby rigidly coupling species abundance to the background model state. This simplification is physically unrealistic, as it neglects the distinct optical properties (e.g., mass extinction efficiencies and single scattering albedos) inherent to different aerosol components. As demonstrated by Tsikerdekis et al. (2021), simply scaling total mass cannot correct structural biases in species-specific optical contributions. Consequently, such DA systems are prone to the misallocation of aerosol mass among different chemical species. For example, the assimilation process may erroneously increase the mass of scattering aerosols to compensate for an underestimation of absorbing species. As highlighted by Yang et al. (2023), relying on a single optical constraint is insufficient to separate the contributions of individual aerosol components. This under-constrained problem ultimately leads to inaccurate vertical profiles of aerosol speciation and introduces subsequent biases in the estimation of aerosol shortwave radiative heating.

Addressing this under-constrained problem is especially important in regions with complex aerosol mixtures, such as East Asia. Characterized by intense anthropogenic emissions and frequent natural dust events, this region exhibits a highly diverse aerosol mixture and complex vertical layering (Li et al., 2018). Because strongly absorbing and scattering aerosols frequently co-exist and overlap vertically, just relying on total extinction constraint here significantly exacerbates the misallocation of aerosol species. Consequently, accurately constraining the component-specific vertical profiles in this region is urgently needed to reduce regional radiative uncertainties and better evaluate atmospheric heating rates.

To bridge this gap, we introduce a data assimilation framework explicitly constrained by vertical profiles of specific aerosol components. While previous studies, such as Tsikerdekis et al. (2021), have successfully performed speciated aerosol assimilation using column-integrated, multi-angle, and polarization measurements to constrain aerosol size and absorption, our approach leverages high-resolution vertical extinction profiles from CALIOP to provide explicit vertical mapping of different aerosol components. By doing so, it overcomes the limitations of traditional bulk-mass approaches and better captures the vertical distribution of individual aerosols, thereby reconstructing the three-dimensional multi-component aerosol structure over East Asia with high accuracy. Beyond structural reconstruction, we quantitatively assess how these optimized vertical profiles govern the radiative response. This study presents a novel pathway for refining component-specific vertical structures via DA, providing critical insights into aerosol-induced atmospheric heating, which is necessary for investigations into regional thermodynamic and dynamic feedbacks. The paper is organized as follows: Sect. 2 details the observational datasets and modelling system, Sect. 3 describes the component-constrained DA methodology, Sect. 4 evaluates the optimized aerosol fields and their radiative implications, and Sect. 5 summarizes the conclusions, discusses the system limitations regarding the background correction and observational uncertainties, and provides future perspectives.

2 Numerical Model and Observational Data

2.1 Forward model

The forward model employed in this study is the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem, version 4.4). The model domain covers East Asia with a 25 km horizontal resolution and utilizes 34 terrain-following sigma levels extending from the surface to 10 hPa. The model physics configuration employs the Thompson microphysics scheme to explicitly account for aerosol-cloud interactions (Thompson et al., 2008). Radiative transfer processes are computed using the Rapid Radiative Transfer Model for GCMs (RRTMG) for both shortwave and longwave bands (Iacono et al., 2008). Boundary layer dynamics are represented by the Yonsei University (YSU) scheme (Hong et al., 2006), coupled with the Unified Noah Land Surface Model (Chen and Dudhia, 2001). Convection is parameterized via the Grell-Freitas ensemble scheme (Grell and Freitas, 2014). The Goddard Chemistry Aerosol Radiation and Transport (GOCART) aerosol module (chem_opt = 300) is used to simulate major tropospheric aerosol species, including sulfate, mineral dust, sea salt, organic carbon (OC), and black carbon (BC) (Chin et al., 2002). Initial and boundary meteorological fields are provided by the ERA5 reanalysis dataset (Hersbach et al., 2020), and the meteorological fields are simulated as a free run without grid nudging. Dust and sea salt emissions are calculated online as a function of surface wind speed. Anthropogenic emissions are mapped from the HTAPv3 mosaic emission inventory for the year 2018 (Crippa et al., 2023), while biomass burning emissions are provided by the Fire INventory from NCAR (FINN) version 2.5 (Wiedinmyer et al., 2023), incorporating a plume-rise parameterization to account for injection heights (Freitas et al., 2007). It should be noted that a portion of the diagnosed model structural deficiencies likely stems from the use of the 2018 HTAPv3 emission inventory for the 2021 study period, which does not account for emission perturbations associated with the COVID-19 pandemic. To capture radiative effects, the model is configured with aerosol-radiation feedback enabled. To prevent chaotic divergence of the meteorological fields over the simulation period while maintaining the aerosol-radiation feedback, a daily re-initialization strategy is employed. Specifically, the meteorological fields are re-initialized every 24 h using ERA5 reanalysis datasets, whereas the aerosol chemical states are continuously cycled and carried forward from the previous cycle. While the aerosol mass simulation follows the original GOCART scheme (Chin et al., 2002), the optical calculations employ the microphysical assumptions inherent to the Community Radiative Transfer Model (CRTM). Aerosols are parameterized as externally mixed spherical particles with lognormal size distributions. For Mie scattering calculations, the effective radius, geometric standard deviation, and complex refractive index of each species are explicitly defined according to the CRTM GOCART-GEOS5 configuration (Lu et al., 2022). To ensure microphysical consistency within the forward operator, the mapping of aerosol mass concentrations from WRF-Chem to CRTM strictly aligns with the GOCART scheme definitions. The effective radii and standard deviations assumed in the CRTM aerosol lookup tables are consistent with those defining the aerosol size bins in the WRF-Chem GOCART module, thereby minimizing internal discrepancies during the mass-to-extinction conversion (Lu et al., 2022; Wang et al., 2022b). These microphysical parameters govern the baseline optical properties and hygroscopic growth, establishing the physical foundation for the observation operator detailed in Sect. 3.1.

2.2 Observation data

2.2.1 CALIOP-MODIS aerosol component retrieval

This study assimilates a global three-dimensional aerosol composition dataset developed by Kudo et al. (2023). This dataset is generated using a synergistic retrieval algorithm that combines attenuated backscatter measurements from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) and passive radiance measurements from the Moderate Resolution Imaging Spectroradiometer (MODIS). It provides vertical profiles of aerosol extinction coefficients, backscatter coefficients, depolarization ratios, and volume concentrations for four specific aerosol components: sea salt (SS), dust (DS), black carbon (BC), and water-soluble (WS).

To generate representative “super-observations” for assimilation, strict quality control and aggregation procedures are applied to the CALIOP-MODIS original orbital retrievals. The data are mapped onto the 25 km WRF horizontal grid by a two-step statistical screening. First, outliers are removed by excluding pixel-level retrievals falling outside one standard deviation of the grid mean. Second, data are retained only for grid cells with sufficient sampling (N>3) and a coefficient of extinction variation below 0.5, thereby filtering out high sub-grid variability (Zhang and Reid, 2006; Schutgens et al., 2017; Cheng et al., 2019b). The screened extinction profiles are then vertically regridded to WRF layers using layer-weighted averaging, requiring a minimum vertical coverage of 20 %. The observational uncertainty ϵ_o is calculated as ${\mathit{ϵ}}_{\mathrm{o}}=\sqrt{{\mathit{ϵ}}_{\mathrm{inst}+\mathrm{retri}}^{\mathrm{2}}+{\mathit{ϵ}}_{\text{rep}}^{\mathrm{2}}}$, where ϵ_inst+retri is the instrumental error and ϵ_rep is the representativeness error. It is important to note that the native pixel-level retrieval uncertainties for speciated aerosol extinction can be substantial, sometimes exceeding 100 % for specific components (Kudo et al., 2023). While spatial aggregation onto the 25 km model grid is performed, these pixel-level errors cannot be assumed to reduce by the inverse square root of the sample size ($\mathrm{1}/\sqrt{N}$) due to the strong spatial correlations inherent in satellite retrievals. Therefore, a baseline of 10 % of the observed aerosol extinction is employed strictly as an empirical tuning parameter for ϵ_inst+retri. In ensemble-based DA frameworks, such empirical tuning is frequently required to yield an optimal Kalman gain and avoid filter divergence. Crucially, the representativeness error (ϵ_rep), derived from the standard deviation of the sub-grid high-resolution retrievals within the 25 km WRF grid cell, is explicitly calculated. In regions with complex aerosol distributions, ϵ_rep typically dominates the error budget (as demonstrated in Fig. S1 in the Supplement), ensuring that the total observation error (ϵ_o) realistically reflects the high uncertainties. It should be noted that ϵ_rep is calculated after the statistical screening steps. This filtering is essential to prevent unphysical retrieval artifacts, including residual cloud contamination and surface reflectance issues, from artificially inflating the observation error. However, it inevitably truncates the natural tails of the sub-grid aerosol distribution. Consequently, the calculated ϵ_rep represents a conservative, lower-bound estimate of the true physical sub-grid variability. The sensitivity of the assimilation performance to this empirical 10 % baseline is further evaluated in Sect. 5.3. It should be noted that this approach likely underestimates the true representativeness uncertainty. As highlighted by Rijsdijk et al. (2025), satellite observations frequently cluster spatially rather than being randomly distributed. Because aerosol fields exhibit strong spatial correlations, a clustered sample captures less sub-grid variability and is less representative of the entire grid cell. While we maintain the present formulation for computational efficiency, accounting for this systematic sampling effect remains an important area for future refinement. Furthermore, the observation error covariance matrix is assumed to be diagonal. This assumption neglects vertical error correlations and cross-covariances between different aerosol species for computational simplicity and to maintain consistency with the variable localization strategy. Finally, hourly aerosol extinction profiles and uncertainties for the four aerosol components are generated at the model resolution for data assimilation. Figure 1 shows the spatial distributions of the gridded AOT integrated from observed extinction coefficients, along with the observation counts for each component used for assimilation during 5–25 March 2021.

Figure 1 Spatial distributions of the gridded Aerosol Optical Thickness (AOT) integrated from observed extinction coefficients, and the corresponding observation counts for each aerosol component used for assimilation (5–25 March 2021).

2.2.2 AERONET

The AErosol RObotic NETwork (AERONET; http://aeronet.gsfc.nasa.gov/, last access: 22 April 2026) provides globally ground-based measurements of aerosol properties (Holben et al., 1998). In this study, AERONET data are used as an independent dataset for the validation of column-integrated aerosol properties. We utilized both AOT and SSA products from the Version 3 (V3) Level 2.0 (quality-assured) dataset (Dubovik and King, 2000; Giles et al., 2019; Sinyuk et al., 2020). Specifically, AOT data are derived from the direct sun algorithm, while SSA retrievals are sourced from the inversion product. AOT at 532 nm is calculated by spectrally interpolating measurements from the 440 and 675 nm channels, while SSA at 532 nm is interpolated using data from the 440 and 675 nm. To match the model temporal resolution, high-frequency observations are averaged within a ±30 min window centered on each hour.

2.2.3 AD-Net

Ground-based lidar measurements from the Asian Dust and Aerosol Lidar Observation Network (AD-Net) (http://www-lidar.nies.go.jp/AD-Net, last access: 22 April 2026; Sugimoto et al., 2003; Shimizu et al., 2004) were used to evaluate the simulated vertical aerosol distributions. We utilized standard AD-Net products, which derive 532 nm aerosol extinction coefficient via the Fernald inversion method (Fernald, 1984). Additionally, the particle depolarization ratio is used to partition total extinction into dust (non-spherical) and spherical components (Shimizu et al., 2004; Nishizawa et al., 2007). To match the model temporal resolution, quality-controlled 15 min observations were averaged hourly (±30 min). For vertical collocation, the observed profiles are linearly interpolated to the model hybrid-sigma levels.

3 Component-resolved LETKF assimilation system

In this study, we developed a component-resolved aerosol data assimilation system with spatial and observational constraints for WRF-Chem (Fig. 2). The system is implemented within a four-dimensional local ensemble transform Kalman filter (4D-LETKF) framework but is configured to perform LETKF analyses in the present experiments. Specifically, the prescribed temporal localization restricts the observations used in each analysis to the current analysis time slot, so the LETKF update effectively corresponds to a three-dimensional LETKF (3D-LETKF), as detailed in Sect. 3.3. While based on a global aerosol assimilation system for the NICAM-SPRINTARS model (Cheng et al., 2019b), this system has been substantially extended to assimilate component-specific vertical extinction profiles rather than total aerosol extinction. To assimilate the multi-species aerosol vertical observations, we develop the new observation operators and implement spatiotemporal and variable localization schemes. Furthermore, severe model underestimations or missed aerosol events often result in near-zero background variance. To address this, we introduce a background correction strategy that maintains sufficient ensemble spread, allowing the system to effectively assimilate the observations and reconstruct unsimulated aerosol events.

Figure 2 Schematic illustrating the component-resolved LETKF system with spatial and observational constraints.

3.1 LETKF system

To generate the background ensemble for the LETKF system, we apply spatiotemporal correlated perturbations to the aerosol emissions rather than perturbing the aerosol concentrations. A total of 20 ensemble members are generated using perturbation factors drawn from lognormal distributions. Based on known emission uncertainties, the means of the perturbation factors are all equal to one, and the standard deviations are set to 350 % for carbonaceous sources (OC and BC) and 130 % for the sulfate precursor SO₂ (Yumimoto and Takemura, 2011). Natural aerosol emissions are also randomly perturbed, assuming uncertainties of 100 % for dust and 500 % for sea salt, consistent with our previous study (Cheng et al., 2019b). Furthermore, since OC, BC, and SO₂ often share common combustion sources, their emission perturbations are assumed to be positively correlated with each other, whereas they remain uncorrelated with the natural aerosol species (Schutgens et al., 2010; Pagowski et al., 2014). Regarding the lateral boundary conditions, the meteorological boundaries are driven by the ERA5 reanalysis. For the chemical boundary conditions, we utilize default climatological profiles rather than time-varying outputs from a global chemical transport model. In the current setup, these chemical boundary conditions are not perturbed. Crucially, for multi-bin species like dust and sea salt in the GOCART scheme, the analysis increments are not simply distributed proportionally to their background mass fractions. Instead, the LETKF independently updates the mass of each size bin based on the flow-dependent background error cross-covariances between the simulated bulk extinction and the individual bin concentrations. This ensemble-driven mechanism allows the assimilation system to dynamically adjust the aerosol size distribution based on the physical uncertainties captured by the ensemble spread.

3.2 Component-specific observation operator

The component-resolved assimilation requires an explicit definition of the aerosol state vector and its mapping to the observed component-specific extinction profiles. For this component-resolved assimilation, the state vector consists of the mixing ratios of 14 GOCART aerosol species: sulfate, hydrophobic and hydrophilic black carbon (BC), hydrophobic and hydrophilic organic carbon (OC), dust in five size bins, and sea salt in four size bins. To link these model states with the component-resolved lidar retrievals, we develop an observation operator that establishes a physical correspondence between the simulated aerosol species and four components of aerosol extinction observations. Unlike traditional operators that aggregate species into bulk fine and coarse modes, our approach explicitly maps individual model variables to their corresponding component observations. Specifically, the simulated mass concentrations of sea salt (M_SS) and dust (M_DU) are mapped directly to their respective extinction coefficients (β_SS and β_DU), and BC mass concentration (M_BC) is aggregated to calculate β_BC. The observed water-soluble extinction (β_WS) is defined as the combined contribution from sulfate (M_SU) and organic carbon (M_OC), formulated as:

$${\mathit{\beta }}_{\text{WS}}(\mathit{\lambda },z)={\mathit{\alpha }}_{\text{OC}}\left(\mathit{\lambda }\right)\cdot M_{\text{OC}}\left(z\right)+{\mathit{\alpha }}_{\text{SU}}\left(\mathit{\lambda }\right)\cdot M_{\text{SU}}\left(z\right),$$

where α represents the mass extinction efficiency (MEE). However, it is important to acknowledge that combining OC and sulfate into a single water-soluble extinction (β_WS) introduces an inherent limitation, as it creates an under-constrained problem at the species level. Since a single optical observation cannot separate two chemically distinct species, the LETKF must partition the analysis increments relying entirely on background ensemble cross-correlations. Because our 20-member ensemble is generated solely by perturbing aerosol emissions under deterministic meteorology, the background cross-correlations between OC and sulfate are primarily driven by the spatial co-location of their emission sources and their subsequent co-transport. While these correlations are physically reasonable for capturing the macro-structure of co-emitted aerosol plumes, we acknowledge that they cannot fully capture the slight differences in aerosol processes between the two species. Consequently, the LETKF modifies the analysis increment based on the background error covariance without a predefined physical preference; rather, the system tends to preferentially adjust the species that exhibits a larger background variance and a higher mass extinction efficiency in the observation operator. To rigorously account for the hygroscopic growth of hydrophilic aerosols, the MEEs are not treated as constants but are dynamically calculated as a function of the model-simulated ambient relative humidity (RH). Specifically, these RH-dependent optical properties are explicitly derived from the Community Radiative Transfer Model (CRTM) GOCART-GEOS5 aerosol lookup tables and spectrally interpolated to the observation wavelength of 532 nm. This dynamic mapping ensures a physically consistent conversion from aerosol mass to optical space under varying meteorological conditions.

However, a fundamental limitation in our observation operator must be acknowledged regarding the microphysical inconsistencies between the forward model and the assimilated retrievals. For the mass-to-extinction conversion in the forward model, the particle sizes and complex refractive indices are strictly defined by the CRTM GOCART-GEOS5 configuration, which treats aerosols as externally mixed pure spherical particles. In contrast, the CALIOP-MODIS retrieval algorithm (Kudo et al., 2023) assumes an internal mixing state for specific components. Most notably, BC is treated using a core-grey shell model, where BC acts as the core internally mixed with water-soluble components as the shell. Consequently, the retrieval employs different assumed size distributions and calculates effective complex refractive indices based on this internal mixing assumption. This fundamental discrepancy in mixing states, sizes, and refractive indices introduces an inherent representativeness error when mapping the simulated externally mixed aerosols to the retrieved internally mixed aerosol extinctions. Nevertheless, it is worth emphasizing that different mixing states primarily affect the aerosol absorption coefficient, whereas their impact on the extinction coefficient is relatively minor. Because our study assimilates the aerosol extinction coefficient rather than the absorption coefficient, the overall impact of this microphysical inconsistency is significantly mitigated. Therefore, despite this theoretical limitation, our assimilation system can still effectively constrain the bulk aerosol mass and vertical structure.

3.3 Localization strategy

We implemented a dual localization framework to mitigate sampling errors inherent in a finite ensemble size and to eliminate spurious correlations. First, spatiotemporal localization restricts the observation influence using a Gaussian function based on the four-dimensional distance:

$$f\left(r\right)=\exp \left(-{\displaystyle \frac{r_{\mathrm{h}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{2}}}}-{\displaystyle \frac{r_{\mathrm{v}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_t^{\mathrm{2}}}}\right),$$

where r_h, r_v, and Δt represent the horizontal, vertical, and temporal distances between the observations and the specific model grid point being analyzed, respectively. The horizontal localization length σ_h, vertical localization length σ_v, and temporal localization scale σ_t are set to 75 km, 1 model layer and 0.01, respectively. The temporal localization scale, σ_t=0.01, is expressed in assimilation time-slot units, with Δt expressed in the same units. This setting acts as a strict temporal filter that only observations assigned to the current analysis time slot contribute to the analysis update, whereas observations in other time slots are effectively excluded. Consequently, although implemented within a 4D-LETKF framework, the LETKF update in the present experiments operates in a 3D-LETKF configuration. The capability to assimilate asynchronous observations is not utilized. Following Miyoshi et al. (2007), we approximate the Gaspari and Cohn (1999) function by truncating Gaussian tails at$r=\mathrm{2}\cdot \sqrt{\mathrm{10}/\mathrm{3}}\cdot \mathit{\sigma }$, discarding observations beyond this cutoff.

Second, we applied component-wise variable localization to prevent unphysical cross-updates between different chemical components. By enforcing a block-diagonal structure on the background error covariance, observations of dust, sea salt, and black carbon are restricted to update only their corresponding mass concentrations. In contrast, water-soluble extinction observations are permitted to simultaneously update both sulfate and OC, consistent with the observation operator. Although physically independent species should theoretically exhibit zero error covariance, a finite ensemble size inevitably introduces sampling errors that lead to spurious cross-correlations. This strict variable localization eliminates these spurious correlations, thereby preserving the physical and chemical consistency of the analysis field.

3.4 Background correction strategy based on spatial and observational constraints

In our system, a 20 member background ensemble is generated by applying spatiotemporal lognormal perturbations to the aerosol emissions (Cheng et al., 2019b). However, severe model underestimations or completely missed aerosol events often result in a near-zero background ensemble mean and spread. It should be noted that these underestimations are not systematic but primarily occur in specific aerosol components (e.g., dust and black carbon), driven by uncertainties in both emissions and meteorological parameterizations such as vertical convection. Because these biases are spatiotemporally heterogeneous, a simple global pre-tuning of emissions is inappropriate. Under such conditions, the standard LETKF would yield a near-zero Kalman gain, effectively ignoring the observations. To mitigate this, we implement a background correction scheme based on spatial and observational constraints prior to the assimilation step. For each analysis grid point, the scheme identifies valid observations within a specified localization radius and applies a distance-dependent localization function, f_i, to weight their contributions. We define a localized weighted bias ratio, R_corr, which quantifies the discrepancy between the observations $y_i^{\mathrm{o}}$ and the background ensemble mean H(x^b)_i within the local domain:

$$R_{\text{corr}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N_{\text{obs}}}\sqrt{f_i}\cdot \frac{H(x^{\mathrm{b}}{)}_i}{y_i^{\mathrm{o}}}}{{\sum }_{i=\mathrm{1}}^{N_{\text{obs}}}\sqrt{f_i}}},$$

where N_obs represents the number of valid observations within the localization radius. To ensure statistical robustness and filter out noise, the correction is triggered only when strict criteria are met: the number of local observations must exceed a minimum threshold (N_obs>3), and the calculated ratio must fall within three standard deviations of the local statistics (3-sigma rule) to exclude outliers. These identified outliers are considered invalid and are completely excluded from both the background correction and the subsequent LETKF assimilation steps. To maintain numerical stability when the background ensemble mean is close to zero (e.g., in the case of completely missed aerosol plumes), two thresholds are applied to the scaling factor R_corr before the correction step: (1) a minimum threshold of $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{15}}$ is set to prevent division by zero; (2) during the 3-sigma outlier check, the lower bound for valid R_corr values is constrained to 0.001. If the calculated R_corr is less than 0.001, the background correction is not applied at that specific grid point. This constraint effectively prevents spurious inflation of the aerosol state when the background concentrations are negligible. In regions where the weighted ratio R_corr falls below an empirical threshold of 0.2, thereby representing severe underestimation, the background ensemble members are rescaled by dividing by this factor. Crucially, when this correction is triggered, the rescaled ensemble directly serves as the analysis field, bypassing the LETKF assimilation step for that specific location. Conversely, in regions where this correction is not triggered, the system proceeds with the standard LETKF assimilation to optimally combine the background and observations. This conditional strategy effectively corrects severe negative biases, where the standard filter would otherwise fail due to near-zero background spread, while preserving the spatial and multivariate consistency of the ensemble.

3.5 Experimental design

To evaluate the impact of data assimilation, two experiments are conducted: a free run without assimilation (FR) and an assimilation run (DA) assimilating the specific aerosol component extinction coefficients. Both experiments spanned from 1 March to 25 March 2021, with the first five days serving as a spin-up period. Furthermore, to mitigate filter divergence and maintain adequate ensemble spread, a fixed multiplicative inflation factor of 1.1 is applied to the background ensemble perturbations, following Sekiyama et al. (2010).

4 Results

This section presents a comprehensive evaluation of the LETKF aerosol component assimilation system, systematically assessing its performance from the reconstruction of aerosol three-dimensional structure to its ultimate impact on radiative properties. Section 4.1 presents a self-validation using the assimilated CALIOP observations to verify the system's ability in reconstructing the three-dimensional distributions of individual aerosol components. To ensure the robustness of this assimilation system, Sect. 4.2 conducts an independent validation using ground-based AERONET and NIES Lidar measurements. This independent evaluation assesses the ability of the assimilation system to reproduce aerosol spatial distributions and temporal evolution, including at locations and times not directly sampled by the assimilated CALIOP observations. Beyond bulk optical depth, Sect. 4.3 further compares the simulated single scattering albedo (SSA) against independent AERONET retrievals, demonstrating the effectiveness of component-specific assimilation in adjusting aerosol scattering and absorption characteristics. Finally, Sect. 4.4 quantifies the impact of optimized aerosol analysis fields on the vertical profiles of the aerosol radiative heating rates. To quantify the model performances, statistical metrics including the bias (BIAS), root mean square error (RMSE), and correlation coefficient (CORR) are calculated between the simulations and observations (Boylan and Russell, 2006; Willmott et al., 2012; Yumimoto et al., 2017).

Figure 3 Density scatter plots comparing the modeled extinction coefficients (km⁻¹) against observations for (a, e) sea salt, (b, f) dust, (c, g) black carbon, and (d, h) water-soluble aerosols. The top row (a–d) represents the FR experiment, while the bottom row (e–h) shows the results from the DA experiment. The color scale indicates the data density. The solid black line represents the 1:1 line while the dashed lines represent the 2:1 and 1:2 lines. Statistical metrics (N: sample size, BIAS: bias, RMSE: root mean square error, CORR: correlation coefficient) are provided in each panel. The statistics are calculated using all valid data points within East Asia (shown in Fig. 1) from the surface up to 12 km altitude.

4.1 Reconstruction of three-dimensional aerosol composition

In this section, we first evaluate the model performance against the assimilated aerosol component extinction coefficients at 532 nm as self-validation (Benedetti et al., 2009). Figure 3 presents scatter density plots comparing the simulated aerosol extinction coefficients from the FR and DA experiments against CALIOP observations for four aerosol components: sea salt, dust, black carbon (BC), and water-soluble (WS) aerosols. As shown in Fig. 3a–d, the FR experiment exhibits a systematic underestimation across all aerosol components. A prominent feature of this negative bias is a high-density region where simulated extinction coefficients approach zero, indicating that the model completely misses aerosol events clearly captured by CALIOP observations. This deficiency is most pronounced for BC (Fig. 3c), where the simulation entirely failed to capture the observed spatio-temporal variations, yielding a notably low CORR of 0.10. Compared with the spatial distributions of the CALIOP observations, it is evident that these elevated BC extinction values are predominantly concentrated over South Asia (Fig. 1). The severe omission of these aerosol plumes in FR experiment demonstrates the background model's inability to represent intense BC emissions and accumulation processes, likely due to significantly underestimated sources in the emission inventories or deficiencies in vertical transport. The component-resolved assimilation system (Fig. 3e–h) significantly enhances model performance, effectively mitigating the systematic underestimation observed in FR experiment across all components. For dust, the DA experiment yields a substantial correction (with the CORR increasing from 0.33 to 0.79 and the RMSE decreasing from 0.090 to 0.054), optimizing the three-dimensional structure of dust aerosols especially in the Taklamakan Desert and Thar Desert (Fig. 1). Similarly, the scatter density distribution of WS aerosols becomes markedly more aligned with the 1:1 line (CORR improving from 0.41 to 0.79 and RMSE decreasing from 0.161 to 0.094), indicating an effective correction of secondary aerosol mass concentration. For sea salt, the assimilation effectively reduces the statistical spread (CORR increasing from 0.33 to 0.69 and RMSE decreasing from 0.033 to 0.022), mitigating random errors driven by uncertainties in sea salt emission parameterizations and marine surface wind fields. For BC, the DA experiment increases the CORR from 0.10 to 0.59 and partially captures the high concentration pollution plumes over South Asia (Fig. 1). However, this improvement remains limited compared to other aerosol species, particularly at high concentrations. This limitation is directly due to the severe negative bias in the background field, exposing a fundamental constraint of ensemble-based assimilation: when the forward model fails to simulate a process due to severe structural errors (e.g., missing emissions or failed vertical lofting), the background ensemble spread approaches zero. Consequently, the assimilation system cannot generate sufficient background error covariance to effectively draw the analysis toward the observations. Although our system includes a correction procedure based on spatial and observational constraints to mitigate severe underestimation, it cannot fully compensate for these fundamental model deficiencies. Therefore, while data assimilation effectively corrects model states given a reasonable background, accurate emission inventories and physical parameterizations remain indispensable for reliably simulating aerosol variations. Furthermore, it is worth noting that the relatively larger BIAS and lower CORR for BC compared to other species (Fig. 3g) can be partially attributed to the inherent discrepancies in microphysical assumptions between the model and the retrieval algorithm. The GOCART model assumes externally mixed pure BC, whereas the CALIOP-MODIS retrieval treats BC as internally mixed with other substances in a core-shell structure. This difference in mixing state significantly alters the assumed optical properties, introducing an unavoidable representativeness error when comparing the simulated and retrieved BC extinctions.

Figure 4 Vertical distribution of extinction coefficients (km⁻¹) for (a) sea salt, (b) dust, (c) black carbon, and (d) water-soluble aerosols. The black solid lines represent the assimilated CALIOP observations. The blue dashed lines and red solid lines denote the results from the FR and DA experiments, respectively. Statistical metrics (BIAS, RMSE, and CORR) calculated for the vertical profiles are listed in each panel. The profiles represent the domain-averaged values over East Asia (shown in Fig. 1).

Figure 5 Comparison of aerosol extinction coefficients (km⁻¹) observed by CALIOP observations, simulated by FR, and DA experiments in four vertical layers (0–1 km, 1–3 km, 3–6 km, and 6–12 km). (a–c) Mean aerosol extinction coefficients for each layer. (d–f) Fractional contribution (%) of each aerosol species to the total extinction. The colors represent different aerosol components: sea salt (blue), dust (brown), black carbon (black), and water-soluble (green).

Figure 4 compares the domain-averaged vertical profiles of aerosol extinction coefficients. The FR experiment systematically underestimates the extinction of all four species across all altitudes. Furthermore, while the FR experiment captures the general vertical distribution pattern of dust, it fails to reproduce the observed vertical structures of the other three components. The DA experiment, however, effectively reconstructs these vertical profiles and substantially reduces the negative biases. Specifically, for dust and water-soluble aerosols, the profiles in DA experiment closely match the observations, yielding CORR of 1.00 and substantially reduced RMSEs. For sea salt, the FR simulation produces a false minimum at approximately 0.5 km instead of the observed near-surface maximum. The DA experiment corrects this profile, properly capturing the low-level accumulation. Regarding BC, the FR profile decreases monotonically with height, completely missing an observed elevated layer peaking at around 6 km. As indicated by the spatial distribution of BC AOT (Fig. 1), this mid-tropospheric peak is primarily driven by intense pollution plumes originating from the Indian subcontinent, which are subsequently lofted into the free troposphere. Chemical transport models frequently fail to reproduce such elevated BC layers due to inherent uncertainties in bottom-up emission inventories and imperfect sub-grid convective parameterizations. Consequently, the presence of these persistent free-tropospheric BC layers over the Indian subcontinent has historically been documented almost exclusively through sparse in-situ profiling using high-altitude balloons and aircraft (Babu et al., 2011; Sunilkumar et al., 2025). The successful reconstruction of this lofted BC plume by the DA experiment (improving the CORR from 0.26 to 0.95) represents a substantial advancement. It is important to distinguish the respective roles of the assimilation components in this reconstruction. As further detailed in Sect. 5.2, the standard LETKF update substantially recovers the vertical shape of the BC profile, whereas the empirical background correction provides the critical magnitude adjustment for the elevated layer where the background ensemble spread approaches zero. However, it should be noted that this represents a self-validation, as the evaluation is performed against the same observations that were assimilated. It underscores the capability of the component-resolved assimilation system to mitigate systematic biases originating from emission deficiencies and to accurately capture complex vertical transport and accumulation processes in the free troposphere. As illustrated in the Fig. S2 in the Supplement, the observations exhibit a distinct and strong lofted BC plume between 15 and 30° N, extending vertically up to approximately 4 km. The FR experiment significantly underestimates the magnitude of this lofted plume, with most of the higher BC extinction values confined near the surface. In contrast, the DA experiment successfully reconstructs the observed lofted plume structure. It captures both the enhanced magnitude and the vertical extent of the BC plume much better than the FR. This cross-sectional analysis clearly demonstrates that the DA not only improves average statistics but also effectively corrects the vertical distribution of aerosols, successfully reproducing the lofted plumes observed in the real atmosphere. Despite the structural improvement, DA experiment still underestimates the peak magnitude of the BC extinction. Interestingly, the FR simulation produces higher BC extinction near the surface compared to the DA experiment (Fig. 4c). This suggests that the failure to capture the 6 km BC peak is not primarily due to underestimated total emissions, but rather attributable to deficiencies in the model's vertical transport mechanisms or an underrepresentation of the biomass burning plume injection heights. Because our assimilation system generates ensemble spread solely by perturbing aerosol emissions while running the meteorology deterministically, the background error covariance lacks flow-dependent information regarding vertical transport uncertainties. Consequently, when the model fails to loft BC to the 6 km altitude, the ensemble spread at this level approaches zero. The DA system can only redistribute mass based on the prescribed emission-driven covariance, limiting its ability to fully correct the altitude of lofted plumes. This limitation explains why the DA experiment captures the shape of the 6 km BC peak but still underestimates its absolute magnitude (Fig. 4c).

Figure 5 further quantifies these vertical discrepancies by integrating the absolute aerosol extinction and fractional contributions across four altitude layers. Consistent with Fig. 4, the FR experiment (Fig. 5b) systematically underestimates the total aerosol extinction at all altitudes compared to CALIOP observations (Fig. 5a). The DA experiment (Fig. 5c) effectively mitigates this negative bias, particularly in the lower-to-mid troposphere (0–6 km), although a slight underestimation persists in the mid-troposphere. Beyond absolute extinction, the fractional analysis (Fig. 5d–f) highlights the capability of data assimilation to correct the vertical aerosol composition. A major deficiency in the FR simulation (Fig. 5e) is the spurious dominance of water-soluble aerosols across all altitude layers, which artificially suppresses the BC and dust fractions. The DA experiment successfully corrects this vertical partitioning. Specifically, in the lower free troposphere (1–3 km), the DA experiment reproduces the observed balanced contribution among dust, BC, and water-soluble aerosols, eliminating the water-soluble overestimation in the FR simulation. In the mid-troposphere (3–6 km), observations indicate that BC becomes the dominant aerosol type (approximately 50 %). The FR simulation misses this feature, maintaining an artificially high water-soluble fraction (up to 78 %). By assimilating the observations, the DA experiment successfully captures the elevated BC fraction and reduces the excessive water-soluble contribution, yielding a composition profile consistent with CALIOP observations. Furthermore, while the FR experiment severely underestimates the dust fraction (remaining below 15 % at all altitudes), the DA experiment effectively reproduces the observed substantial dust contributions (20 %–35 %) throughout the vertical column. Overall, these results demonstrate that the DA method can effectively constrain both the total aerosol loading and the vertical partitioning of different aerosol species.

4.2 Spatiotemporal assessment using independent AERONET and Lidar observations

To evaluate the assimilation performance against unassimilated datasets, we conducted an independent validation using ground-based AERONET observations. Figure 6 presents the density scatter plots and the frequency distribution of AOT deviations for all available AERONET stations during the study period. The FR experiment exhibits a systematic underestimation of AOT, characterized by a negative BIAS of −0.278, an RMSE of 0.521, and a CORR of 0.573. The DA experiment effectively mitigates this negative bias, reducing the BIAS to −0.208 and the RMSE to 0.488, while improving the CORR to 0.626. The frequency distribution of AOT deviations further confirms this improvement. The DA experiment narrows the error distribution, yielding a higher frequency of near-zero deviations (e.g., the fraction within the $-\mathrm{0.10}<\text{deviation}<\mathrm{0.10}$ range increases from 23.36 % to 28.10 %) and notably reducing the frequency of large positive deviations (>0.5) from 21.58 % to 16.70 %.

Figure 6 Density scatter plots comparing the modeled AOT against AERONET observations for (a) the FR and (b) the DA experiments. The color scale indicates the data density. The solid black line represents the 1:1 line while the dashed lines represent the 2:1 and 1:2 lines. Statistical metrics (N: sample size, BIAS: bias, RMSE: root mean square error, CORR: correlation coefficient) are provided. (c) Frequency distribution of AOT deviations (Modeled minus Observed). The blue and red lines denote the results from the FR and DA experiments, respectively. The percentages indicate the proportion of samples falling within specific deviation ranges.

The spatial distribution of the assimilation impact is illustrated in Fig. 7. While observations show high aerosol loadings over South Asia, the FR experiment displays widespread negative biases and high RMSE values across most sites. Figure 7d and e highlight the percentage improvements of the DA experiment relative to the FR simulation in terms of absolute BIAS and RMSE, respectively. Statistically, the DA experiment yields lower RMSE values at 42 out of 59 sites (71.19 %) and reduces the absolute BIAS at 54 sites (91.53 %).

Figure 7 Spatial distributions of observed AOT and model evaluation metrics against AERONET observations. (a) The observed mean AOT averaged over the study period. (b) Mean Bias (BIAS) and (c) Root Mean Square Error (RMSE) for FR experiment. The colored circles represent the locations of observation sites, with the color scale indicating the magnitude of the variable. (d) The percentage improvement of the DA experiment relative to the FR simulation in terms of absolute BIAS (defined as $\frac{\left(\right|{\text{BIAS}}_{\text{FR}}|-|{\text{BIAS}}_{\text{DA}}\left|\right)}{\left|{\text{BIAS}}_{\text{FR}}\right|}\times \mathrm{100}\mathit{\%}$), and (e) the percentage improvement in RMSE (defined as $\frac{({\text{RMSE}}_{\text{FR}}-{\text{RMSE}}_{\text{DA}})}{{\text{RMSE}}_{\text{FR}}}\times \mathrm{100}\mathit{\%}$).

Time series analysis at four representative AERONET stations (Fig. 8) further demonstrates the temporal consistency of the DA improvements. Regionally, DA successfully corrects the FR's underestimation of high-AOT events in source regions (e.g., Beijing-CAMS and Dhaka_University) while maintaining strong agreement with ground observations in downwind areas (e.g., Fukuoka and Gangneung_WNU). This confirms that assimilation effectively constrains both local aerosol concentrations and long-range transport. Crucially, these improvements highlight the capability of the assimilation system to overcome the limitations of the satellite's narrow swath and sparse sampling. Through horizontal localization, the observational information is effectively spread spatially to update adjacent grid cells. Furthermore, the optimized initial conditions are continuously propagated spatiotemporally, driven by atmospheric transport and model dynamics. Consequently, the assimilation benefits extend well beyond specific CALIOP footprints and overpass times, yielding a domain-wide enhancement in the spatiotemporal representation of aerosols.

Figure 8 Time series comparison of AOT between model simulations and AERONET observations at four selected sites during 5–25 March 2021. The panels show results for (a) Beijing-CAMS, (b) Dhaka University, (c) Fukuoka, and (d) Gangneung-WNU. The blue and red solid lines represent the time evolution of AOT from the FR and DA experiments, respectively, while the black open circles denote the corresponding AERONET observations. Statistical metrics, including the Root Mean Square Error (RMSE) and Correlation Coefficient (CORR), are inset within each panel for both experiments.

Furthermore, we performed another independent validation of the vertical aerosol structures using the ground-based AD-Net lidar network. Figure 9 compares the simulated time-height cross-sections of aerosol extinction coefficients with these ground-based lidar measurements at the AD-Net Fukue site (32.75° N, 128.68° E) from 5 to 25 March 2021. Observations reveal a complex vertical structure featuring intense aerosol layers predominantly confined below 1 km, accompanied by episodic elevated transport plumes extending up to 2–4 km (e.g., during 9–11 and 13–16 March). The FR experiment severely underestimates the extinction magnitude across all layers, missing both the intense near-surface aerosol loadings and the elevated plumes. In contrast, the DA experiment significantly reduces these negative biases, successfully reproducing the high extinction in the lower troposphere and the vertical evolution of the transport episodes. Notably, as CALIOP tracks are sparse over Japan (Fig. 1), this improvement at the downwind Fukue site is not driven by local data assimilation. Instead, it results from the optimized three-dimensional aerosol initial fields over the upstream Asian continent. Unlike column-integrated AOT assimilation, assimilating speciated CALIOP profiles provides essential vertical and compositional constraints. The accurate reproduction of the vertical structure at Fukue confirms that optimizing upstream initial conditions effectively improves the simulation of subsequent long-range transport and downwind aerosol distributions.

Figure 9 Comparisons of total observed aerosol extinction coefficients from AD-Net (a) and simulated coefficients for FR (b) and DA (c) experiments from 5 to 25 March 2021 at Fukue site (32.75° N, 128.68° E). (d) The time-averaged vertical profiles of aerosol extinction coefficients for the period.

4.3 Improvement in single scattering albedo and absorption

To further evaluate the assimilation performance on aerosol optical properties, we conducted an independent validation of the simulated single scattering albedo (SSA). Figure 10 compares the simulated SSA with independent AERONET observations. Interestingly, while the model background generally underestimates the black carbon fraction over the broader domain (as indicated by the CALIOP comparisons), the FR experiment exhibits a negative SSA bias at the AERONET sites. This suggests the local overestimation of the absorbing aerosol fraction in the column-integrated background at these specific locations. This contrasting bias between the local overestimation at AERONET sites and the regional underestimation compared to CALIOP can likely be attributed to the imbalance in the modeled aerosol lifecycle. Since AERONET sites are generally located closer to anthropogenic emission sources than the broader spatial coverage of CALIOP, this suggests that the emissions of absorbing aerosols might be overestimated, while their atmospheric lifetimes are too short. Consequently, these aerosols are heavily concentrated near the sources but fail to be adequately transported to broader regions. The DA experiment effectively corrects this negative bias, reducing the RMSE from 0.091 to 0.059. Furthermore, DA experiment significantly increases the fraction of data points with absolute errors below 0.05 (from 30.21 % to 60.28 %), demonstrating improved simulations of aerosol absorption properties. This improvement highlights the advantage of assimilating speciated aerosol observations. Traditional assimilation of total aerosol extinction primarily scales the total aerosol mass, which largely preserves the background aerosol composition ratio and thus struggles to correct such localized SSA biases. In contrast, assimilating CALIOP speciated extinctions adjusts the 3D distribution of individual components. Through this spatially varying correction, the DA experiment successfully optimizes the relative proportions of scattering and absorbing aerosols at the AERONET sites, thereby improving the simulated SSA. It should be noted that the modeled SSA in the DA experiment is slightly underestimated compared to AERONET. This is an expected behavior inherited from the assimilated observations, as the synergistic CALIOP-MODIS retrieval algorithm tends to slightly overestimate the BC fraction, leading to a lower SSA relative to AERONET (Kudo et al., 2023). As shown in Fig. 11, by independently adjusting individual aerosol components based on CALIOP speciated observations, the DA experiment successfully reduces both SSA BIAS and RMSE at 33 out of 36 AERONET sites (91.67 %), yielding mean improvement rates of 33.56 % and 29.82 %, respectively.

Figure 10 Same as Fig. 6 but for single scattering albedo (SSA).

Figure 11 Same as Fig. 7 but for SSA.

Figure 12 presents the time series of simulated and observed SSA at four representative sites. At South Asian (Gandhi_College and Dhaka_University) and Southeast Asian (Mandalay_MTU) sites, the FR experiment consistently underestimates SSA. As shown in Fig. 1, these regions are characterized by high loadings of BC and water-soluble aerosols, primarily driven by intense anthropogenic emissions and spring biomass burning. Although the domain-averaged background generally underestimates the BC fraction (as shown in Fig. 4), the persistent negative SSA bias at these specific locations indicates a local underestimation of the scattering-to-absorbing aerosol ratio. By assimilating CALIOP speciated observations, the DA experiment effectively constrains the concentrations of individual aerosol species, correcting this local imbalance. Consequently, DA significantly increases the simulated SSA, reducing the RMSE from 0.071 to 0.034 at Gandhi_College and from 0.060 to 0.037 at Mandalay_MTU. Similarly, at the downwind Fukuoka site in East Asia, aerosols are a complex mixture influenced by the long-range transport of dust and anthropogenic pollutants. The DA experiment also successfully corrects the underestimated SSA, reducing the RMSE from 0.100 to 0.052. The occasional sharp drops at Gandhi_College and Dhaka_University are primarily driven by the data assimilation process during specific satellite overpasses. When the satellite captures a dense plume of highly absorbing aerosols (like BC or dust) and the assimilation using this data, the system instantaneously adjusts the scattering-to-absorbing aerosol ratio, leading to a sudden drop in the simulated SSA at that specific time step. These results demonstrate that the speciated assimilation dynamically optimizes the aerosol composition, leading to a more accurate representation of aerosol optical properties across both source and downwind regions. It is worth noting that traditional assimilation of total AOT lacks direct observational constraints on vertical speciation, relying heavily on the model's background covariance to distribute increments, which can sometimes struggle to correct the aerosol composition ratios. Consequently, pure AOT assimilation struggles to correct structural biases in aerosol optical properties like SSA. In contrast, by explicitly assimilating component-specific extinction profiles, our DA system independently adjusts the mass of scattering and absorbing species, thereby successfully correcting the local scattering-to-absorbing ratio and significantly improving the simulated SSA.

Figure 12 Same as Fig. 8 but for SSA.

4.4 Radiative heating response to optimized aerosol vertical profiles

Figure 13 illustrates the mean vertical profiles of the shortwave (SW) heating rate along the CALIOP tracks. To isolate the instantaneous radiative impact of the assimilated aerosols, the DA heating rates are diagnosed by updating the model restart files with the analysis aerosol fields and integrating the radiation scheme for a single time step. Compared to the FR experiment, the DA experiment produces a distinct warming effect in the lower to middle troposphere (2–6 km). The track-averaged DA-minus-FR heating-rate difference reaches a maximum of approximately 0.3 K d⁻¹ at 4–5 km (Fig. 13b). As indicated by the vertical aerosol distributions in Fig. 4, this mid-level heating peak can be physically interpreted as the superposition of the elevated BC layer (peaking around 6 km) and the lower-level dust layer (peaking near 1 km). Since the shortwave heating rate is determined by both the aerosol absorption capacity and the available downward solar radiation flux, the combination of these absorbing aerosols culminates in the strongest heating response at 4–5 km. It is important to acknowledge how the discrepancies in aerosol mixing state assumptions between the retrieval algorithm and the forward model affect the diagnosed heating rates. The CALIOP-MODIS retrieval assumes an internal core-grey shell mixing state for black carbon (BC) (Kudo et al., 2023), which yields the SSA at 532 nm ranging from 0.44 (dry conditions) to 0.64 (wet conditions at RH=80 %). In contrast, the GOCART optical scheme in CRTM assumes externally mixed, pure BC particles (Lu et al., 2022), which possess a much lower SSA of approximately 0.20 due to the lack of a scattering shell. When the data assimilation system constrains the total extinction to match the observations, the absorption coefficient is determined by the fraction (1−SSA). The ratio of the absorption fraction between the internally mixed retrieval and the externally mixed model, defined as $R_{\text{abs}}=\frac{(\mathrm{1}-\text{SSA}{)}_{\text{int}}}{(\mathrm{1}-\text{SSA}{)}_{\text{ext}}}$, ranges from 0.45 to 0.70. Because R_abs<1, applying the externally mixed optical scheme in the model probably overestimates the absorption fraction per unit extinction. Figure 14 presents the relationship between the assimilation-induced increments in extinction coefficient and SW heating rate. The samples are categorized into three regimes based on their Single Scattering Albedo (SSA). The scatter plot and corresponding linear regressions reveal a clear, quantitative dependence of heating efficiency on aerosol absorption properties. For highly absorbing aerosols (SSA<0.85), a small increment in extinction leads to a sharp increase in the SW heating rate, as indicated by the steepest slope. As the SSA increases to intermediate values ($\mathrm{0.85}<\text{SSA}<\mathrm{0.92}$), the slope decreases, reflecting a moderate heating response. For predominantly scattering aerosols (SSA>0.92), the heating efficiency per unit extinction is the lowest among the three categories. This distinct stratification confirms that atmospheric SW heating is not only determined by total aerosol loading, but also fundamentally governed by the proportion of absorbing components (e.g., black carbon and dust) that directly convert intercepted solar radiation into thermal energy.

Figure 13 Track-averaged vertical profiles of shortwave (SW) heating rate (K d⁻¹). (a) Mean profiles for the FR and DA experiments. (b) Mean profile of the paired DA-minus-FR heating-rate differences. Both panels use the same CALIOP-track sampling and averaging procedure.

Figure 14 Response of shortwave radiative heating to aerosol extinction perturbations derived from data assimilation. The samples are categorized by SSA into three groups: SSA<0.85 (vermilion circles), $\mathrm{0.85}<\text{SSA}<\mathrm{0.92}$ (green triangles), and SSA>0.92 (blue squares). The solid, dash-dot, and dashed lines represent linear regressions for each category, respectively.

To investigate the physical mechanisms of aerosol radiative effects across different source regions, Fig. 15 presents a specific case study at 06:00 UTC on 7 March 2021, illustrating the latitude-height cross-sections of speciated aerosol extinction and the corresponding shortwave (SW) heating rates. The spatial distribution of aerosols exhibits strong latitudinal gradients due to varying regional emissions and transport patterns. Compared to CALIOP observations, the FR experiment fails to reproduce the vertical distribution of aerosol species, while the DA experiment effectively corrects these vertical structures. The differences between the DA and FR experiments (DA-FR) reveal a strong dependence of the SW heating rate on the distribution of black carbon (BC). For instance, north of 40° N, although the assimilation decreases the dust extinction, it significantly increases the BC concentration. Because BC is highly absorbing, its increase dominates the local radiative response, leading to a net positive SW heating anomaly in this region. Similarly, at lower latitudes (15–23° N), the DA experiment introduces a substantial increase in BC at elevated layers (2–6 km). The enhanced SW heating is spatially consistent with this BC increment, with local instantaneous DA-minus-FR heating-rate increases reaching approximately 3.0 K d⁻¹. The spatial pattern of the enhanced SW heating perfectly mirrors this BC increment, confirming that BC is the primary driver of the mid-level atmospheric warming. Conversely, a prominent negative SW heating anomaly (reduced heating) appears in the lower troposphere (below 2 km) at these lower latitudes (15–23° N). This cooling effect is driven by two mechanisms: first, the dense layer of absorbing BC aloft acts as a radiative shield, significantly attenuating the downward solar radiation that reaches the lower atmosphere; second, the simultaneous reduction of both BC and scattering water-soluble aerosols in the lower layers further diminishes local absorption. Ultimately, this vertical radiative heating structure, strong warming aloft and reduced heating below, can decrease the environmental temperature lapse rate and enhance the static stability of the lower troposphere. This highlights the importance of accurately simulating the vertical distribution of individual aerosol species for radiative forcing assessments.

Figure 15 Latitude-height cross-sections of aerosol component extinctions and SW heating rates along the CALIOP track at 06:00 UTC 7 March 2021. (a–l) Comparison of extinction coefficients (units: km⁻¹) for Dust, BC, and Water-soluble species among CALIOP observations, FR, and DA experiments including the difference (DA minus FR). (n–p) Similar to upper panels but for the SW heating rate (units: K d⁻¹).

5 Conclusion and Discussion

5.1 Conclusion

Accurate quantification of aerosol direct radiative forcing (DRF) remains a significant challenge in climate modeling, largely due to persistent uncertainties in the vertical distribution of distinct aerosol components. Traditional data assimilation (DA) methods relying on column-integrated optical depth or bulk extinction often fail to constrain the specific aerosol composition, leading to misrepresentations of aerosol absorption and scattering properties. To address this under-constrained problem, we developed a component-resolved data assimilation system based on the local ensemble transform Kalman filter (LETKF) with spatial and observational constraints within the WRF-Chem model. By assimilating CALIOP-MODIS synergistic retrievals of species-specific extinction profiles, this system explicitly optimizes the three-dimensional mass distributions of individual aerosol components (dust, sea salt, black carbon, and water-soluble aerosols) over East Asia.

Our evaluations demonstrate that the component-resolved DA experiment substantially improves the spatiotemporal representation of complex aerosol mixtures. The free-running (FR) experiment exhibited systematic underestimations of aerosol loadings and completely missed elevated black carbon (BC) plumes. In contrast, the standard LETKF effectively constrained the simulated vertical structural shapes, while the integrated background correction successfully restored the magnitude of completely missed elevated black carbon (BC) plumes, yielding analysis fields that are highly consistent with the assimilated CALIOP observations. However, we acknowledge two limitations in this evaluation. First, comparing the analysis field against the assimilated CALIOP observations serves as a self-validation, rather than an independent validation. Second, without independent profile data for specific aerosol species, we cannot fully verify the exact distribution of the four individual components. A major advancement of this system is its ability to capture persistent free-tropospheric BC layers over the Indian subcontinent. Such features are frequently missed by chemical transport models due to emission inventory deficiencies and have historically required sparse in-situ balloon or aircraft profiling to detect. Independent validations against ground-based AERONET and AD-Net lidar observations confirmed that the assimilation not only reduced the biases and root mean square errors (RMSEs) of column-integrated AOT, but also accurately captured the vertical evolution of long-range transport episodes. Furthermore, by independently adjusting the mass of scattering and absorbing species, the DA system effectively corrected local biases in the single scattering albedo (SSA), highlighting its unique capability to optimize fundamental aerosol optical properties. The optimized component-specific vertical profiles profoundly modulate the atmospheric shortwave (SW) radiative heating. Our results reveal that the heating efficiency is highly sensitive to the vertical partitioning of absorbing aerosols, particularly BC. The assimilation of elevated BC plumes induced a pronounced mid-tropospheric warming, accompanied by a reduction in lower-tropospheric heating due to the attenuation of downward solar radiation.

5.2 Discussion on the DA system limitations and the bypass mechanism

While the component-resolved assimilation system demonstrates significant improvements, it is crucial to objectively distinguish the contributions of the background correction mechanism from the pure LETKF assimilation. It should be acknowledged that the background correction scheme acts as a simplified nudging approach that bypasses the rigorous uncertainty weighting of the standard Kalman filter. Consequently, this direct rescaling becomes less trustworthy when observational uncertainties are large. Although we implemented a 3-sigma outlier check and required a minimum number of local observations to mitigate the impact of noisy retrievals, ignoring observational uncertainty during the correction step remains a limitation.

To demonstrate the true capability of the LETKF, we have conducted a new assimilation experiment (DA_NoBypass) as pure assimilation. Spatially and temporally, the bypass mechanism was triggered in approximately 28.6 % of the total assimilated grid points across all assimilation cycles. The assimilated grid point is defined as the total number of grid cells containing valid CALIOP observations across all assimilation cycles. Disaggregated by species, it triggered in 43.6 % of cases for dust, 17.8 % for water-soluble aerosols, 24.5 % for BC, and approximately 42.9 % for sea salt. Vertically, the triggers were unevenly distributed across different altitudes, predominantly occurring in the free troposphere. As shown in Fig. 16, the DA_NoBypass experiment significantly improves the CORR across different aerosol species compared to the FR experiment. For instance, the CORR for Black Carbon (BC) dramatically improved from 0.26 (FR) to 0.79 (DA_NoBypass). Similar structural improvements are also observed in sea-salt (CORR improved from 0.87 to 0.93) and water-soluble (from 0.90 to 0.99). The increase in CORR in Fig. 16 indicates improved agreement in the shape of the domain-averaged vertical profiles; it does not, by itself, establish improved temporal or horizontal spatial variability. Point-wise statistics across the full collocated sample are therefore additionally presented in Fig. S3 in the Supplement.

Figure 16 Vertical distribution of extinction coefficients (km⁻¹) for (a) sea salt, (b) dust, (c) black carbon, and (d) water-soluble aerosols. The black solid lines represent the assimilated CALIOP observations. The blue dashed, red solid, and green solid lines denote the results from the FR, DA, and DA_NoBypass experiments, respectively. Statistical metrics (BIAS, RMSE, and CORR) calculated for the vertical profiles are listed in each panel. The profiles represent the domain-averaged values over East Asia (shown in Fig. 1).

However, while DA_NoBypass successfully captures the structural variations, it still underestimates the magnitude of strong aerosol layers compared to the DA experiment. The DA experiment further reduces the RMSE and BIAS, and yields even higher CORR (e.g., reaching 0.95 for BC). Therefore, the extreme local improvements in completely unsimulated plumes highlight a fundamental limitation of ensemble-based DA that empirical corrections remain indispensable when background spread collapses. Conversely, the independent AERONET and AD-Net validations presented in Sect. 4.2 provide a much more robust and conservative estimate of the system's true predictive skill. To verify the spatiotemporal skill of the pure LETKF, point-wise statistics for DA_NoBypass across all valid space-time-height observations are presented in Fig. S3, confirming genuine flow-dependent improvements beyond just the domain-averaged vertical shape. Furthermore, to rigorously isolate the dynamical propagation of the LETKF from the forced local agreement of the bypass, we evaluated the DA_NoBypass experiment against independent AERONET AOT observations (Fig. S4 in the Supplement). The results show that DA_NoBypass still significantly reduces the RMSE and BIAS relative to the FR experiment, confirming that the flow-dependent covariance can improve the downstream aerosol transport, although the DA experiment yields the best overall performance. Because these sites are often far from CALIOP overpass tracks, the modest but solid improvements observed there reflect the genuine dynamical propagation of the assimilated fields rather than the forced local agreement of the bypass mechanism.

5.3 Sensitivity to observational uncertainty

As discussed in Sect. 2.2.1, the 10 % baseline instrumental and retrieval error applied in our DA experiment serves strictly as an empirically tuned parameter, as spatial aggregation cannot simply reduce the native retrieval errors by a factor of $\mathrm{1}/\sqrt{N}$ due to strong spatial correlations. To rigorously evaluate the sensitivity of the assimilation performance to the prescribed observation error covariance, we conducted an additional experiment (DA_Obserr) over the entire study period. In this setup, the instrumental and retrieval error (ϵ_inst+retri) baseline was inflated fivefold, from 10 % to 50 %.

Figure 17 presents a self-validation of the domain-averaged vertical profiles of aerosol extinction coefficients, comparing the FR, DA, and DA_Obserr experiments against the assimilated CALIOP observations. As expected from data assimilation, inflating the observation error variance decreases the Kalman gain, thereby reducing the observational weight during the analysis update. Consequently, the analysis increments in DA_Obserr are slightly constrained, drawing the optimized profiles closer to the model background relative to the DA experiment. For instance, the RMSE for BC slightly increases from 0.031 in DA to 0.033 in DA_Obserr.

Figure 17 Same as Fig. 16 but for DA_Obserr experiment.

Nevertheless, the assimilation system exhibits substantial robustness despite a fivefold increase in the prescribed baseline error. Because the representativeness error (ϵ_rep) frequently dominates the total observation error budget in regions with complex aerosol distributions (as shown in Fig. S1), inflating the instrumental baseline to 50 % does not linearly explode the total observation error fed into the filter. The DA_Obserr experiment successfully reconstructs the complex vertical layering of aerosols that the FR completely fails to capture. The structural improvements remain remarkably consistent; notably, the CORR for BC maintains a high value of 0.94 (compared to 0.95 in DA and 0.26 in FR). This indicates that the primary added value of assimilating speciated vertical profiles lies in providing critical constraints on the vertical distribution and chemical partitioning of aerosols. Even when the observations are assigned lower confidence, the flow-dependent background error covariance effectively propagates the vertical gradient information to mitigate severe structural biases in the model background. It is important to note that the bypass mechanism remained active in the DA_Obserr experiment. Because the background correction criteria (R_corr) do not incorporate the prescribed observation-error covariance, these specific bypass updates are insensitive to the inflated observation error. Therefore, the relatively weak sensitivity of the overall assimilation performance to the 50 % error baseline is attributable not only to the robustness of the flow-dependent covariance and the dominance of representativeness error, but also to the continued operation of this error-insensitive background correction pathway. Therefore, our conclusions regarding the optimization of component-specific vertical structures are highly robust and not overly sensitive to the empirical specification of the observation error baseline.

5.4 Future perspectives

Furthermore, while the component-resolved assimilation system represents an advancement, several physical limitations should be acknowledged. First, the GOCART optical scheme treats aerosols as externally mixed spheres, whereas the assimilated CALIOP-MODIS retrieval assumes internal mixing. This discrepancy introduces representativeness errors in mass extinction efficiency and absorption. Recent studies demonstrate that particle size and mixing state strongly modulate single scattering albedo and radiative forcing, with biases up to −3.0 W m⁻² in industrial regions (Tiwari et al., 2025; Guan et al., 2026a). Importantly, achieving accurate single scattering albedo requires more than mass adjustment alone, particle size and mixing state must also be accounted for (Tiwari et al., 2026), including changes in the spectral signature of absorption across wavebands (Guan et al., 2026b). Second, the elevated BC plumes captured in our analysis are not solely emission-driven, both orographic lofting and local buoyancy-induced factors impact transport from Southeast Asia into the free troposphere play critical roles (Cohen et al., 2017; Wang et al., 2021). Third, as demonstrated by Tiwari et al. (2025), jointly constraining BC mass and number can help distinguish between fresh emission hotspots and aging-dominated regions, while the additional datasets used in that study (TROPOMI UV, AERONET multi-wavelength SSA) are already available over East Asia and could be assimilated to further refine size and mixing-state estimates. Addressing these uncertainties will be essential for reducing biases in aerosol direct radiative forcing.

While this study emphasizes the critical role of component-resolved vertical structures in determining instantaneous radiative heating, the subsequent thermodynamic and dynamic feedbacks need further investigation. Future studies should explore how these component-driven heating rate anomalies perturb atmospheric temperature profiles and regional circulation fields. In particular, the enhanced static stability induced by absorbing aerosols aloft can significantly suppress planetary boundary layer development. Investigating these aerosol-radiation-boundary layer interactions, and their two-way feedback mechanisms on the accumulation of surface pollutants and cloud formation, will be a crucial next step toward a more comprehensive understanding of aerosol-climate interactions. Additionally, future work will explore the assimilation of high-temporal-resolution aerosol component retrievals from ground-based networks, such as AD-Net and SKYNET, to further constrain the diurnal variations of aerosol vertical structures. Furthermore, recent advancements in passive remote sensing, such as the retrieval of black carbon using near-ultraviolet bands from GOSAT-2/CAI-2 imagers (Nakajima et al., 2026), provide unprecedented spatial coverage of absorbing aerosols. Assimilating these novel high-resolution passive retrievals over land areas, particularly over China and India, will be a highly promising direction to further constrain regional aerosol radiative forcing.