JRAero : the Japanese Reanalysis for Aerosol v 1 . 0

A global aerosol reanalysis product named the Japanese Reanalysis for Aerosol (JRAero) was constructed by the Meteorological Research Institute (MRI) of the Japan Meteorological Agency. The reanalysis employs a global aerosol transport model developed by MRI and a 2-dimensional variational data assimilation method. It assimilates maps of aerosol optical depth (AOD) from MODIS onboard Terra and Aqua satellites every 6 hours and has a TL159 horizontal resolution (approximately 1.1° × 1.1°). This paper describes the aerosol transport model, the data assimilation system, the observation 10 data, and the set-up of the reanalysis and examines its quality with AOD observations. Comparisons with MODIS AODs that used for the assimilation showed that the reanalysis showed much better agreement than the free run (without assimilation) of the aerosol model and improved underand overestimation in the free run, thus confirming the sanity of the data assimilation system. The reanalysis had a root mean square error (RMSE) = 0.05, a correlation coefficient (R) = 0.96, a mean fractional error (MFE) = 23.7%, a mean fractional bias (MFB) = 2.8%, and an 15 index of agreement (IOA) = 0.98. The better agreement of the first guess, compared with the free run, indicates that aerosol fields obtained by the reanalysis can improve short-term forecasts. AOD fields from the reanalysis also agreed well with monthly averaged global AODs obtained by the Aerosol Robotic Network (AERONET) (RMSE = 0.08, R = 0.90, MFE = 28.1%, MFB = 0.6%, and IOA = 0.93). Site-by-site comparison showed that the reanalysis was considerably better than the free run; RMSE was <0.10 at 86.4% of the 181 AERONET sites, 20 R was >0.90 at 40.7% of the sites, and IOA was >0.90 at 43.4% of the sites. However, the reanalysis tended to have a negative bias at urban sites (in particular, megacities in industrializing countries) and a positive bias at mountain sites, possibly because of insufficient anthropogenic emissions data, the coarse model resolution, and the difference in representativeness between satellite and ground-based observations. 25


Introduction
Airborne aerosols are tiny particles that are globally distributed in Earth's atmosphere from the troposphere to the stratosphere (including over urban areas, oceans, deserts, and forests) and that affect various aspects of human society and the Earth system.In the short term and most obviously, aerosol particles can degrade visibility and damage aviation and transport (Wilkinson et al., 2012;Prata and Tupper, 2009).The World Health Organization (WHO), which has reported country estimates of air pollution exposure and its health impact (WHO, 2016), has suggested that 6.5 million deaths (11.6 % of all global deaths) may be associated with indoor or outdoor air pollution, and 92 % of the world's population lives in places where air quality levels do not meet the WHO Ambient Air quality guideline of an annual mean PM 2.5 (particulate matter with a diameter of less than 2.5 µm) concentration of less than 10 µg m −3 .Impacts on the Earth system include effects on ocean biogeochemistry and climate.For instance, 500 Mt of dust particles settles to the oceans and supplies iron, a vital element for ocean productivity (Shao et al., 2011).Aerosols also modify the radiation balance and influence the albedo and properties of clouds through both direct and indirect effects (IPCC, 2013).
In recent decades, our understanding of the aerosol life cycle and aerosol impacts has been advanced by innovation and progress in both observational and numerical modeling techniques.Among various observation techniques, Published by Copernicus Publications on behalf of the European Geosciences Union.
remote sensing using passive and active sensors has revealed both the spatial distribution and temporal evolution of aerosols at regional and global scales.The Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua satellites is an example of a passive sensor.MODIS has provided aerosol optical properties (AOPs), including aerosol optical depth (AOD) and the Ångström exponent, globally since 2000 (Remer et al., 2005;Levy et al., 2007;Zhang and Reid, 2006).More recently, Himawari 8, a geostationary meteorological satellite launched on 7 October 2014, has been providing full-disk images of AOPs every 10 min (Bessho et al., 2016;Kikuchi et al., 2017;Yumimoto et al., 2016).The Aerosol Robotic Network (AERONET; http://aeronet.gsfc.nasa.gov/) is a global ground-based network of passive instruments, including sun photometers, that has provided AOPs at several wavelengths for more than 20 years (Holben et al., 1998(Holben et al., , 2001;;Li et al., 2014).Backscatter lidar is an active-type aerosol observation sensor.The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIPSO), the first polarization lidar in orbit, has provided continuous global measurements of vertical aerosol distributions since 2006 (Winker et al., 2010;Liu et al., 2008).The European Aerosol Research Lidar Network (EARLINET; https://www.earlinet.org;Matthais et al., 2004) and the Asian Dust Network (AD-Net; http://www-lidar.nies.go.jp/AD-Net;Sugimoto et al., 2008) are regional groundbased lidar networks, and the Micro-Pulse Lidar Network (MPLNET; http://mplnet.gsfc.nasa.gov) is a global groundbased lidar network.Although satellites and ground-based networks have remarkably increased the amount of observed data that are available, observational coverage is still limited and spatially and temporally uneven.Passive satellite-borne sensors cover wide swaths across large regions, but their retrievals are usually vertically integrated AOPs that lack information on vertical distributions, and they can be obtained only under daytime and clear-sky conditions.Active sensors have the capability of measuring vertical aerosol profiles both under clouds (by ground-based lidar) and above clouds (by space-based lidar).However, their field of view is quite narrow compared to that of passive sensors so it is difficult to use lidar to cover large regions.
Numerical modeling of the aerosol life cycle has advanced considerably in the last decade.Various global aerosol transport models have been developed, and multimodel intercomparison projects have been carried out (e.g., AeroCom, Kinne et al., 2006; International Cooperative for Aerosol Prediction (ICAP), Sessions et al., 2015).Global aerosol transport models have been developed by numerous weather prediction centers, research institutes, and universities (Kinne et al., 2006;Sessions et al., 2015, and references therein).An aerosol transport model can provide three-dimensional, gridded aerosol distribution data at regular time intervals.However, multimodel intercomparison projects have shown that large differences and uncertainties due to insufficient emissions data, poor parameterization of aerosol processes (e.g., transport, chemistry, settling, and deposition), and errors in meteorological fields remain in the models.
Data assimilation, which is the integration of observation data into a numerical model, is one way of overcoming these shortcomings.At an early stage in the development of data assimilation methods, Collins et al. (2001) and Wang et al. (2006) assimilated satellite measurements by using an optimal interpolation method and a Newtonian nudging scheme.Hakami et al. (2005) and Yumimoto et al. (2007Yumimoto et al. ( , 2008) ) attempted to apply a four-dimensional variational method (a so-called advanced data assimilation method) to inverse modeling of black carbon (BC) and dust aerosols with ground-based observations and regional models.To date, measurements obtained by various observation platforms, including MODIS (Dai et al., 2014;Huneeus et al., 2012;Wang et al., 2012;Zhang et al., 2008), CALIPSO (Sekiyama et al., 2010;Zhang et al., 2011Zhang et al., , 2014)), Himawari mination of the initial and boundary conditions of not only aerosol transport models but also numerical weather forecasting models, climatological analyses of aerosols and their climate effect, satellite measurement retrievals (e.g., Yokota et al., 2009), use as truth data for observation system simulation experiments (e.g., Yumimoto, 2013;Zoogman et al., 2014), and use as input data for epidemiologic studies of PM 2.5 (e.g., Atkinson et al., 2014;Kloog et al., 2011).However, at this time, aerosol reanalysis is still at an early stage of development.
In this study, we conducted a global aerosol reanalysis named the Japanese Reanalysis for Aerosol (JRAero) for the period 2011-2015 and validated its quality with AOD measurements.The valued-added MODIS AOD observations provided by the NRL and the University of North Dakota (NRL-UND MODIS AODs) were assimilated into the MASINGAR mk-2 (Model of Aerosol Species IN the Global AtmospheRe), a global aerosol transport model developed at JMA's Meteorological Research Institute (MRI), by using a two-dimensional variational data assimilation system.This paper constitutes a comprehensive report on JRAero.The aerosol model, the data assimilation method, and the observation data used in this study are outlined in Sect.2.1-2.3, and the setup of the reanalysis is presented in Sect.2.4.Section 3 focuses on the sanity check and evaluation of the reanalysis product with MODIS AOD and independent observation data.We present our conclusions in Sect.4, and we offer remarks on future directions in the development of this reanalysis in Sect. 5.

