The prototype NOAA Aerosol Reanalysis version 1.0: description of the modeling system and its evaluation

In this study, we describe the first prototype version of global aerosol reanalysis at the National Oceanic and Atmospheric Administration (NOAA), the prototype NOAA Aerosol Reanalysis version 1.0 (pNARA v1.0) that was produced for the year 2016. In pNARA v1.0, the forecast model is an early version of the operational Global Ensemble Forecast System-Aerosols (GEFS-Aerosols) model. The three-dimensional ensemble-variational (3D-EnVar) data assimilation (DA) system configuration is built using elements of the Joint Effort for Data Assimilation Integration (JEDI) framework being developed at the Joint Center for Satellite Data Assimilation (JCSDA). The Neural Network Retrievals (NNR) of aerosol optical depth (AOD) at 550 nm from the Moderate Resolution Imaging Spectroradiometer (MODIS) instruments are assimilated to provide reanalysis of aerosol mass mixing ratios. We evaluate pNARA v1.0 against a wide variety of Aerosol Robotic Network (AERONET) observations, against the National Aeronautics and Space Administration's Modern-Era Retrospective Analysis for Research and Applications 2 (MERRA-2; Gelaro et al., 2017; Randles et al., 2017; Buchard et al., 2017) and the European Centre for Medium-Range Weather Forecasts' Copernicus Atmosphere Monitoring Service Reanalysis (CAMSRA; Inness et al., 2019), and against measurements of surface concentrations of particulate matter 2.5 (PM_2.5) and aerosol species. Overall, the 3D-EnVar DA system significantly improves AOD simulations compared with observations, but the assimilation has limited impact on chemical composition and size distributions of aerosols. We also identify deficiencies in the model's representations of aerosol chemistry and their optical properties elucidated from evaluation of pNARA v1.0 against AERONET observations. A comparison of seasonal profiles of aerosol species from pNARA v1.0 with the other two reanalyses exposes significant differences among datasets. These differences reflect uncertainties in simulating aerosols in general.

1 Introduction

Aerosols affect earth's energy balance through the absorption and scattering of solar and terrestrial radiation (Mitchell, 1971). Aerosols also affect the weather and climate through their indirect effects on cloud microphysics, reflectance, and precipitation (Twomey, 1974, 1977; Albrecht, 1989; Jones et al., 1994; Ackerman et al., 2000; Lohmann et al., 2000). In recent decades, the impacts of aerosol direct and indirect radiative effects on numerical weather prediction (NWP) have been studied extensively (Tompkins et al., 2005; Pérez et al., 2006; Mulcahy et al., 2014; Bozzo et al., 2017; Nowottnick et al., 2018). Aerosol particles can also provide the surface area for deposition of gas phase chemicals and subsequently affect the sulfuric oxidation and production of new aerosols (Andreae and Crutzen, 1997). The deposition of mineral dust into the ocean affects marine productivity (Chen et al., 2007; Shao et al., 2011). Aerosols also have fundamental impacts on human health by contributing to respiratory, cardiovascular, and allergic diseases (Pöschl, 2005).

In the past decade, several operational weather prediction centers have produced global aerosol reanalyses following data assimilation (DA) methodologies used for meteorological reanalyses. The National Aeronautics and Space Administration's Global Modeling and Assimilation Office (GMAO) developed an offline (i.e., where the transport is driven using an already reanalyzed dataset) aerosol reanalysis based on the meteorological reanalysis, Modern-Era Retrospective Analysis for Research and Applications (MERRA; Buchard et al., 2015). Subsequently, NASA's GMAO produced the version 2 of MERRA (MERRA-2), which provides both aerosol and meteorological reanalyses by assimilating aerosol and meteorological observations concurrently (Gelaro et al., 2017; Randles et al., 2017; Buchard et al., 2017). The European Centre for Medium-Range Weather Forecasts (ECMWF) produced a global reanalysis of atmospheric composition, including aerosols and trace gases, known as the Copernicus Atmosphere Monitoring Service (CAMS) interim reanalysis (CAMSiRA) (Inness et al., 2013; Flemming et al., 2017), and updated the chemistry and aerosol modules to CAMS reanalysis (CAMSRA) (Inness et al., 2019). The U.S. Navy Research Laboratory (NRL) generated an aerosol reanalysis (Lynch et al., 2016) using the NRL Aerosol Analysis and Prediction System (NAAPS) (Rubin et al., 2016). The Meteorological Research Institute (MRI) of the Japan Meteorological Agency (JMA) also produced a global aerosol reanalysis product named the Japanese Reanalysis for Aerosol (JRAero; Yumimoto et al., 2017).

Aerosol reanalysis represents a uniform and continuous best estimate of the true aerosol state in the atmosphere. It is obtained by constraining model forecasts with observations in the process of data assimilation (Lahoz and Schneider, 2014). Consequently, the reanalysis provides the best available information on temporal and spatial variability of aerosols and can serve as a basis for assessing their impacts on the whole earth system. For instance, Bozzo et al. (2017) incorporated the 11-year aerosol climatology from CAMSiRA to improve the aerosol direct radiative effect in the global forecast model. Benedetti and Vitart (2018) used CAMSiRA data to initialize the aerosol fields in a prognostic aerosol experiment to assess the impact of direct radiative effect on subseasonal forecasts. Bender et al. (2019) included MERRA-2 as a reference dataset to study the aerosol–cloud–radiation interaction. Some other studies have investigated the regional long-term variability of aerosol activities with one or multiple reanalysis datasets (e.g., MERRA-2 and CAMSRA) (Kalita et al., 2020; Cao et al., 2021; Xian et al., 2022). Furthermore, aerosol reanalyses can be used as benchmarks in evaluating a newly developed system (Huang et al., 2023, and this study).

The model is a key component in production of aerosol reanalysis. For the consideration of computational constraints, complex chemistry involving aerosols needs to be reduced by employing parameterizations which only capture the most essential processes. For instance, the Goddard Chemistry Aerosol Radiation and Transport (GOCART) utilized in MERRA-2 assumes no interaction between aerosol species (i.e., external mixing) and considers sulfur oxidation with prescribed climatology for OH, NO₃, and H₂O₂ (Chin et al., 2000, 2002; Colarco et al., 2010). In contrast, CAMSRA (Inness et al., 2019) considers a more comprehensive chemistry, such as the parameterization of the secondary organic aerosols (not considered in CAMSiRA). In NAAPS, the primary and secondary organic aerosols are preprocessed at the initialization stage. The Model of Aerosol Species in the Global Atmosphere (MASINGAR mk-2), which is used to produce JRAero, also simplifies the processes involving production of secondary organic aerosols. Besides the simplification of chemistry, the interaction between aerosols and atmospheric physics is usually neglected. Both MERRA-2 and CAMSRA radiatively coupled the prognostic aerosols but did not consider the cloud–aerosol interaction (Gelaro et al., 2017; Inness et al., 2019). The climate model of JRAero actively coupled the aerosol information from MASINGAR mk-2 to the radiation and the two-moment bulk cloud microphysics schemes (Yukimoto et al., 2012).

Assimilated observations are a critical component of a reanalysis product. For aerosols, these are typically retrievals of aerosol optical depth (AOD) at 550 nm. Because AOD is a measure of extinction of light over a whole atmospheric column, it carries a limited amount of information on the chemical composition and sizes of particles and also cannot provide information on the vertical distribution in the atmosphere. Consequently, the resulting analysis is largely determined by the prior information contained in the model state at a given time. This provides very limited ability for AOD DA to correct aerosol speciation, size, and vertical distributions in the atmosphere.

In 2011, the National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Prediction (NCEP) implemented the NOAA Environmental Modeling System (NEMS) Global Forecast System (GFS) Aerosol Component, version 1 (NGAC v1; Lu et al., 2016) with GOCART parameterization for dust-only forecasts. Later, NCEP upgraded the system to NGAC v2, which provided forecasts for all aerosol species in GOCART parameterization (Wang et al., 2018; Bhattacharjee et al., 2018). In September 2020, NGAC was replaced by the Global Ensemble Forecast System (GEFS)-Aerosols model which relies on the Finite-Volume Cubed-Sphere (FV3) dynamical core and GOCART parameterization (Zhang et al., 2022). Currently, the aerosol forecasting system at NCEP does not include DA, despite previous studies showing significant impact of AOD observations on aerosol initial conditions and the subsequent forecasts (Benedetti et al., 2009; Pagowski et al., 2010, 2014; Schwartz et al., 2014; Choi et al., 2020). A novel three-dimensional ensemble-variational (3D-EnVar) DA system based on elements of the Joint Effort for Data Assimilation Integration (JEDI) has been developed at NOAA to provide the best estimates of atmospheric aerosol distributions by assimilating AOD at 550 nm from the Visible Infrared Imaging Radiometer Suite (VIIRS) (Huang et al., 2023).

In this study, we used and designed a specific JEDI-based 3D-EnVar DA configuration to produce the prototype NOAA Aerosol Reanalysis version 1.0 (pNARA v1.0). The reanalysis for 2016 was produced as a pilot product and its evaluation is presented here. We evaluated the year-long reanalysis against independent observations and the gridded aerosol products from MERRA-2 and CAMSRA. This paper is structured as follows: Section 2 describes the observations used in this study. Sections 3 and 4 describe the forward model and the DA system, respectively. Section 5 reports the results of the evaluations. Section 6 concludes the study and provides perspectives for future work.

2 Observations

pNARA v1.0 assimilates the Neural Network Retrieval (NNR) AOD derived from observations of the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard Terra and Aqua satellites (Castellanos and da Silva, 2017; Randles et al., 2017). Unlike the standard retrieval method that relies on a physical model (e.g., a radiative transfer model), the NNR approach retrieves the AOD based on a well-trained model conducted by a neural network machine learning method. The training dataset contains MODIS level-2 reflectances, glint, solar and sensor angles, cloud fraction, and albedo derived using GEOS-5 surface wind speeds or climatology over ocean or land, respectively, targeting AOD observations from the Aerosol Robotic Network (AERONET). The NNR AOD is available at wavelengths over ocean: 470, 500, 550, 660, and 870 nm; over dark land: 440, 470, 500, 550, 660, and 870 nm; and over bright land: 440, 470, 500, 550, and 660 nm. Only MODIS NNR AOD at 550 nm has been assimilated in MERRA-2 (Randles et al., 2017) and pNARA v1.0 (this study).