Description of the data assimilation system
The data assimilation system consists of an aerosol model, a data assimilation module, and observation data.We describe each component in the following subsections.The setup of the reanalysis and the observation dataset used for an independent validation are also presented.

Overview of MASINGAR mk-2
The global aerosol transport model MASINGAR mk-2 (Yukimoto et al., 2012) is the first major update to the original MASINGAR (Tanaka et al., 2003), which was developed by MRI and JMA.MASINGAR mk-2 is coupled online with an atmospheric general circulation model (AGCM; MRI-AGCM3; Yukimoto et al., 2012) through a general-purpose coupler (Scup; Yoshimura and Yukimoto, 2008), and it is capable of treating the major tropospheric aerosol components, BC, organic carbon (OC), mineral dust, sea-salt, and sulfate aerosols, and their precursors.JMA started to use the original MASINGAR for operational sand and dust storm forecasting in January 2004, and it changed over to MASINGAR mk-2 in November 2014.MASINGAR mk-2 serves as a member of the ICAP multi-model ensemble (MME) (Sessions et al., 2015).The main features of the update were (1) replacement of the coupled AGCM, (2) improvement in the treatment of stratospheric sulfate aerosols, (3) renovation of emission schemes for dust and sea-salt aerosols, and (4) redesign of calculation of AOPs.Here, we provide an overview of MASINGAR mk-2 in which we focus on the updated features.

Physical processes
MASINGAR mk-2 includes advection, convective, diffusive transport, emissions, chemical reaction, and removal processes.In the model, the continuity equation of the mixing ratio of the ith aerosol component at the time step n + 1 (x i (t n+1 )) is solved by successively applying operators, as follows: where A, D, C, and T , respectively, denote operators associated with advection, eddy diffusion, convective transport, and transformations, due to emissions, depositions, and chemical reactions.
The model employs a semi-Lagrangian advection scheme (Staniforth and Côté, 1991) for advection (A).The scheme allows the use of a much longer time step compared to a Eulerian advection scheme while conserving numerical stability and accuracy.In the scheme, the upstream point is first searched with horizontal and vertical wind fields and then the mixing ratio is obtained by a three-dimensional interpolation.The interpolation includes a correction for overshooting or undershooting the advection fluxes to ensure mass conservation and a non-negative value.
Eddy diffusion (D) is defined as where K z is the vertical eddy diffusion coefficient.The coefficient for aerosols is assumed to be the same as that for water vapor and is provided by MRI-AGCM3 through the coupler.MRI-AGCM3 calculates the diffusion coefficient by the turbulence model of Nakanishi andNiino (2006, 2009), which is an improved version of the MY scheme (Mellor andYamada, 1974, 1982) called the MYNN scheme.Convective transport (C) is estimated by using updraft or downdraft mass fluxes calculated by the mass-flux-type cumulus parameterization scheme developed by Yoshimura et al. (2015) and embedded in MRI-AGCM3.Further details of convective transport are given by Yukimoto et al. (2012).
The transformation process (T ) includes gravitational settling, dry and wet depositions, emissions, and chemical reactions.Gravitational settling is estimated from the terminal velocity (V s ), which is estimated on the basis of Stokes' law as follows: www.geosci-model-dev.net/10/3225/2017/Geosci.Model Dev., 10, 3225-3253, 2017 where C c is the Cunningham slip-flow correction, ρ a is air density, ρ p and r p are the density and radius of the aerosol particles, g is gravitational acceleration, and µ is the viscosity of air.The model assumes that the particles are spherical.Ginoux (2003) investigated the effect of nonsphericity on gravitational settling and suggested that the introduction of nonsphericity to the modeling of gravitational settling does not significantly improve dust modeling.Dry deposition at the ground surface removes aerosol particles from the atmosphere.To calculate the dry deposition velocity (V d ), the model employs the resistance analog model (Seinfeld and Pandis, 2016), in which the dry deposition velocity is expressed as the inverse of the sum of resistances, where r a and r b are the aerodynamic and quasi-laminar resistance, respectively.Tanaka et al. (2003) have described derivation of the resistances in detail.MASINGAR mk-2 considers two wet deposition processes, namely, in-cloud and below-cloud scavenging.Incloud scavenging is parameterized following Giorgi and Chameides (1986), and the loss rate ic in s −1 is estimated as where F c , T c , β, and t are the fraction of the cloud with precipitation, the duration of the precipitation, the frequency of conversion of cloud water to rainwater, and the model time step, respectively.These parameters are derived from the rainwater formation rate and the cloud amount provided by the AGCM.The loss rate due to below-cloud scavenging bc is evaluated as bc = where λ bc and P are the below-cloud scavenging coefficient and precipitation intensity, respectively.To obtain λ bc , the collision efficiency between aerosol particles and raindrops is calculated by considering Brownian diffusion, interception, and inertial impaction (Slinn, 1984;Seinfeld and Pandis, 2016).Tanaka and Chiba (2005) describe the below-cloud scavenging procedure in detail.Gravitational settling and dry and wet deposition of mineral dust and sea-salt aerosols are particle size dependent.The model also includes a process for the re-emission of particles due to evaporation of raindrops.The fraction of re-emitted particles is assumed to be proportional to the amount of evaporated rainwater.