For evaluation of pNARA v1.0 and the free model run, we used the version 3 level-2.0 AOD retrievals from AERONET (Giles et al., 2019), particulate matter 2.5 (PM_2.5) measurements from the OpenAQ dataset (https://openaq.org; last access: 3 May 2023), and speciated aerosol surface concentration measurements from the Interagency Monitoring of Protected Visual Environments (IMPROVE) network (Hand et al., 2019; http://vista.cira.colostate.edu/Improve/improve-program/; last access: 3 May 2023). AERONET is a ground-based remote sensing aerosol network that uses a sun photometer to measure direct solar irradiance and provides the AOD retrievals at 340, 380, 440, 500, 670, 870, 940, and 1020 nm. In the past more than 25 years, the network has expanded to more than 600 stations. The manual quality control used in version 2 leads to a significant delay in the level-2.0 database. AERONET version 3 fully automated the quality control algorithm to screen out cloud contaminated data, which eliminates the manual efforts and reduces the processing time for quality-assured data (i.e., level 2.0).

OpenAQ is an open-source platform collecting air quality measurements from various sources globally. In this study, we only use government- and research-institute-maintained (i.e., “reference grade”) stations, which are majorly over the USA and Europe, to evaluate the performance of pNARA v1.0. Besides PM_2.5 measurements, the OpenAQ database also ingests the measurements of gas phase pollutants, such as O₃, NO₂, and SO₂, which can be utilized as an independent dataset for trace gas evaluation in future studies. The IMPROVE network was established as the visibility network over the USA in 1985. It derives visibility metrics through the measurement of speciated aerosol mass concentrations, including anions of sulfate, nitrate, nitrite, chloride, and organic and elemental carbonaceous species.

3 Model description

In pNARA v1.0, we use an early Common Community Physics Package (CCPP) version of the operational Global Ensemble Forecast System-Aerosols (GEFS-Aerosols) model. GEFS-Aerosols provides prognostic mass mixing ratios of five bins dust, five bins sea salt, hydrophobic and hydrophilic black and organic carbon, and sulfate. It relies on GOCART parameterization implemented with the GFS physics package and Finite-Volume Cubed-Sphere (FV3) dynamical core (Lin, 2004; Harris et al., 2021). It handles sub-grid transport and wet scavenging of aerosols in the atmospheric physics module. Other chemical processes, such as emissions, chemical reactions, dry deposition, and settling, are handled by the chemical module driven by meteorological fields from the atmospheric component of the model. Note that the formation of secondary organic aerosols and nitrates are not considered in the GEFS-Aerosols. The model has been recently updated for operations to reduce some biases (Zhang et al., 2022), but the most recent version was not available at the time of the execution of this study.

In GEFS-Aerosols, the background fields of OH, H₂O₂, and NO₃ used in the parameterization of simplified sulfur chemistry in GOCART are updated with a monthly mean climatology from the 2015 version of NASA Global Modeling Initiative's (GMI) chemical model. The model also includes the 1D plume rise module adapted from the High-Resolution Rapid Refresh (HRRR) Smoke model to improve vertical distribution of smoke emissions. The biomass burning emissions use version 3 of Blended Global Biomass Burning Emissions Product (GBBEPx v3; Zhang et al., 2019). This parameterization blends the Quick Fire Emissions Dataset (QFED) used in MERRA-2 (Darmenov and da Silva, 2015) and daily emissions derived by the hotspots and the fire radiative power observations from polar-orbiting and geostationary satellites. The anthropogenic emissions are based on the inventories from the Community Emissions Data System (CEDS). The sea salt scheme has been updated to the recent version of the GOCART scheme (Colarco et al., 2010). The dust scheme has been updated to FENGSHA (Dong et al., 2016; Zhang et al., 2022), which means the “wind-blown dust” in Mandarin Chinese.

4 Data assimilation system

To produce pNARA v1.0, we used a three-dimensional ensemble-variational (3D-EnVar) aerosol data assimilation system using components from the Joint Effort for Data Assimilation Integration (Huang et al., 2023). JEDI is primarily developed at the Joint Center for Satellite Data Assimilation (JCSDA) but receives significant contributions from partner and sponsor agencies (i.e., NOAA, NASA, U.S. Navy, U.S. Air Force, and U.K. Met Office). It aims to provide an integrated and unified DA framework for earth system prediction applications and reduce redundant efforts across research and operational communities. (More information about JEDI can be found on the JCSDA website (https://www.jcsda.org/jcsda-project-jedi; last access: 3 May 2023.)

Figure 1 Schematic of the 3D-EnVar data assimilation system for pNARA v1.0. Note that IODA stands for Interface for Observational Data Access and UFO stands for Unified Forward Operator.

The forward observation operator that converts model variables (mass mixing ratios of aerosols and pressure layer depths) to AOD observations relies on aerosol specific extinction coefficients interpolated from look-up tables (i.e., “AodLUTs” operator in JEDI). Tabulated values were obtained at NASA GMAO from the Mie theory for spherical particles (Wiscombe, 1980) and from the T-matrix approach for non-spherical dust (Meng et al., 2010).

In this application, the 3D-EnVar aerosol DA system assimilates MODIS NNR AOD retrievals at 550 nm to the GEFS-Aerosol model every 6 h. For MODIS NNR AOD, the cloud-affected data was screened and trained with AERONET measurements during the neural network training (Randles et al., 2017). We also thinned the horizontal resolution from 10 to 50 km to reduce the correlation between observations. The result is a reanalysis of mass mixing ratios of five GOCART aerosol species, with sea salt and dust distributed over five size bins. Figure 1 displays the 3D-EnVar aerosol DA system. The 3D-EnVar system used for pNARA v1.0 is a combination of a local ensemble transform Kalman filter (LETKF; Bishop et al., 2001; Hunt et al., 2007) using 40 members and a 3D-variational (3D-Var; Lorenc, 1986) solver that produces the control analysis using the ensemble information to obtain the background error covariance matrix. The LETKF analyses are re-centered on the 3D-Var analysis to ensure consistency between the two analyses (Hamill and Snyder, 2000; Lorenc, 2003; Buehner, 2005). The control analysis constitutes the reanalysis output. The forecasts and analyses of control and ensemble members are conducted at C96 (∼100 km) resolution. Given the coarse resolution, this study focuses on the larger temporal and spatial scale performance. Note that no other observations are assimilated since meteorological fields in the control and in the ensemble are driven by analyses from the global data assimilation system (GDAS) at NCEP.

To address model biases and spread deficiency of the ensemble, a scheme to scale and perturb source emissions was devised. Scaling factors are derived for dominant aerosols based on AOD deficits over regions. The scale factor of 2 is selected for all aerosol species to increase the emission rate for reducing the model biases. Perturbations to emissions of ensemble members represent spatially and temporally correlated patterns following the approach in the stochastically perturbed parametrization tendencies (SPPT; Palmer et al., 2009) scheme. With applying SPPT to emissions, it can enhance the variance between ensemble members (i.e., larger ensemble spread), and thus the influences of observations to the analysis are increased. This approach and a method to obtain scaling factors and amplitudes of emission perturbations are detailed in Huang et al. (2023).

Figure 2 Probability density plots of modeled AOD at 500 nm with respect to AERONET retrievals for the free model run and pNARA v1.0 during four seasons. Axes are in logarithmic scale.

5 Evaluation

In this section, we divide the year 2016 into four seasons to investigate the performance of pNARA v1.0. Winter includes December, January, and February (denoted as DJF); spring includes March, April, and May (denoted as MAM); summer includes June, July, and August (denoted as JJA); and fall includes September, October, and November (denoted as SON).

5.1 Comparison against AERONET

Figure 2 illustrates the global seasonal comparison of AOD at 500 nm from pNARA v1.0 and the free model run with respect to AERONET, and the statistics are also provided. This figure lists absolute and relative biases, R² correlations (also known as coefficients of determination), and, for the reader's convenience, correlation coefficients (R) for each season and the whole year for the free model run and pNARA v1.0. It is followed by Fig. 3 where time series display the bias and the R² correlation scores. Because the standard AOD at 550 nm is not directly observed in AERONET, we chose to perform the evaluation of AOD at 500 nm to avoid interpolation and log-linearization errors in AERONET observations. The free model run considerably and systematically underestimates the AOD throughout the year. Among the four seasons, the largest bias between model and measurements occurs in winter (i.e., DJF). Compared with the free model run, our reanalysis, pNARA v1.0, has a substantially better agreement with AERONET throughout the whole year in terms of bias and correlation. We note that the statistical scores vary seasonally as they are influenced by the density of AERONET observations that are used in the evaluation. (The highest and lowest density of observations occurs in the summer and winter, respectively, over North America and Europe.)

Figure 3 The time series of the monthly mean of (a) biases and (b) R² correlation for the free model run (in blue) and pNARA v1.0 (in red) with respect to AOD at 500 nm measurements from AERONET. Both are calculated based on all the paired points within each month.

Figure 4 Probability density plots of 440–870 nm Ångström exponent of the free model run (a) and pNARA v1.0 (b) vs. AERONET for 2016.

As can be seen in Fig. 3, the performance of the free model run in terms of bias is the best during summer months and the worst from January to March. For pNARA v1.0, the negative absolute biases are much reduced (on average by about 0.05), though not enough to remove the overall negative bias of the reanalysis. pNARA v1.0 also has significantly improved R² correlation throughout the whole year (by about 0.2). Statistics of the free model run and reanalysis vary throughout the year nearly in parallel demonstrating the crucial role of the forecast model in data assimilation. In other words, data assimilation cannot drastically improve the quality of simulations if the forecast model is seriously deficient.

Figure 4 displays the probability density plots of 440–870 nm Ångström exponent (AE) for the free model run and the reanalysis against the AERONET observations throughout 2016. The AE is calculated by Eq. (1). It provides a measure of relative extinction of light by aerosols at different wavelengths and primarily reflects the size distribution of particles (Schuster et al., 2006):

$$\begin{array}{}\text{(1)}& \mathrm{AE}\left({\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}}\right)=-{\displaystyle \frac{\log \left(\mathrm{AOD}\left({\mathit{\lambda }}_{\mathrm{1}}\right)/\mathrm{AOD}\left({\mathit{\lambda }}_{\mathrm{2}}\right)\right)}{\log \left({\mathit{\lambda }}_{\mathrm{1}}/{\mathit{\lambda }}_{\mathrm{2}}\right)}},\end{array}$$

where λ₁ and λ₂ are 440 and 870 nm, respectively.

In general, the differences between the free model run and the reanalysis are insubstantial, and both correlate poorly with observations. The lack of an improvement in the reanalysis demonstrates that assimilating AOD at 550 nm alone only minimally impacts size distribution and/or the composition of aerosols. Saide et al. (2013) and Tsikerdekis et al. (2021) have shown that a more realistic representation of aerosols can be produced by assimilating multi-wavelength retrievals of AOD, fine-mode fraction AOD, and single scattering albedo. We hope more aerosol retrieval products from future missions, such as the Plankton, Aerosol, Cloud, Ocean Ecosystem (PACE) mission (https://pace.gsfc.nasa.gov; last access: 3 May 2023), can introduce more information into the analysis system in the near future.

Figure 5 Probability density plots of 440–675 nm absorption Ångström exponent vs. scattering Ångström exponent for AERONET (a) and pNARA v1.0 (b) for 2016. OC: organic carbon; EC: elemental carbon.

In Fig. 5, probability density plots of 440–675 nm absorption AE (AAE) versus Scattering AE (SAE) are matched for AERONET Almucantar retrievals (Sinyuk et al., 2020) and pNARA v1.0. For brevity, a similar scatter plot for the free model run is omitted since it is only marginally different from the latter. We obtained the model's AAE and SAE with Eqs. (2a) and (2b):

$$\begin{array}{}\text{(2a)}& {\displaystyle }\mathrm{AAE}\left({\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}}\right)=-{\displaystyle \frac{\log \left(\mathrm{AAOD}\left({\mathit{\lambda }}_{\mathrm{1}}\right)/\mathrm{AAOD}\left({\mathit{\lambda }}_{\mathrm{2}}\right)\right)}{\log \left({\mathit{\lambda }}_{\mathrm{1}}/{\mathit{\lambda }}_{\mathrm{2}}\right)}},\text{(2b)}& {\displaystyle }\mathrm{SAE}\left({\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}}\right)=-{\displaystyle \frac{\log \left(\mathrm{SAOD}\left({\mathit{\lambda }}_{\mathrm{1}}\right)/\mathrm{SAOD}\left({\mathit{\lambda }}_{\mathrm{2}}\right)\right)}{\log \left({\mathit{\lambda }}_{\mathrm{1}}/{\mathit{\lambda }}_{\mathrm{2}}\right)}},\end{array}$$

where absorption AOD (AAOD) and scattering AOD (SAOD) at wavelengths λ₁ (equal to 440 nm) and λ₂ (equal to 675 nm) are calculated using the Eqs. (3a) and (3b). The calculations use single scattering albedo (SSA) and extinction coefficient (Ext) at corresponding wavelengths (λ), mixing ratios (q) for each aerosol species (s), dry air density (ρ), and layer thickness (dz) for each pressure level (p):

$$\begin{array}{}\text{(3a)}& {\displaystyle }\begin{array}{rl}\mathrm{AAOD}\left(\mathit{\lambda }\right)& ={\sum }_p\left[{\sum }_s(\mathrm{1}-{\mathrm{SSA}}_{s,\mathit{\lambda }})\times {\mathrm{Ext}}_{s,\mathit{\lambda }}\times q_{s,p}\right]\\ & \times \mathit{\rho }\left(p\right)\times \mathrm{d}z\left(p\right)\end{array}\text{(3b)}& {\displaystyle }\begin{array}{rl}\mathrm{SAOD}\left(\mathit{\lambda }\right)& ={\sum }_p\left[{\sum }_s{\mathrm{SSA}}_{s,\mathit{\lambda }}\times {\mathrm{Ext}}_{s,\mathit{\lambda }}\times q_{s,p}\right]\\ & \times \mathit{\rho }\left(p\right)\times \mathrm{d}z\left(p\right).\end{array}\end{array}$$

Based on Cazorla et al. (2013), the whole plane in the left panel of Fig. 5 was further divided into separate sections representing different types of dominant aerosol species. These authors associate AAE with the representation of chemical composition and SAE with the representation of particle sizes. Figure 5 reveals important deficiencies in the representation of optical properties of aerosols in the model. Because of the dependence of absorption on refractive indices of aerosols, we believe that the strong lack of variability in AAE (vertical axis in Fig. 5) in pNARA v1.0 compared with AERONET is explained by assumptions and simplifications in GOCART parameterization. These are, in the order of importance, the external mixing of the aerosols in the model, uniform mineralogy of dust sources across the globe, unaccounted variety of organic aerosols, and prescribed size distributions of particles. A separate investigation would be required to ascertain the impacts of the above factors on the realism of simulations. Given the importance of radiation absorption by atmospheric aerosols in stratification and clouds, our exposure of the GOCART shortcomings identified above highlight uncertainties in modeling meteorology–chemistry interactions in weather and climate models that rely on this or similar schemes.

Figure 6 Comparison of seasonal mean biases of AOD at 550 nm of pNARA v1.0 vs. CAMSRA and MERRA-2.

5.2 Comparison against MERRA-2 and CAMSRA reanalyses

In this section, the AOD and aerosol vertical distribution in pNARA v1.0 are evaluated against MERRA-2 and CAMSRA. To do the direct comparison between these three datasets, we post-processed the aerosol mass mixing ratio in pNARA v1.0 and MERRA-2 from its model levels to the same pressure levels as CAMSRA based on the level thickness and surface pressure. For any levels below the terrain, the surface value is assigned, which is the same approach as in CAMSRA.

Figure 6 shows the seasonal mean biases of AOD at 550 nm of pNARA v1.0 with respect to CAMSRA and MERRA-2. Overall, pNARA v1.0 is closer to MERRA-2. This can be attributed to the similarity of aerosol modules (i.e., GOCART) and the same FV3 dynamical cores (Lin, 2004; Harris et al., 2021) in the models used in NOAA and NASA reanalyses. However, pNARA v1.0 generates substantially smaller dust plumes over the Sahara overall and especially during summer. Compared with CAMSRA, pNARA v1.0 shows significant discrepancies over oceans throughout 2016. It also has lower aerosol concentrations over India, especially during the winter. The AOD in pNARA v1.0 is usually lower over eastern China with the largest differences occurring in the spring. These instances of lower aerosol concentrations could be attributed to the fundamental differences of aerosol modules between pNARA v1.0 and CAMSRA. Throughout 2016 pNARA v1.0 generated more dust over the Gobi Desert compared with both CAMSRA and MERRA-2. Occasionally, concentrations of dust over this area in pNARA v1.0 seem unrealistically high. This suggests the need for further improvements to the model parameterization of dust uplift, which is determined by meteorological and land surface conditions.

Figure 7 The areas selected for comparison of vertical profiles of aerosols.

Figure 8 Total dust mass mixing ratio profiles averaged over the North Atlantic Ocean (NATL) (a) and North Africa and the Middle East (NAFRME) (b) during JJA.

While the size distributions of dust and sea salt are identical in MERRA-2 and pNARA v1.0, they are different from those in CAMSRA. Approximate conversions between the models would be possible, but not straightforward, since the bin sizes and distribution parameters for these two aerosol types differ quite significantly (e.g., dust bins in MERRA-2: 0.1, 1.0, 1.8, 3.0, 6.0, and 10.0 µm, and in CAMSRA: 0.03, 0.55, 0.9, and 20.0 µm; sea salt bins in MERRA-2: 0.03, 0.1, 0.5, 1.5, 5.0, and 10.0 µm in dry conditions, and in CAMSRA: 0.03, 0.5, 5.0, and 20.0 µm at relative humidity 80 %). Therefore, in the following we chose to compare the total mass mixing ratios of these aerosol species for all the reanalyses. Comparison of vertical profiles of various aerosols was performed over several geographic areas where different aerosols are expected to dominate (Fig. 7). Here, for illustration purposes only, we compare the dust aerosols over North Africa and the Middle East (NAFRME) and the North Atlantic Ocean (NATL; transported dust), the carbonaceous aerosols over equatorial Africa and South Africa and the surrounding tropical Ocean (SAFRTROP) and Siberia (RUSC2S), the anthropogenic aerosols over East Asia (EASIA), and the sea salt aerosols over the Southern Ocean (SOCEAN).

Figure 8 displays the vertical profiles of total dust mass mixing ratios from the free model run, pNARA v1.0, MERRA-2, and CAMSRA, over NAFRME and NATL during the summer. In general, the four datasets show variable discrepancies in dust depending on the region of interest. Over the NATL region the dust profiles are similar to each other, while over the NAFRME region GOCART-based models generate higher dust mass mixing ratios near the surface (in particular with the MERRA-2 dataset). Where the Saharan dust profiles are uplifted over the NAFRME region the free model run creates a more vertically mixed profile, while pNARA v1.0 and MERRA-2 show higher values near the surface. By assimilating the similar observation dataset (i.e., MODIS NNR), the 3D-EnVar system corrects the analysis toward MERRA-2.

Figure 9 Same as Fig. 8, but for hydrophilic organic and black carbon (OCPHILIC and BCPHILIC) and hydrophobic organic and black carbon (OCPHOBIC and BCPHOBIC) aerosols over South Africa and the surrounding tropical Ocean (SAFRTROP; top) and Siberia (RUSC2S; bottom) during JJA.

.

Figure 10 Same as Fig. 8, but for sulfate (a, c) and total dust (b, d) aerosols over East Asia (EASIA) during DJF (a, b) and MAM (c, d)

Figure 9 shows a comparison of profiles of carbonaceous aerosols over SAFRTROP and RUSC2S during JJA. As illustrated in Fig. 6, high values of AOD systematically occur over these areas in 2016. For CAMSRA, concentrations of all carbonaceous aerosols are significantly higher near the surface when compared with the other simulations. Interestingly, much higher concentrations of hydrophilic organic and black carbons (OC and BC) exist for CAMSRA at mid-levels. A comparison between the free model run and pNARA v1.0 shows that DA introduced more carbonaceous aerosols over RUSC2S, while there are no considerable differences over SAFRTROP. Comparison of all these simulations indicates that there exist notable discrepancies between vertical profiles of carbonaceous aerosols over areas where extensive wildfires occurred. These discrepancies can be attributed to different parameterizations of biomass burning emissions in the models. For instance, GEFS-Aerosols utilized the GBBEPx biomass burning emissions (Zhang et al., 2019), which is conducted based on QFED (Darmenov and da Silva, 2015) and the fire radiative power observations from polar-orbiting and geostationary satellites. CAMSRA used biomass burning emissions from the Global Fire Assimilation System (GFAS; Kaiser et al., 2012), which shows lower emission compared with QFED (Pan et al., 2020). Also, the smoke plume rise is considered in GEFS-Aerosols but not in MERRA-2 and CAMSRA. These aerosols attest to large uncertainties that exist in parameterizing wildfires.

Figure 10 illustrates vertical profiles of sulfate and dust aerosols over East Asia (EASIA) during DJF and MAM. The results are consistent with the negative AOD biases of pNARA v1.0 over EASIA shown in Fig. 6. The AOD deficit over EASIA in pNARA v1.0 compared with the other reanalyses occurs as a consequence of considerably lower concentrations of both sulfate and dust aerosols.

Figure 11 Same as Fig. 8, but for sea salt, dust, sulfate, and hydrophobic OC over the Southern Ocean during JJA (top) and SON (bottom).

Figure 12 Probability density plots of modeled PM_2.5 with respect to OpenAQ measurements for the free model run (a) and pNARA v1.0 (b). Evaluation period: May–December 2016.

In Fig. 6, pNARA v1.0 displays significant negative AOD bias over oceans with respect to CAMSRA and smaller bias over oceans with respect to MERRA-2. We chose the area marked as SOCEAN in Fig. 7 to investigate reasons for the differences between the reanalyses. In Fig. 11, we show vertical profiles of mixing ratios of selected aerosol species over this area during JJA and SON. Among the three reanalyses sea salt aerosols were lifted to higher elevation in CAMSRA resulting in larger AOD values over the oceans. AOD assimilation in pNARA v1.0 leads to a better agreement with CAMSRA and MERRA-2 compared with the free model run, but the sea salt loading is still lower during JJA. It is worth mentioning that there exist divergences between reanalyses in profiles of mass mixing ratios of other aerosol species, though the values are small. MERRA-2 has the most dust while CAMSRA has the most hydrophobic OC. Both MERRA-2 and CAMSRA have more sulfate aerosols than pNARA v1.0. Shapes of vertical profiles of mass mixing ratios of dust and sulfate vary widely between the reanalyses.

In conclusion, we note that assimilation of AOD leads to convergence of pNARA v1.0 towards both MERRA-2 and CAMSRA. This statement holds for AOD and generally for vertical profiles of aerosols. However, assimilation of an integrated quantity, such as AOD over a single wavelength (here 550 nm), has little impact on the “shape” of vertical profiles, and values of mass mixing ratios usually appear to be scaled proportionally throughout the depth of the atmosphere. Most importantly, there exist marked differences in seasonal vertical profiles of aerosol species between the three reanalyses.

5.3 Comparison against in situ measurements of surface aerosol concentrations

In the following we present an evaluation of pNARA v1.0 against measurements of surface PM_2.5 concentrations from OpenAQ and concentrations of aerosol species from IMPROVE.

Figure 13 Probability density plots of speciated aerosol surface concentration of pNARA v1.0 with respect to IMPROVE measurements for (a) elemental carbon (EC), (b) organic carbon (OC), (c) sea salt, and (d) sulfate.

Hourly measurements of surface concentrations of PM_2.5 collected in the global OpenAQ database are available from late April 2016 onwards. First, we note that in situ measurements display high spatial and temporal variability, which reflects the origin and subsequent evolution of this species. The horizontal resolution of our model (about 100 km), which also affects accuracy of the representation of terrain topography, is far too coarse for the results to be compatible with such measurements. Also, these measurements are obtained close to the ground, which is significantly lower than the bottom level of our model (on average about 20 m). This makes the calculation of the concentration at the surface dependent on the atmospheric stratification in the surface layer. Because of uncertainties involved, here we simply multiply the mixing ratio of model PM_2.5 expressed in Eq. (4) by moist air density obtained from the model diagnostic 2 m temperature, humidity, and surface pressure:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathrm{PM}}_{\mathrm{2.5}}& =\mathrm{BC}\mathrm{1}+\mathrm{BC}\mathrm{2}+\mathrm{OC}\mathrm{1}+\mathrm{OC}\mathrm{2}+\mathrm{SULF}+\mathrm{DUST}\mathrm{1}\\ & +\mathrm{DUST}\mathrm{2}\times \mathrm{0.38}+\mathrm{SEAS}\mathrm{1}+\mathrm{SEAS}\mathrm{2}+\mathrm{SEAS}\mathrm{3}\\ & \times \mathrm{0.83}.\end{array}\end{array}$$

In this equation acronyms denote the following aerosol species: BC – hydrophilic and hydrophobic black carbonaceous, OC – hydrophilic and hydrophobic organic carbonaceous, SULF – sulfate, DUST – dust in the two smallest size bins, and SEAS – sea salt in the three smallest size bins. Finally, we note that AOD, which represents a column-integrated aerosol quantity, may poorly correlate with surface values alone since aerosol-rich layers often occur at raised elevations. Nevertheless, in our opinion, such evaluation statistics not only provide valuable information on the reliability of global model results for PM_2.5 forecasting and analysis but also on the relevance of assimilating AOD for such purposes. The statistics were calculated for 00:00, 06:00, 12:00, and 18:00 UTC and are presented in probability density plots in Fig. 12. Model bias with respect to in situ PM_2.5 measurements coincides with its bias against in situ AERONET AOD retrievals. The correlation with PM_2.5 measurements is much poorer compared with the correlation with AERONET AOD retrievals. Overall and modest improvements in statistics can be noted for the reanalysis compared with the model free run. The spatial density of OpenAQ measurements varies considerably over the globe and it is the highest over North America, Europe, and East Asia. To account for the geographical variability, we performed tests with spatial thinning of data, but sensitivity of the statistics to this procedure was minor.

Measurements of daily averaged (beginning at midnight local standard time) concentrations of selected aerosol species by the IMPROVE network occur once in 3 days and are limited to North America and a single site in South Korea. The sites are located remotely from industrial and population centers. For comparison with the model following guidance from Hand at al. (2019), measured concentrations of organic carbonaceous species and sulfate were scaled by 1.8 and 1.375, respectively. Six-hourly model mixing ratios of species were weighted appropriately to account for the local time of the measurements and, as for comparison with PM_2.5 from OpenAQ, multiplied by 2 m moist air density. Probability density plots for the carbonaceous species, sea salt, and sulfate for the reanalyses, which differ little from the free model run, appear in Fig. 13. Poor performance of the model and lack of improvement from the assimilation can be noted. Because of the scarcity of systematic measurements of concentrations of individual aerosol species, or lack thereof, outside of North America, our statistics only reflect performance of the model over this geographical area. It would be, however, overly optimistic to expect that the model skill significantly improve elsewhere.