Aerosol components
The standard version of MASINGAR mk-2 represents aerosols with 10 externally mixed size distributions (i.e., mineral dust, sea salt, and five carbonaceous and three sulfate categories).The model represents mineral dust and seasalt aerosols by discrete size bins and assumes lognormal size distributions for the other aerosol components.
For mineral dust particles, the model uses a size bin method that logarithmically divides the particle size range from 0.2 to 20 µm into 10 size classes.The volumetric mean diameters of dry particles in each size bin are 0.271, 0.430, 0.681, 1.08, 1.71, 2.71, 4.30, 6.81, 10.8, and 17.1 µm.The dust particle density is assumed to be 2.65 g cm −3 following Tegen and Fung (1994).The dust emission scheme was completely replaced in the update.Now, the dust emission flux is estimated by using the wind erosion model developed by Shao et al. (1996) instead of by using the empirical approach of Gillette (1978).The dust emission flux from the ith dust size bin F i is estimated as where C d , A, and F 0i are a global tuning parameter, the erodible area fraction, and the dust emission flux, respectively, estimated by the wind erosion model.C d is set to 0.001 based on the study by Tanaka and Chiba (2005).The erodible area fraction depends on ground surface conditions and is represented by where A v , A s , A w , A l , and A t denote the land-cover factors for vegetation, snow, water, land use, and soil type, respectively.Equation ( 8) means that larger cover with vegetation, snow, or water suppresses dust emissions.Vegetation cover A v is a function of the leaf area index (LAI) and is set to 1 when LAI is larger than a threshold value (1.2), following Lunt and Valdes (2002).The area fraction of snow cover is used for A s .Thus, when the whole area of a model grid cell is covered by snow, the model estimates zero dust emissions from that grid cell.The value of LAI and the area fraction of snow cover are provided by the hydrology, atmosphere, and land (HAL) land surface model embedded in MRI-AGCM3.
The area fraction of inland water (i.e., oceans, rivers, and lakes) is used for water cover A w , and dust emission from grid cells covered by inland waters is suppressed.The landuse type is used to identify potential erodible surfaces.For instance, a grid cell covered by "evergreen broadleaf forests" is not a potential erodible surface, so its land-use factor A l = 0, whereas grid cells covered by "sand desert" have A l = 1.The land class type (LCT) database provided by the US Geological Survey (USGS; http://www.usgs.gov) is used to obtain A w and A l .Surfaces covered by Lithosol are excluded as a possible dust source (i.e., A t = 0), and soil texture data from the Food and Agriculture Organization (FAO) are used to calculate A t .
The wind erosion model estimates the friction velocity, threshold friction velocity, and saltation flux and then calculates the dust emission flux F 0i .The saltation flux of a saltat-ing particle of diameter D s is calculated as follows: Here, u * t (D s ) is the threshold friction velocity of a saltating particle of diameter D s , and u * s is the friction velocity.C s is a coefficient that depends on the saltating particle of diameter D s .ρ a and g denote the air density and the gravitational constant, respectively.The threshold friction velocity is calculated using the empirical formula of Shao and Lu (2000), and the effects of soil water are calculated in accord with Fécan et al. (1999).To consider the frozen soil effect (Kang et al., 2013), the soil water factor of Fécan et al. (1999) is enhanced by 50 % when the soil temperature drops below 0 • C. The roughness length for momentum transfer and soil water content provided by HAL are used for the estimation of u * t (D s ) and u * s .Finally, the wind erosion model calculates the dust emission flux F 0i from the saltation flux using the energy-based dust emission scheme proposed by Shao et al. (1996): Here, u * t, d is the threshold friction velocity of a dust particle, and β is an empirical function of the particle diameter of a saltating particle (D s ) and a dust particle (D d ): β = 10 −5 [1.25 ln (D s ) + 3.28] exp (−140.7Dd + 0.37) .
Carbonaceous aerosols are classified as BC or OC, and BC and OC are further divided into hydrophobic or hydrophilic states.It is assumed that 80 % of BC and 50 % of OC are emitted from both anthropogenic sources and biomass burning in a hydrophobic state and the rest in a hydrophilic state.Hydrophobic BC and hydrophobic OC both become hydrophilic through aging processes, and the conversion rate is based on an e-folding time of 1.2 days (Chin et al., 2002).Although MASINGAR mk-2 does not calculate secondary organic aerosol production explicitly, this process is represented implicitly by giving OC production from terpene using emission data.Emission amount of terpene is provided by the emission inventory.OC produced from terpene is assumed to be secondary organic aerosol and treated as hydrophilic.It is assumed that only hydrophilic species are removed by wet deposition processes.The size distributions of BC and OC are represented by lognormal distributions with a number-equivalent geometric mean radius of 0.0118 and 0.0212 µm, respectively, and standard deviations of 2.0 and 2.2, respectively, under dry conditions (Hess et al., 1998;Chin et al., 2002).The particle density of BC is assumed to be 1.25 g cm −3 , and that of OC is assumed to be 1.8 g cm −3 .
MASINGAR mk-2 has a sulfur chemistry model that treats eight major sulfur compounds, including sulfur dioxide (SO 2 ), sulfate (SO 2− 4 ), dimethyl sulfide (DMS), hydrogen sulfide (H 2 S), carbonyl sulfide (OCS), and methane sulfonic acid (MSA), and it includes nine gas-phase reactions and two aqueous-phase reactions.Tanaka et al. (2003) and Yukimoto et al. (2011) present details of the sulfur chemistry model.It is well known that OCS contributes to the formation of stratospheric sulfate aerosols (Turco et al., 1980).The chemical process including OCS was newly added to MASINGAR mk-2 to improve simulation of stratospheric sulfate aerosols.MASINGAR mk-2 can separately treat sulfate aerosols originating from anthropogenic, biogenic, and volcanic sources.The size distribution of the sulfate aerosols is also represented by a lognormal size distribution with a mean radius of 0.07 µm and a standard deviation of 2.03 under dry conditions (Hess et al., 1998;Chin et al., 2002), although different size distributions can be given for each sulfate component.The particle density of sulfate is assumed to be 1.7 g cm −3 .
MASINGAR mk-2 receives meteorological fields (e.g., wind fields, air temperature, surface pressure, clouds, precipitation, and the eddy diffusion coefficient) and surface conditions (e.g., soil water content and LAI) from the AGCM for the calculation of the advection, diffusion, deposition, and emission processes.In turn, the AGCM receives mixing ratio and deposition fluxes of the aerosol components from MASINGAR mk-2 via the coupler for the simultaneous calculation of the direct and indirect radiative effects of aerosols in its simulation.Yukimoto et al. (2012) provide a detailed description of the calculation of radiative effects.

Calculation of aerosol optical depth
The extinction coefficient at a given wavelength of aerosol component l in the kth vertical layer is calculated as follows (e.g., Tegen and Lacis, 1996;Chin et al., 2002): www.geosci-model-dev.net/10/3225/2017/Geosci.Model Dev., 10, 3225-3253, 2017 Here, Q l is the extinction efficiency factor, which is the cross-section-weighted mean extinction efficiency over a given particle size distribution of the lth aerosol component, x l,k denotes the mass concentration of the lth aerosol component in the kth vertical layer, ρ a l is the particle mass density of the lth aerosol component, and r eff l is the effective radius (cross-section-weighted mean radius over a given particle size distribution) of the lth aerosol component.Q l is calculated on the basis of Mie scattering theory for a homogeneous sphere at a given wavelength using the modeled aerosol size distributions (see Sect. 2.1.2) and the complex refractive index, obtained from the software package Optical Properties of Aerosols and Clouds (Hess et al., 1998).
Under humid conditions, all quantities of hygroscopic components in Eq. ( 13) change with RH because they take up water.The factors for hygroscopic growth with RH of the hygroscopic components (i.e., sulfate, hydrophilic OC and BC, and sea-salt aerosols) are taken from Chin et al. (2002).The complex refractive indices of the hygroscopic components are basically determined by volume-weighted averaging of the complex refractive index of water and that of each dry component.In the calculation of the extinction coefficient, the sulfate component is assumed to be ammonium sulfate and the mass concentration of the sulfate component is increased by the ammonium sulfate-to-sulfate molecular ratio to compensate for the absence of ammonium in the model.The extinction coefficient of organic aerosols is estimated by replacing the OC mass concentration with the organic matter (OM) mass concentration using an OM-to-OC factor of 1.4 (White and Roberts, 1977;Japar et al., 1984;Russell, 2003).MASINGAR mk-2 calculates the aerosol extinction coefficients at wavelengths of 550 and 870 nm.The total AOD is derived by integration of α l,k in all aerosol components and all model vertical layers as follows: Here, L, K, and z k are the number of aerosol components, the number of model vertical layers, and the box height between the upper and lower boundaries of the kth vertical layer, respectively.We used the modeled AOD value at the typical visible wavelength of 550 nm in this study.

Assimilation method
We previously developed an aerosol data assimilation system (MASINGAR/LETKF) based on MASINGAR mk-2 and an ensemble-based assimilation technique called the local ensemble transform Kalman filter (LETKF; Hunt et al., 2007;Miyoshi and Yamane, 2007).In fact, Yumimoto et al. (2016) have reported successful results with MASINGAR/LETKF in the assimilation of products from MODIS and Himawari 8.However, the computational cost of MASINGAR/LETKF is quite high due to the necessity of ensemble simulation, and it is still unrealistic for development of the long-term reanalysis product.Therefore, in this study, for the initial development of a reanalysis product, we developed an aerosol data assimilation system based on a sequential variational method.
The assimilation system (MASINGAR/2D-Var) uses a twodimensional variational method (2D-Var) for the core of the assimilation system.NAVDAS-AOT (Zhang et al., 2008(Zhang et al., , 2014) ) and the NAAPS reanalysis also employ 2D-Var (Lynch et al., 2016).Experience and knowledge obtained during the development of MASINGAR/LETKF, and by conducting experiments with that system, were utilized in the development of MASINGAR/2D-Var.
The cost function in a three-dimensional variational method (3D-Var) is generally defined as where x denotes the vector of modeled aerosol mass mixing ratios.The suffix f represents the forecast (a priori), and y o denotes a vector that contains observations used for the assimilation.H is an observation operator that transforms model variables to observation space.In this study, we used the NRL-UND MODIS AODs (τ o ) as the observational constraint (y o ).Therefore, the observation operator includes the conversion of the aerosol mass mixing ratio into AODs (described in Sect.2.1.3)and the interpolation of model space into observation space.The first term on the right-hand side of Eq. ( 15) represents the departure of the analysis (a posteriori) aerosol mass mixing ratio from the forecast ratio weighted by the background error covariance matrix for the aerosol mass mixing ratio (P).The second term on the righthand side represents the difference between the modeled and observed AODs weighted by the observation error covariance matrix (R).To search for the optimal solution, which minimizes the cost function, the gradient vector of Eq. ( 15) is calculated as follows: Here, H T is an adjoint of H .For the assimilation of two-dimensional observations such as the NRL-UND MODIS AODs, we reduced the computational cost by degrading the assimilation system from three-dimensional to two-dimensional and analyzed the modeled AOD (two-dimensional variable) instead of the modeled aerosol mixing ratio (three-dimensional variable).The cost function (Eq.15) is redefined as where τ denotes the modeled AOD vector, and P τ is the background error covariance matrix for AOD.H I is the interpolation into observation space and a linear operator.The gradient vector of Eq. ( 16) is derived as follows: Because the observation operator H I is a linear operator, the cost function (Eq.17) becomes a quadratic scalar function.
Therefore, the minimum of the cost function is obtained for τ = τ a , which makes ∇J τ = 0.The analysis increment of AOD is derived from Eq. ( 18) as Suffix a represents the analysis (a posteriori), τ o − H I τ f is the innovation, and K is the Kalman gain, defined as follows: Because the MODIS AOD is the product of a columnintegrated aerosol optical property and does not provide any information about the vertical profile or the aerosol components, we allocate the analysis increment of AOD to the aerosol mass mixing ratio while keeping the shape of the vertical profile and the rate of each aerosol component in the forecast (a priori) aerosol fields.First, the analysis increment of AOD is derived for each aerosol component.The analysis increment of the AOD for the lth aerosol component is calculated as where τ f l is the forecast AOD of the lth aerosol component.Then, we distribute the analysis increment of AOD of each aerosol component to the aerosol mixing ratios in each vertical layer.The analysis increment of the mixing ratio of the lth aerosol component and the kth vertical layer is derived as where x f l,k and α f l,k represent the forecast mixing ratio and the extinction coefficient of the lth aerosol component and the kth vertical layer, respectively.In this distribution, we assume that the hygroscopic growth rate is unchanged.Finally, we obtain the analysis mixing ratio as follows: MASINGAR/2D-Var is developed with expansibility that allows the assimilation of three-dimensional observations in the future update.
We introduced a localization technique used in LETKF to the system that divides the model space into local regions using a prescribed localization scale.The localization technique solves the analysis increment of AOD (δτ a ) at each model grid with observations included in the local region independently, reduces spurious error covariance with distance, and enables parallel processing to be used to reduce computational cost.
We define the background error covariance between model grids m and n (i.e., the (m, n) element of P τ ) as where σ m τ and σ m τ are the background error standard deviations at model grids m and n, respectively.C m,n is a smooth weighting function whose value becomes smaller when the distance between model grids m and n becomes larger.The function prevents a spurious large covariance between distant model grids.We use the second-order auto-regressive (SOAR) approximation (Daley and Barker, 2001) following Zhang et al. (2008) to calculate the value of C m,n : where R m,n denotes the distance (km) between model grids m and n, and L is the horizontal error correlation length.Zhang et al. (2008) calculated the spatial correlation between satellite observations and model forecasts (1 • × 1 • horizontal resolution) as a function of distance and found that the SOAR with L set to 200 km can fit the correlation and when the distance was more than 1000 km, the spatial correlation decreased to less than 0.05.On the basis of their results, we set the localization scale and the horizontal error correlation length to 1000 and 200 km, respectively.
The background error standard deviation was derived from the MASINGAR mk-2 simulation without data assimilation (free run).We collected AOD values within ±15 days of the targeted hour, and then calculated their mean value (τ FR ) and standard deviation ( σFR ).For instance, for 12:00 UTC on 25 July 2014, AOD values at 12:00 UTC on 15 July-9 August were gathered (i.e., 31 AOD values) and used for the calculation of the mean and standard deviation.We calculated the background error standard deviation as follows: Here, τ f is the forecast AOD.The fraction on the righthand side of Eq. ( 26) indicates that the standard deviation is normalized by the mean value.(a) Whole period  fields.To introduce the horizontal structure of AOD field at the assimilation time into the background error covariance, the normalized standard deviation is multiplied by the forecast AOD.

Assimilation dataset: the NRL-UND MODIS AODs
We employed the MODIS Level 3 AOD product provided by the NRL and UND (https://earthdata.nasa.gov/earthobservation-data/near-real-time/nrt-value-added-modisaerosol-optical-depth-product-available)as the assimilation dataset.The NRL-UND MODIS AOD product was produced for the purpose of aerosol data assimilation and is based on the NASA operational MODIS Level 2 Collection 5 AOD dataset (Remer et al., 2005;Levy et al., 2007).These data (i.e., Dark Target AODs) have been subjected to extensive quality assurance (QA) and quality check (QC) procedures (Zhang and Reid 2006;Shi et al., 2011;Hyer et al., 2011).
The QA and QC procedures included (1) a stringent filter to reduce outliers, eliminate cloud contamination, and exclude bad conditions when aerosol detection is likely to be inaccurate, and (2) empirical corrections to reduce systematic biases over land and ocean.Then, the data are aggregated onto a 1 • × 1 • grid to reduce random spatial variation in AOD values with additional buddy checks.The product is provided every 6 h (i.e., 00:00, 06:00, 12:00, 18:00 UTC).The QA and QC procedures and the validation of the processed data have been reported by Zhang and Reid (2006), Shi et al. (2011), andHyer et al. (2011).Although the cutoff latitudes for the product are 40 • S over water in the Southern Hemisphere and 80 • N in the Northern Hemisphere, in this study, the AOD product observed from 60 • S to 60 • N was assimilated.We assumed that the observation error covariance matrix (R) was diagonal and assigned uncertainty of AOD provided by the NRL-UND MODIS AOD product to the diagonal component of the observation error covariance matrix.The uncertainty includes empirical estimation of observation error and representativeness error based on variability of the L2 dataset (Zhang and Reid, 2006).
The horizontal distribution of the total number of NRL-UND MODIS AOD data assimilated in the reanalysis run (RA) of MASINGAR mk-2 during the entire reanalysis period and seasonally is shown in Fig. 1.Little or no observation data are available from regions covered by bright deserts (e.g., the Sahara, the Kalahari Desert of southern Africa, the Arabian Peninsula, the Middle East, central Asia, inland deserts of China, and the inland desert of Australia).Cloud coverage also affects the availability of observation data.For instance, the Intertropical Convergence Zone, which is frequently covered by clouds, has relatively fewer observation data.During winter in the Northern Hemisphere (i.e., December-February), little or no observation data are available from higher latitudes because of snow or ice cover.There are no AOD observations over the ocean south of 40 • S because of the cutoff used by the QA/QC procedure.This cutoff aims to filter out cloud contamination in the Southern Ocean (Toth et al., 2013).The time series of the number of the NRL-UND MODIS AOD data assimilated in RA shows that there are more observation data during the boreal summer than during the boreal winter (Fig. 2a).No observation data were available from 06:00 UTC on 1 October 2013 to 09:00 UTC on 9 October 2013.On average, about 1400 observations were assimilated at each analysis time.We assimilated about 2 000 000-2 900 000 data each year.During the whole reanalysis period, about 13 000 000 observation data points were available.

Evaluation dataset: the AERONET AOD
We employed the AERONET product (version 2, level 2.0 (quality-assured)) for the independent evaluation of the reanalysis.Various studies of AOP retrieval and aerosol simulation and assimilation have used the AERONET product for validation (Kinne et al., 2006;Levy et al., 2013;Sessions et al., 2015;Lynch et al., 2016;Rubin et al., 2016).Because the AERONET product does not include AODs measured at 550 nm, we used the Ångström law to derive AODs at 550 nm from the AODs and Ångström exponents measured at multiple wavelengths (340-870 nm).Then, the AERONET AODs were averaged into 1 h bins and paired with the model results.We used observation data from 277 AERONET sites situated between 60 • S and 60 • N from which observations covering more than 5 % of the reanalysis period were available (Fig. 3).Roughly speaking, we obtained observation data from about 190 sites each month, but the number of available data varied seasonally and reached a maximum in boreal summer (Fig. 4).In the most recent year (i.e., 2015), the numbers of both available stations and data were smaller, because it takes time for level 2.0 (quality-assured) data to be established; the level 2.0 data become available only after final calibration has been applied and manual data inspection has been completed.To compare the modeled AODs with the AERONET AODs, the gridded modeled values were linearly interpolated to the AERONET site locations.

Experiment setup
The configuration of the reanalysis is summarized in Table 1.nate system.At the three lowest levels, the vertical resolution (layer thickness) was about 100, 200, and 300 m.The time step of the aerosol model was 900 s, and we used hourly model output data for the evaluation.The AGCM has the same spatial resolution and time step as MASINGAR mk-2.The exchange of meteorological and aerosol variables between the AGCM and MASINGAR mk-2 through the coupler is performed at every model time step (i.e., every 900 s).
At the assimilation step, the exchange follows the assimilation.This means that the AGCM receives the analyzed aerosol fields from MASINGAR mk-2.In this study, the operational global analysis provided by JMA (GANAL/JMA; JMA, 2002) at 6 h intervals was used for nudging the AGCM.
The nudging scheme is applied to the horizontal wind components and air temperature.The nudging term is applied at each time step of the integration by temporarily interpolating the variables by linear interpolation.We used anthropogenic and biomass burning emissions of sulfur dioxide, BC, and OC from the MACCity (MACC/CityZEN EU projects) emission inventory (http://accent.aero.jussieu.fr/MACC_metadata.php) and the Global Fire Assimilation System (GFAS) dataset (http://www.gmes-atmosphere.eu/about/project_structure/input_data/d_fire).This MASINGAR mk-2 setup, excepting the model resolution and the assimilation setting, is identical to the operational forecasting setup of the ICAP MME.
Figure 5 shows a schematic diagram of the reanalysis procedure.For this study, we performed two runs: a model simulation without data assimilation (free run; hereafter FR) and the reanalysis run (RA), in which NRL-UND MODIS AODs were assimilated every 6 h.The reanalysis period was from 1 January 2011 to 31 December 2015 (5 years) with a spinup period of 3 months (October-December 2010).We also compared the 6 h forecast AODs from the analyzed state (first guess; hereafter FG) with the NRL-UND MODIS AODs.
3 Evaluation results

Chi-square test
We used the chi-square test (χ 2 ), which evaluates the balance between the innovation and the background and observation error covariances, to evaluate the long-term stability of the assimilation performance (Ménard et al., 2000;Miyazaki et al., 2012Miyazaki et al., , 2015)).In this study, the chi-square value was defined as follows:  No observation data were available from gray areas (also see Fig. 1a).
Here, m is the number of observation data.When the background and observation error covariances properly balance the innovation, χ 2 has the ideal value of 1. χ 2 > 1 indicates overconfidence (i.e., underestimation) of the model or observation errors.
The time evolution of χ 2 during RA (Fig. 2b) shows that it decreased in the first 2 months of the reanalysis period and then remained approximately constant with an average value of 0.30 and a standard deviation of 0.077.This result confirms that the assimilation performance was stable.The relatively large χ 2 value in October 2013 was caused by the lack of assimilation data at that time.In almost all cases, χ 2 was less than the ideal value of 1, which implies persistent overestimation of the background or observation errors with respect to the innovation.We performed an additional 2-year assimilation experiment (2011-2012) in which the background error covariances were uniformly decreased by 60 %, and found that the average χ 2 value was 0.56 with a standard deviation of 0.16.Although the average was increased in this experiment, almost all of the χ 2 values were still lower than the ideal value, which implies persistent overestimation of the observation errors.Both RA and FG from the additional experiment obtained worse scores in the validations with MODIS and AERONET AODs than the standard experiment.For the χ 2 value, the additional experiment shows much larger variation (standard deviation).These results imply that although there were the persistent overesti-mates of background and observation errors, they were wellbalanced and stable in the standard experiment.

Sanity check by MODIS AODs
To validate the data assimilation system, we compared the AOD fields from FR, FG, and RA with the NRL-UND MODIS AODs.Better agreement of RA AODs, compared to FR AODs, confirms the accuracy of the data assimilation system, because NRL-UND MODIS AODs were used as an observational constraint.The FG performance is an indication of whether the analyzed aerosol fields can improve the short-term forecasting until the next analysis (until the next MODIS AOD is available).
The RA AOD spatial distribution showed quite good agreement with the NRL-UND MODIS AOD distribution (Fig. 6a and b).The distribution of the 5-year averaged difference (RA AOD minus FR AOD) (Fig. 6c) shows that, in general, assimilation increased AODs over the central Pacific Ocean, implying that in FR, MASINGAR mk-2 underestimated sea-salt aerosols.Large negative and positive differences over Canada, Siberia, and Indonesia resulted from improved simulation of carbonaceous aerosols from forest fires.Over the Sahel (south of the Sahara) and its downwind regions (e.g., the Atlantic Ocean), mineral dust particles were increased by assimilation.In contrast, over other dust-dominant regions (central China, Australia, the Persian Gulf, and Argentina), AODs (mainly mineral dust particles) were decreased.The negative difference around the Mediterranean Sea indicating that the model overestimate dust particles transported from the Sahara.Over industrializing areas, such as India and the eastern coast of China, assim-ilation increased anthropogenic pollutants (i.e., sulfate and carbonaceous aerosols), filling the gap between observed and FR AODs.The large positive differences around central Africa and the Gulf of Guinea reflect increased carbonaceous aerosols due to forest fires.Figure 6d exhibits the distribution of the increment (RA AOD minus FG AOD) that is derived from the 5-year average of the modifications by the assimilations (i.e., δτ a in Eq. 19).The increment shows lower amplitudes and different distributions in several regions (particularly in the downwind regions of aerosol sources) compared to the difference (Fig. 6c).It is because that the difference is affected by transport of the modifications after the assimilations and the increment (FG AOD) takes into account accumulation of previous assimilations.The effect by the accumulation also appears as much better statistics of FG AOD (e.g., lower root mean square error (RMSE) and mean fractional bias (MFB)) compare to FR AOD (Fig. 7).
Temporal evolution of the RMSE, linear Pearson correlation coefficient (R), mean fractional error (MFE), MFB, and index of agreement (IOA) for FR, FG, and RA AODs are shown in Fig. 7 (formulations of these statistical measures are given in Appendix A).For FR, these statistical measures (blue lines and dots) showed seasonal variations reflecting the seasonal cycle of the observation coverage (shown in Fig. 1) and the aerosol distribution (e.g., springtime Asian dust storms).Occasional large RMSEs and low R values for FR AODs (e.g., July 2011 and March 2015) were caused mainly by large-scale dust storms or forest fires.Values of the RA AOD statistical measures (red lines and dots) are much better than the FR values.In most cases, the RMSEs are less than 0.06 (the 80th percentile of RMSE for RA is 0.057), and R and IOA values are larger than 0.9 (20th percentile of R and IOA for RA is 0.91 and 0.94, respectively).The 80th percentile of MFE for RA AOD is 26.6 %.The MFB time series showed that assimilation considerably reduced the bias found in FR.During 96.4 % of the reanalysis period, MFB for RA AODs is within ±10 %.In most cases, RA AODs meet the model performance goals for particulate matter and light extinction (MFE ≤ +50 % and MFB ≤ ±30 %) suggested by Boylan and Russell (2006).Moreover, all FG AOD (green lines and dots) statistical measures are better than the FR AOD statistics.This result means that aerosol fields obtained by the reanalysis improved short-term forecasting until the next analysis (until the next MODIS AOD is available).
Scatter plots of FR, FG, and RA AODs versus NRL-UND MODIS AODs for the whole reanalysis period (Fig. 8) and by season (Fig. 9) show no seasonal dependency of RA.However, MFB values are relatively larger in boreal spring and summer, perhaps reflecting larger AOD values over the Northern Hemisphere during those seasons (Table 2).RA AODs are much more aligned with the 1 : 1 line in the scatter plots compared to both FR and FG AODs, and the RA statistical measures are much better, both seasonally and for the whole analysis period (RMSE ≤ 0.05, R ≥ 0.95, MFE ≤ 25.1 %, MFB = 1.3-4.1 %, and IOA ≥ 0.97; Table 2).The scatter diagrams show that RA AODs tended to be underestimated with respect to NRL-UND MODIS AODs.In fact, the mean bias (averaged value of RA AODs minus NRL-UND MODIS AODs) is slightly negative (−0.007).However, MFB is slightly positive (2.8 %).This discrepancy implies that although RA underestimated relatively large AODs (> 0.5), smaller AODs (< 0.5) were overestimated.The FG AODs also show better agreement than FR AODs with the NRL-UND MODIS AODs throughout the reanalysis period (Table 2) and meet the model performance goal.The frequency distribution of FR AOD deviations (observed AODs minus modeled AODs) (Fig. 10) is symmetric and mound-shaped, similar to a Gaussian distribution, and supports the assumption of the data assimilation that Geosci.Model Dev., 10, 3225-3253, 2017 www.geosci-model-dev.net/10/3225/2017/  the background errors follow a Gaussian distribution.However, the requirement of an unbiased background error is not strictly met (mean bias between FR AODs and NRL-UND MODIS AODs is −0.032).FG and RA AOD deviations shows squeezed distributions, and their peaks are closer to 0 than the FR peak; 92.0 % (79.6 %) of RA (FG) AOD deviations are within ±0.05.

Evaluation with 1 h binned data
We compared the modeled AODs with the 1 h binned AERONET AODs during the whole reanalysis period (Fig. 11) and during each season (Fig. 12) in scatter plots.In general, the FR AOD distribution is squashed vertically, an indication that FR AODs generally underestimated AERONET AODs.This underestimation is also reflected in negative MFB values (Table 3).In boreal spring (Fig. 12c) and summer (Fig. 12e), FR AODs in regions where AERONET AODs were between 0.0 and 1.0 showed positive biases, which were caused by the overestimation of carbonaceous aerosols from forest fires in Canada and mineral dust aerosols over central China, Australia, and the Persian Gulf mentioned in Sect.3.2.RA AODs (Figs. 11b, 12b, d,  f, and g) are more narrowly distributed along the 1 : 1 line, an indication that assimilation improved under-and overestimates in the FR AODs.The statistical measures of RA for the whole reanalysis period (Table 3) meet the model performance goals.The MFE and MFB values obtained in boreal summer compared better to those obtained in other seasons.This improvement might reflect the relatively larger number of NRL-UND MODIS AODs available in boreal summer (see Table 2).

Evaluation with monthly averaged data
We also evaluated monthly averaged AODs, because monthly AOD data are often used in climate studies (Lynch et al., 2016).Monthly averaged AODs were derived from paired 1 h binned AERONET and modeled AODs.In general, the monthly averaged RA AODs (Table 4) showed better agreement than the 1 h binned AODs (Table 3) with the AERONET AOD data; in particular, RMSE and MFB of RA were about 43 and 86 % lower, respectively.As in the pairwise comparison with the 1 h binned AODs, FR AODs underestimated the monthly averaged AERONET AODs (Fig. 14a and negative MFB value in Table 4), whereas RA AODs reduced or eliminated the underestimation.As a result, the RA MFB value, 0.6 %, is quite good (Table 4).Moreover, in the monthly data, 74.0 % (89.4 %) of RA AOD deviations from AERONET AODs were within ±0.05 (±0.10).
Time series of the statistical measures calculated by using the monthly averaged AERONET AODs (Fig. 15) show that, in general, RA AOD performance is much better than FR AOD performance throughout the reanalysis period, except for MFB in December 2015.The only statistical measure that shows seasonality is MFE.RMSE of RA is almost always (58 months during 5 years) lower than 0.10; 53.3 % (32 months) of R values and 73.3 % (44 months) of IOA values are greater than 0.90; 80.0 % (48 months) of MFE and MFB values are less than 34 % and within ±10 %, respectively; and all months except for the initial month of the reanalysis (i.e., January 2011) meet the model performance goals.
Statistical measures for RA at each AERONET site are shown in Fig. 16     Fig.17 as a reference.In these figures, the 181 AERONET sites with more than 36 monthly averaged data are plotted.In general, RA shows much better agreement globally than FR with the AERONET data.In RA, RMSE < 0.01, R > 0.90, and IOA > 0.90 at 86.4 % (154 sites), 40.7 % (74 sites), and 43.4 % (79 sites) of the 181 sites, respectively.MFE is less than 50 % at 181 sites (91.2 %), and MFB values at 81 sites (44.5 %) are within ±15 %.Consequently, at 44.5 % (81 sites) of the 181 AERONET sites, RA meets the model performance goal.RMSEs are relatively large at sites in central Africa, India, southeast Asia, and eastern coastal China, but the other statistical measures (e.g., R and MFE) in these regions are not particularly worse than those in the other regions.Relatively higher AOD values (see Fig. 6) are responsible for the larger RMSE values.At Beijing (39.977 • N, 116.381 • E; Fig. 18b), RMSE (0.28) and the negative MFB value (−49.3 %) are large, and IOA is relatively low (0.65), but the correlation coefficient is relatively high (R = 0.80).We also compared time series of monthly averaged AERONET and modeled AODs with the time series of monthly averaged NRL-UND MODIS AODs there in Fig. 18.It should be noted that the monthly averaged MODIS AODs are derived by using all available data, and neither the AERONET nor the modeled AODs are paired.At Beijing, although RA successfully captured the temporal variation observed by AERONET, its AOD values were almost always lower than the observed values.A similar negative MFB value when AOD levels are high has been observed for other megacity AERONET sites in industrializing countries, including Xi-anHe (39.754 • N, 116.962 • E; about 56 km from Beijing), Kanpur, India (26.513 • N, 80.232 • E; Fig. 18h), and Mex-ico_City, Mexico (19.334 • N, 99.182 • W).The statistics for these sites are summarized in Table 5. Averaged RA AODs are lower than averaged AERONET AODs, and MFB values are negative (−22.0 to −53.7 %), but correlation coefficients are relatively high (R = 0.78-0.79).Multi-model intercomparison studies (Kinne et al., 2006;Sessions et al., 2015) have pointed out that aerosol models having negative biases for high-AOD events is a common problem.Insufficient anthropogenic emissions data and model resolution for megacities are plausible reasons for the negative biases.Zheng et al. (2015) evaluated the impact of heterogeneous chemistry with regional chemical transport model in eastern and central China (urban and industrialized area of China), and suggested the significant role of heterogeneous chemistry in regional haze formation.While the current version of MASINGAR mk-2 includes the nine gas-phase and two aqueous-phase reactions of the sulfate chemistry, the implementation of the heterogeneous chemistry reactions is under development.The absence of the heterogeneous chemical productions may partly explain the negative bias.Assimilation partly reduced the negative bias (see Fig. 18b and h), although the amount of improvement was limited.The probability of a successful retrieval can be reduced during high-AOD events (Lynch et al., 2016); thus, fewer available satellite observations over megacities during high-AOD events may also account for the negative biases in RA AODs.Rubin et al. (2017) applied an ensemble-based assimilation method to NAAPS and found that flow-dependent error covariances estimated by ensemble simulations utilized the AERONET AOD efficiently and brought better analyses at Beijing and Kanpur.Sophistication of the background error covariance and assimilation of additional observations have the potential to improve the analyses in the megacities.
At Lulin (23.469One characteristic shared by these sites is a relatively low AOD level (0.025-0.056), and another shared feature is that all four sites are situated at high elevation in mountainous areas (Table 5).The coarse model resolution might obscure the effect of local terrain and result in this positive bias.In addition, the NRL-UND MODIS AODs also overestimated the AERONET AODs because of the difference in the rep-  resentativeness of the observations.This fact can explain the limited improvement by the assimilation.
In addition to the sites discussed above, we show time series of monthly averaged AERONET and modeled AODs of other sites (Fig. 18), which we selected as representative of urban areas (Goddard Space Flight Center -GSFC, Moldova, and Osaka) or because they show the influence of African dust (IER_Cinzana, Capo_Verde, and Ragged_Point), Arabian dust (Mezaira), African smoke (Ascension_Island), southeast Asian smoke (Chiang_Mai_Met_Sta), or South American smoke (Rio_Branco).Cart_Site was selected as representative of a clean area.Nine of these sites were also selected for validation of the ICAP-MME (Sessions et al., 2015).At the urban sites (Fig. 18f, k, and i), assimilation improved the negative bias in FR AODs, and RA AODs showed much better agreement with the observations than FR AODs, in particular at GSFC (Fig. 18f), although a slight negative MFB (−1.9 %) remained in RA AODs.At IER_Cinzana (Fig. 18g), situated in an African dust source region, RA AODs success-fully captured the observed temporal variation (R = 0.93, IOA = 0.89), but a negative MFB (−32.7 %) remained after the assimilation.The availability of fewer NRL-UND MODIS AOD observations over bright surfaces (also see Fig. 1) can explain the limited improvement achieved by assimilation.At Capo_Verde (Fig. 18c) and Ragged_Point (Fig. 18m), which are downwind of an African dust source, the RA results showed nearly perfect agreement with observations, though the sites are on opposite sides of the Atlantic Ocean.At Mezaira (Fig. 18j), which is influenced by Arabian dust, AERONET captured the AOD peaks in boreal summer.FR considerably overestimated the AOD values of the peaks, and assimilation improved the overestimation, even though no MODIS data are available over this AERONET site.MODIS observations near the site (i.e., over the Persian Gulf; see Fig. 1) contributed to this improvement, and the overall performance of RA was good (Fig. 18j).At Ascension_Island (Fig. 18a), which is influenced by African smoke, assimilation improved the persistent negative bias.In RA AODs, RMSE = 0.03 and R = 0.92, and the negative MFB (−3.0 %) was much smaller than that in FR AODs.At Chiang_Mai_Met_Sta, which is influenced by southeast Asian smoke (Fig. 18e), RA AODs reproduced the seasonal variation (peaks in spring) seen in the AERONET AODs (R = 0.96), but they underestimated peak AOD values (MFB = −12.9%).The NRL-UND MODIS peak AOD values were also lower than the AERONET peak AODs, because larger spatial variability included in the dense aerosol events filtered the high AOD values out through the QA process.This fact can explain the negative bias of the peak AOD values in the RA AODs.At Rio_Branco (Fig. 18n), which is influenced by South American smoke, RA AODs successfully captured both the temporal variation and the winter peaks (R = 0.93), but peak AOD values were slightly underestimated (MFB = −1.5 %).At the clean site (Cart_Site; Fig. 18d), assimilation complemented the underestimates in FR AODs, and the RA performance was good (RMSE = 0.02, R = 0.94, MFE = 15.3 %, MFB = −2.5 %, and IOA = 0.96).

Conclusions
A global aerosol data assimilation system was developed based on a global aerosol transport model, MASINGAR mk-2, and a two-dimensional variational data assimilation method.The assimilation system was used to produce an aerosol reanalysis product named the Japanese Reanalysis for Aerosol (JRAero) for the period 2011-2015 through the assimilation of NRL-UND MODIS AOD data with a horizontal resolution of TL159 (about 1.1 • × 1.1 • ).In this paper, we have outlined the data assimilation system and presented the general specifications of JRAero.We have also presented results of our validation of the reanalysis AODs with AERONET AODs.
Chi-square test results confirmed that the stability of the assimilation performance was good throughout the reanalysis period.In almost all cases, the chi-square value was lower than the ideal value of 1; this result implies persistent overestimation of the background and observation errors with respect to the innovation.
Comparison of the reanalysis results with the NRL-UND MODIS AODs showed that assimilation improved both negative and positive biases in the free run (FR) of MASIN-GAR mk-2: negative biases over oceans caused by underestimation of sea-salt aerosols, gaps in carbonaceous aerosols from forest fires in Canada, Siberia, Indonesia, and central Africa, discrepancies in dust source and downwind regions (e.g., the Sahel, Atlantic Ocean, Mediterranean Sea, central China, Australia, and the Persian Gulf), and underestimates in industrializing areas (in particular, India and eastern China) were compensated by the assimilation.The reanalysis AODs showed quite good agreement with the NRL-UND MODIS AODs (RMSE = 0.05, R = 0.96, MFE = 23.7 %, MFB = 2.8 %, and IOA = 0.98), confirming the accuracy of the assimilation system.FG AODs also showed better agreement with the NRL-UND MODIS AODs than the FR AODs (although worse compared to the reanalysis AODs).This result indicates that aerosol fields provided by the reanalysis are capable of substantially improving short-term forecasting till the next analysis (the next MODIS AOD is available).
The statistics of the reanalysis AODs in a comparison with monthly averaged AERONET AODs (RMSE = 0.08, R = 0.90, MFE = 28.1 %, MFB = 0.6 %, and IOA = 0.93) were even better than those in the comparison with 1 h binned AERONET AODs.In site-by-site comparisons, the reanalysis performance was also better than the FR performance.
At 86.4 % of the 181 AERONET sites, the RMSE of RA was < 0.10; at 40.7 % of sites, R was > 0.90; and at 43.4 % of sites, IOA was > 0.90.However, the reanalysis AODs tended to underestimate the observed AODs at urban sites (in particular, megacities in industrializing countries), possibly because anthropogenic emissions data and model resolution were insufficient.At high-elevation mountain sites (2500-4200 m), persistent positive biases were found in the reanalysis and improvement by assimilation was limited.The coarse model resolution, which likely obscures the effect of local terrain, and the difference in representativeness between satellite and ground-based observations can explain the overestimation and limited improvement.

Future directions
To enhance the current version's quality and address its problems, we propose the following for the next version: 1. NRL-UND MODIS AODs (or Dark Target AODs) are unavailable over deserts.These regions without observational constraints are responsible for the limited improvement of RA obtained near dust source regions.Inclusion of AOD data retrieved by the Deep Blue algorithm, which is able to complement Dark Target AODs by retrieving AODs over bright arid land surfaces (Hsu et al., 2006;Sayer et al., 2013), could improve the reanalysis quality in regions around deserts.In the current version, we assimilated two-dimensional maps of AODs (vertically integrated aerosol optical property) and assumed that the shape of vertical profile before assimilation was unchanged after assimilation.However, vertical profiles of aerosols likely affect their transport, deposition, and climate effects.A space-borne lidar, CALIPSO (Winker et al., 2010), has been providing continuous measurements of aerosol vertical distributions over the globe since 2006.The use of lidar data would allow the vertical profiles of the reanalysis to be adjusted.Furthermore, the inclusion of information on the size distribution (e.g., Ångström exponent) would raise the quality of the reanalysis.
2. The limited temporal coverage of MODIS AOD (once a day) might cause temporal jumps or discontinuities in the reanalysis.Lynch et al. (2016) suggested that the importance of model tuning in particular when there are areas not covered by assimilated observations.To obtain a better coverage of assimilation data is important for the future development.Smoother techniques (e.g., four-dimensional variational method; Yumimoto andTakemura, 2013, 2015) and ensemble Kalman smoother (Schutgens et al., 2012) are also useful.
remained in the reanalysis.A plausible reason is the coarse model resolution, which is insufficient to resolve high-AOD events around megacities and local terrain effects in mountainous areas.Therefore, we plan to rerun the reanalysis with a finer resolution and check the performance of the model.Re-examination of the background error by an ensemble-based method (Yumimoto et al., 2016) and assimilation of additional observations (e.g., the AERONET AOD) have the potential to improve the analyses at the megacity sites (Rubin et al., 2017).
4. In the next version, we will employ the JRA-55 reanalysis dataset (Kobayashi et al., 2015;Harada et al., 2016), the most recent meteorological reanalysis product developed by JMA, as nudging data in the AGCM.The use of these data is expected to enable more accurate simulation of aerosol transport and deposition and to allow longer integration of the reanalysis.
5. Updating to more recent anthropogenic emissions inventory data will improve the negative biases, especially for industrializing countries where anthropogenic emissions are rapidly increasing.
6.In addition to the updates of the model resolution and input data, physical and chemical processes in the model need to be sophisticated.Newly developed dust emission (e.g., Kang et al., 2014), optical calculation, microphysics parameterization, and wet deposition schemes (e.g., Oshima et al., 2009;Oshima and Koike, 2013) will be applied in the next version.Aerosol microphysical and optical parameters (e.g., size distribution, hygroscopicity, and refractive index) will also be updated by taking advantage of recent measurements.
7. The chi-square test results imply model and observation errors are relatively large with respect to the innovation.Re-examination of these errors is required.Insights from studies using the ensemble-based assimilation system (Yumimoto et al., 2016) should make it possible to provide the background error covariance matrix with a more appropriate amplitude and structure.In the current version, the background error was in proportion to the forecast AOD (Eq.26) and became small where the model did not predict aerosols.Therefore, the analysis could not reproduce aerosol events that satellites observed but the model failed to predict (e.g., dust storms and biomass burning) due to the small background error.The ensemble-based estimate of the background error considering uncertainty in emissions will bring better analysis for this situation.
8. In the present paper, the reanalysis was validated with only the AOD observations.Further validation is also needed.Comparison with vertical AOD distributions obtained from space-and ground-based lidars would contribute to improvement of the vertical profiles of the reanalysis and become a good prearrangement of assimilation with the lidar data.Validation using mass concentrations observed in situ and by aircraft campaigns should effectively improve the model's emission, transport, diffusion, and deposition processes.Furthermore, intercomparison with other reanalysis products (e.g., CAMSiRA, MERRA-2, and NRL NAAPS reanalyses) will provide insight into various aspects of not only the data assimilation system but also the aerosol model and help us increase their sophistication.

Figure 2 .
Figure 2. (a) Time series of the number of 6-hourly NRL-UND MODIS AOD data.(b) Time series of the 5-day moving average of the chi-square value (χ 2 ; thick black line) and its standard deviation (gray shading).

Figure 3 .
Figure 3. AERONET sites used in this study.Red circles denote the AERONET sites used in Fig. 18.

Figure 4 .
Figure 4. Time series of the number of monthly AERONET AOD data (gray shading) and the number of AERONET sites (thick black line) used in the evaluation.

Figure 5 .
Figure 5. Schematic diagram of the reanalysis procedure.A map showing NRL-UND MODIS AODs is also shown at each assimilation time.

Figure 10 .
Figure 10.Frequency distributions of deviations (observed AODs minus modeled AODs) from the NRL-UND MODIS AODs.The percentages of deviations between −0.05 and +0.05 are shown at the top, and the percentages less than −0.5 (bottom left) or greater than +0.5 (bottom right) are shown in parentheses.

Figure 13 .
Figure 13.Frequency of deviations (observed AODs minus modeled AODs) from the 1 h binned AERONET AOD.The percentages of deviations between −0.05 and +0.05 (−0.10 and +0.10) are shown at the top left (top right), and the percentages less than −0.5 (bottom left) or greater than +0.5 (bottom right) are shown in parentheses.

Figure 15 .
Figure 15.Time series of (a) RMSE, (b) R, (c) MFE, (d) MFB, and (e) IOA for FR (blue) and RA (red), validated by monthly averaged AERONET AODs.The number of AERONET sites used at each time point for these calculations is shown by the thick black line in Fig. 4.

Table 1 .
Configuration of the reanalysis.

Table 2 .
Statistical measures of the NRL-UND MODIS AODs versus FR, FG, and RA AODs.

Table 3 .
Statistical measures of the 1 h binned AERONET AODs versus FR and RA AODs.

Table 4 .
Statistical measures of the monthly averaged AERONET AODs versus FR and RA AODs.

Table 5 .
Statistical measures of the monthly AERONET AODs versus FR and RA AOD at megacity and mountain AERONET sites.China (28.365 • N, 86.948 • E), CASLEO, Argentina (31.799 • S, 69.306 • W), and Mauna_Loa, Hawaii (19.539 • N, 155.578 • W) (Table 5).All of these sites have large MFE and MFB values; when MFE is equal to MFB, the modeled AODs are always larger than the observed AODs.
(Sessions et al., 2015)averaged AERONET and modeled AODs.The AERONET, FR, and RA AODs are shown with black, blue, and red lines, respectively.Gray shading denotes standard deviations of the AERONET AODs.Monthly averages of NRL-UND MODIS AODs over each AERONET site are shown by circles.RMSE, R, MFE, MFB, and IOA values for RA are also shown.Locations of sites are shown by red circles in Fig.3.Beijing, Capo_Verde, Cart_Site, Chiang_Mai_Met_Sta, GSFC, Kanpur, Moldova, Ragged_Point, and Rio_Branco are included among the sites selected for validation of the ICAP-MME(Sessions et al., 2015).