6 Conclusions

This study documents the development and the evaluation of pNARA v1.0, which is the first global aerosol reanalysis product at NOAA. pNARA v1.0 is made available for distribution at https://esrl.noaa.gov/gsd/thredds/catalog/retro/global_aerosol_reanalysis/catalog.html (last access: 3 May 2023). To produce pNARA v1.0, the GEFS-Aerosols model was used to forecast aerosol mass mixing ratios. The GEFS-Aerosols relies on GOCART parameterization coupled with the FV3 dynamical core using GFS physics. The reanalysis containing 3D aerosol mass mixing ratios for 2016 was generated in a process of a 3D-EnVar data assimilation using a specifically designed JEDI configuration. The system assimilated neural network retrievals of AOD at 550 nm from MODIS instruments onboard Terra and Aqua satellites. We evaluated the AOD against AERONET observations, the reanalyses from NASA GMAO's MERRA-2, and ECMWF's CAMSRA, and measurements of surface concentrations of PM_2.5 and selected aerosol species.

The evaluation against AERONET observations shows that the assimilation of AOD retrievals at 550 nm improves the overall agreement with AOD at 500 nm, especially in the spring and summer. However, the comparison of the Ångström exponent indicates that the assimilation of the single-wavelength AOD induces minimal improvements to the speciation and the size distributions of aerosols. Furthermore, the significant underdispersion in the probability density plot of the absorption AE vs. the scattering AE along the ordinate (Fig. 5) indicates what we believe is a shortcoming of the GOCART parameterization. We suspect that this shortcoming stems primarily from the assumption of the external mixing of aerosols (i.e., no interaction among aerosol species) in this scheme.

In terms of AOD at 550 nm, pNARA v1.0 shows close proximity to MERRA-2, while it has significant negative biases against CAMSRA. pNARA v1.0 shows lower amounts of dust aerosol over the Sahara and Middle East. The GEFS-Aerosols model and consequently our reanalysis display higher dust concentrations over the Gobi Desert than both reanalyses. These deficiencies over deserts call for improvement to model parameterization of dust uplift. For biomass burning areas, GOCART-based reanalyses, MERRA-2, and pNARA v1.0 share similar characteristics, while CAMSRA shows significantly different amounts and ratios of carbonaceous aerosol species. Over East Asia, pNARA v1.0 displays significantly lower amounts of aerosols than CAMSRA and MERRA-2. This is due to the lower concentrations of sulfate aerosols throughout the year and dust during winter and spring. The discrepancies between the reanalyses indicate that there exist significant differences in parameterizations of anthropogenic, biomass burning, and wind-driven dust uplift between the three models.

In this study, we demonstrated that the reanalyses of AOD at 550 nm produced at NASA (Buchard et al., 2017), ECMWF (Inness et al., 2019), and NOAA (this study) show marked differences globally. (See Table A1 for further comparison among reanalyses.) Also, seasonally averaged vertical profiles of aerosol species vary significantly between the three reanalyses. This pertains both to the chemical composition of aerosol mixtures and the shapes of the vertical profiles of species. The differences between the reanalyses may arise from distinct model numerics, physical and chemical parameterizations, boundary conditions as sources of tracers, and approaches in data assimilation. These observations combined with the poor representation of the absorption spectrum (Fig. 5) by simple aerosol schemes, such as GOCART, lead us to the conclusion that our assessment of the state of atmospheric aerosols and their radiative impacts is hardly adequate to allow detailed forecasts of stratification and clouds with aerosol-sensitive physical parameterizations. For instance, Bozzo et al. (2017) and Mulcahy et al. (2014) have demonstrated the importance of including more realistic aerosol states in the parameterization of aerosol–cloud–radiation interactions. For reasons noted above, we believe that at the current state of aerosol science, large uncertainties in simulating aerosol–meteorology interactions in the weather and climate models exist.

Comparisons with measurements of surface concentrations of PM_2.5 from OpenAQ show limited skill of the model and reanalysis in this task. (See Table A2 for further comparison among reanalyses.) Performance of the model and reanalysis were poor when the results were compared with surface concentrations of carbonaceous, sulfate, and sea salt aerosols. Given the limitations of the model that were listed in the section dedicated to its evaluation against measured concentrations of aerosols at the surface, such results may not be unexpected but are nevertheless disappointing since they suggest that our global reanalysis has a limited value for those health and epidemiological studies in which chemical composition of aerosols is considered. It is possible that extending chemical parameterizations to include the formation of secondary organic aerosols would improve model forecasts (Fan et al., 2022).

Future space missions (e.g., PACE which we mentioned above) and new algorithms (e.g., Zhou et al., 2021) promise to provide novel retrievals of multi-wavelength solar and lunar AOD, fine-fraction AOD, single scattering albedo, and retrievals of aerosol layer height. These developments should enhance scope of evaluations, pose additional constraints in data assimilation, and eventually lead to better aerosol forecasts and reanalyses.

The main goal of this paper is to present our initiative to produce the first-ever global aerosol reanalysis at NOAA. As discussed above, our reanalysis has deficiencies that will be addressed in turn. Huang et al. (2023) outlined the potential enhancements to our assimilation approach, estimates of observation and model errors, and systemic correction of model biases. We will report our advances in the future.

Appendix A

Here we provide the spatial distribution of the temporal statistics from the experiments (Figs. A1 and A2) and the comparison of the three reanalysis datasets (Tables A1 and A2).

Figures A1 and A2 display the temporal biases, R², and RMSE of AOD at 500 nm against the AERONET retrievals (Fig. A1) and the PM_2.5 measurements against OpenAQ (Fig. A2) from both the free model run and pNARA v1.0 analyses. For AOD at 500 nm, pNARA v1.0 shows smaller biases, better correlation, and reduced RMSE in most places. However, over Asian regions (e.g., India and Southeast Asia), the improved correlation and similar RMSE indicate that the reanalyses from pNARA v1.0 provide better temporal variation but not magnitude. In terms of the PM_2.5, Fig. A2 shows no discernible differences between two experiments. This indicates that the DA system with AOD retrievals barely helps the PM simulations in GEFS-Aerosols. Note that stations over India are substantially underestimated compared with stations over the USA and Europe. This suggests that the PM_2.5 in this area involves more complicated chemical processes, a situation which is not resolved by the simple chemistry in the GOCART model. Tables A1 and A2 show the station averaged temporal statistics for biases, RMSE, and R² from pNARA v1.0, CAMSRA, and MERRA-2 with respect to AOD at 550 nm (Table A1) from AERONET and PM_2.5 from OpenAQ (Table A2). For AOD at 550 nm, we horizontally interpolate the products publicly available from each dataset to AERONET stations. The resolution of pNARA v1.0 is 0.5 by 0.5^∘, MERRA-2 is 0.5 by 0.625^∘, and CAMSRA is 0.5 by 0.5^∘. We also derived the AOD at 550 nm from the companion AERONET AOD retrievals at other wavelengths for the comparison. Table A1 shows that the CAMSRA has smaller biases among the three datasets, while pNARA v1.0 has smaller RMSE and larger R² values. (Note that the overall R² values for the whole dataset of AOD at 550 nm from pNARA v1.0, CAMSRA, and MERRA-2 are 0.768, 0.725, and 0.730, respectively.) For PM_2.5, the derivation provided by GMAO (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/FAQ/; last access: 20 January 2024) is applied to the MERRA-2 product. Similar to AOD at 550 nm, the horizontal interpolation was applied to MERRA-2 and CAMSRA, whereas it was done on a native model grid for pNARA v1.0. Table A2 shows that CAMSRA overestimates the surface PM_2.5, whereas MERRA-2 and pNARA v1.0 both underestimate these observations. Like the comparison of AOD at 550 nm, pNARA v1.0 compares favorably in RMSE and R² to the other reanalyses with respect to the temporal variation in PM_2.5. These two tables illustrate differences between the reanalyses. We note that a fairer comparison would result if all the reanalyses were available on native grids. Also, the overarching goal of this study was not to stratify various reanalyses in terms of accuracy. Instead, it was to show how significant are the differences that exist between them with respect to the mixing ratios of various aerosol species and the implications that this fact carries for modeling of the impacts of aerosols on weather and climate.

Figure A1 The horizontal distribution of temporal biases (a, d), R² correlation (b, e), and RMSE (c, f) of AOD at 500 nm against AERONET measurements of 2016.

Figure A2 The horizontal distribution of temporal biases (a, d), R² correlation (b, e), and RMSE (c, f) of PM_2.5 (µg m⁻³) against OpenAQ measurements. Evaluation period: May–December 2016.

Table A1 Comparison of station averaged temporal statistics for bias, RMSE, and R² correlation of pNARA v1.0, CAMSRA, and MERRA-2 against AERONET AOD at 550 nm over four seasons and the whole year of 2016.

Table A2 Comparison of station averaged temporal statistics for bias, RMSE, and R² correlation of pNARA v1.0, CAMSRA, and MERRA-2 against OpenAQ PM_2.5 measurements over May–December 2016.