Improved Satellite Retrieval of Tropospheric NO 2 Column Density via Updating of Air Mass Factor ( AMF ) , Part I : Case Study of Southern China

Improving air quality and reducing human exposure to unhealthy levels of airborne chemicals are important global missions, particularly in China. Satellite remote sensing offers a powerful tool to examine regional trends in NO2, thus providing a direct measure of key parameters that strongly affect surface air quality. To accurately resolve spatial gradients in NO2 concentration using satellite observations and thus understand local and regional aspects of air quality, a priori input data at sufficiently high spatial and temporal resolution to account for pixel-to-pixel variability in the characteristics of the land and atmosphere are required. In this paper, we adapt the Berkeley High Resolution product (BEHR v3.0A, v3.0B and v3.0C) and meteorological outputs from the Weather Research and Forecasting (WRF) model to describe column NO2 in southern China. The BEHR approach is particularly useful for places with large spatial variabilities and terrain height differences such as China. We retrieved tropospheric NO2 vertical column density (TVCD) within part of southern China, for four seasons of 2015, based upon satellite datasets from Ozone Monitoring Instrument (OMI). Retrieval results are validated by comparing with MAX-DOAS tropospheric column measurements conducted in Guangzhou. BEHR retrieval algorithms are more consistent with MAX-DOAS measurements than OMI-NASA retrieval, opening new windows into research questions that require high spatial resolution, for example retrieving NO2 vertical column and ground pollutant concentration in China and other countries.


Introduction
Air pollution is a worldwide issue.Atmospheric pollutants can have devastating effects on human health and ecosystems, including contributions to long term human diseases such as cancer and cardiovascular malfunctions [1,2].Poor air quality also causes deterioration of forest ecosystems and direct damage to trees and plants [3,4].In particular, nitrogen oxides (NOx = NO + NO2) are trace gases that have both anthropogenic and natural sources, such as the burning of fossil fuels, emissions from motor vehicles, lightning and biomass burning [5][6][7][8].These trace gases have important effects on air quality and the production of ozone and aerosols within the atmosphere [9,10].NOx regulates surface ozone levels, initiates acid rain, and contributes to enhanced aerosol formation (e.g.NH4NO3) and affects the global climate through chemistry of the greenhouse gases, ozone and methane [11,12].
In recent years, the Chinese government has implemented control policies and plans to reduce the amount of trace gases emitted by fossil fuel burning and vehicle emissions.These plans include clear targets and goals for the reduction of energy consumption, and emission of trace gases (such as NO2 and SO2) in Five-Year Plans (FYPs) 12 (2011-2015) and 13 (2016-2020) [13,14].These documents include targeted reduction percentages for each specific pollutant at both national and provincial levels.However, it is difficult to assess the effectiveness of control measures practiced at the provincial level, especially for NO2.This is mainly due to the lack of a complete ground monitoring network for the whole of China.Only Beijing, Guangzhou, Nanjing, Shanghai and a few more developed cities in China had ground monitoring stations for air quality research studies before 2013 [15], so there are no long-term and large-scale monitoring datasets from which to assess changes in pollutant emissions throughout previous years.To conduct in-situ NOx and NO2 measurements, some Chinese local governments follow the example of many countries and have used chemiluminescence techniques for measuring ground surface NO2 concentrations since 2010 [16,17], and the number of monitoring stations within the "Guangdong-Hong Kong-Macao-Pearl River Delta Regional Air Monitoring Network" has increased from 13 to 18 since its establishment in September 2014 [18].Satellite datasets and remote sensing techniques monitor the column density of pollutants within the lower troposphere, and ground pollutant concentrations are derived from inverse formulations [19] and nonlinear statistical models [20].This data provides a long record (more than 20 years) that can be used to evaluate changes in emissions and chemistry.
Emission changes occur on local and regional scales with gradients of change that are often small compared to standard approaches to interpreting satellite observations.In this paper, we combine detailed high resolution land-use of China with similarly high resolution meteorological and chemical information derived from the Weather Research and Forecasting (WRF) and Community Multiscale Air Quality (CMAQ) models to develop high spatial resolution NO2 columns, using the Berkeley High Resolution (BEHR) [21][22][23] approach in deriving an air mass factor (AMF) and applying the algorithm to the main part of southern China.Then, by combining the adjusted AMF spatial patterns with available satellite datasets and remote sensing techniques, we retrieve the tropospheric NO2 VCD for the different seasons of 2015.We compare different versions of BEHR and OMI-NASA retrieval results with MAX-DOAS tropospheric column measurement datasets obtained in Guangzhou, and find that BEHR VCDs show better agreement with the independent MAX-DOAS measurements.
This article is structured as follows: Section 2 includes the description of satellite retrieval.Section 3 describes our study areas, the model and satellite datasets and the set-up of a common grid for easier manipulation and comparison of the big datasets required for our analysis.Section 4 presents the results and comparison of VCDs based upon different model outputs as well as satellite products and retrieval algorithms.Section 5 includes more details regarding improvement and correction of BEHR v3.0B to v3.0C, and provides a statistical comparison of these retrieval methods.Raw measurements in Guangzhou are used for validation in Section 6, which verifies that the revised retrieval algorithm is an improvement compared with previous methods.Numerical accuracies of each NO2 satellite retrieval method are discussed and highlighted for further research studies.Finally, Section 7 summarizes the main findings and concludes with some possible future research direction.

Principles of NO2 remote sensing
In 2004, OMI was launched with satellite EOS-Aura to measure solar radiation backscattered by the Earth with the help of 2D charge-coupled device (CCD) detectors [24,25].It passes the equator every day around 13:30 LT, with a spatial resolution of 13 × 24 km 2 at its nadir and a spectral resolution of 0.5 nm.After raw datasets from satellites are collected, the reflected sunlight is fit using Differential Optical Absorption Spectroscopy (DOAS) [25], based on Beer-Lambert law, that the light intensity at a specific wavelength can be expressed as in Equation (1).
The absorbance can be separated into a component that varies slowly with wavelength ( ()), and another varies quickly with wavelength (the exponential part), with  being the absorption cross section of single trace gas i (temperature and wavelength dependent), and   () stands for the concentration element in trace gas i.By using slowly varying component as baseline, the slant column density (SCD) can be computed for individual trace gases by multistep removal of instrumental noise and effects like Rayleigh scattering [26].Afterwards, the stratospheric contribution of SCD is removed, with methods in [24,27] as the two main approaches, and AMF is estimated from the a priori NO2 profile, scattering weights and radiance cloud fraction obtained from the regional model, chemistry transport model and satellite products.AMF is a necessary quantity that converts tropospheric SCD into vertical column density (VCD) to account for atmospheric light scattering [28].The AMF is calculated as the ratio of modeled slant to vertical column densities (shown in Equation ( 2)), which is then used to convert the measured slant column density to a vertical column density by Equation ( 4).
Here SCDtrop represents the tropospheric slant column density, W is a scattering weight from TOMRAD, NO2p represents the a priori NO2 profile from the CMAQ simulation, and we assume "0" and  to be the surface layer and tropospheric height respectively.

BEHR algorithm for NO2 column retrieval
This paper adopts NASA Standard Product v3.0 [30] and a modified version of the BEHR product [21].The methods of spectral fitting and stratospheric subtraction in BEHR algorithm are essentially the same as for OMI-NASA retrieval, as described in Section 2.1.However, the BEHR algorithm uses higher resolution NO2 profiles, surface reflectance and terrain pressure data comparing with the NASA algorithm [31].For this domain, these NO2 profiles are simulated by WRF-CMAQ, with its configuration listed in Table 1 of Section 3.2.Thus, NO2 components obscured by clouds in the atmosphere are also included and an accurate spatial tropospheric VCD map is obtained [32].Throughout this paper, we only use pixels that have Cloud Fraction ≤ 0.2, tropospheric AMF > 10 -6 , and satisfy Quality Flags constraints [31], and reject all low quality pixels.

Calculation of AMF and Tropospheric VCD
The BEHR algorithm is described in detail in [21].Briefly, AMF is calculated as an integral contributed by scattering weight and shape factor, as shown in Equations ( 5) and (6), where p0 refers to the surface pressure for the clear-sky AMF and the cloud pressure for the cloudy-sky AMF, ptrop represents tropopause pressure, and () =  ()(1 − 0.003() + 0.66) represents the scattering weight that depends on several meteorological factors;  () obtained from look-up tables, () taken from WRF, and () is the shape factor that depends on NO2 vertical profile ().
For cloudy-sky conditions, the cloud radiance fraction f is taken into consideration.The tropospheric AMF for general atmospheric conditions is written as the linear combination of clear-sky AMF and cloudy-sky AMF, shown in Equation (7).AMF = ∫ ()() (5) Combining Equations ( 5) and (7), and let  () and  () be the scattering weights of the clear and cloudy parts of the concerned pixel, the BEHR AMF can be evaluated by Equation (8), based on Section 2 of [21], for which the parameters are well-defined in previous paragraphs.
The BEHR retrieved tropospheric VCD is calculated by dividing the OMI-NASA tropospheric SCD by AMFBEHRtrop, as shown in Equation (9) [30,33].
2.2.2 Differences in BEHR v3.0A, v3.0B and v3.0C Retrieval The AMF calculation requires several a priori inputs, including surface pressure and reflectivity, and NO2 vertical profiles.Compared to the OMI-NASA product, the values used in the BEHR AMF calculation are ingested at higher resolution and averaged to the OMI pixels: the surface reflectivity is derived from the MODIS MCD43D products at 1 km resolution; terrain elevation is calculated from the 1 km GLOBE terrain database [21].Vertical profiles of NO2, temperature, pressure, and geopotential are obtained from model simulations; for the BEHR product described in this paper covering the southern China domain, these profiles are simulated at 9 km resolution with the WRF-CMAQ model.BEHR v3.0A used a fixed 200 hPa tropopause and calculated surface pressure from the GLOBE terrain height, assuming a 7.4 km scale height.However, in BEHR v3.0B, the calculation of surface pressure makes use of the method described in [34] to combine the GLOBE terrain height with WRF surface pressure, while the tropopause pressure is calculated based on the WMO thermal tropopause definition [35], with temperature profiles obtained from WRF model.Discontinuities in the tropopause pressure are masked and filled in by an iterative interpolation algorithm, as the tropopause pressure is assumed to vary smoothly in space.Finally, the tropopause height is calculated from interpolated thermal tropopause pressure [32].The differences between BEHR v3.0A and v3.0B retrieval are described in details in [21] and Changelog in [36].
The implementation of calculating tropopause pressure in BEHR v3.0B has been found to be sensitive to the spacing of the levels within temperature profiles, therefore we test an experimental branch of BEHR, v3.0C that modifies how the tropopause pressure, which defines the upper limit for AMF calculation, is computed.BEHR v3.0C uses the tropopause pressure available in the NASA Standard Product [24,30].

CMAQ Tropospheric VCD Simulation
A priori NO2 profiles are necessary for running BEHR.A CMAQ modeling system driven by WRF outputs simulates the NO2 concentration (in ppbv) of each pixel within each of the 26 vertical layers to provide the corresponding profiles.CMAQ simulation results are based on corresponding WRF outputs, together with the Sparse Matrix Operator Kernel Emissions (SMOKE) modeling system, which provides emission data and initial and boundary conditions [37,38], however lightning NOx emission is not well configured in CMAQ simulation.Each of these 26 vertical layers has their respective sigma values, pressure, height and thickness, so we can calculate the simulated tropospheric VCD based on these data.If we assume that the spatial resolution of our retrieved domain is 9 × 9 km 2 , Ci is the concentration (ppmV) of layer i that consists of the sum of the concentrations of all valid pixels within the i th layer, hi is the thickness (in m) of the i th vertical layer, and the room temperature is (273.15+ 25) K, the tropospheric NO2 VCD of a particular timeslot can be calculated by Equation (10), with the aid of the ideal gas equation (PV = nRT).Here 12.187 represents the inverse of Universal Gas Constant (R), 46.01 is the molecular weight (g/mol) of NO2, 7.64 × 10 is the factor that converts μg/m into molecules/cm 2 , and the term ∑ 9000 ℎ  represents total tropospheric NO2 column concentration in μg.

Region of Interest
We focus on OMI satellite data for the part of southern China (including Dongguan, Foshan, Guangzhou, Huizhou, Jiangmen, Shenzhen and Zhuhai), with longitude from 109°E -117.5°E,latitude from 19.5°N to 25.5°N, as shown in Figure 1.Most of these cities are located in the coastal region of southern China, which is characterized by complex topography and a variety of geographical features, including hills and islands.In particular, Shenzhen and Hong Kong are close to the sea and have steep gradients in land surface [39], so sea-land breezes cause some changes in meteorological and climatic conditions in these areas [40,41].The region shown in Figure 1 is intentionally chosen to investigate and test the accuracy and validity of satellite remote sensing algorithms and chemical transport models in retrieving tropospheric NO2 VCD, under conditions where high spatial resolution is expected to be important.

Datasets
BEHR is an algorithm that focuses on achieving high spatial resolution from satellite remote sensing of NO2.The required inputs to drive the BEHR algorithm in retrieving the vertical column of NO2 are listed in Table 1 below.For the NO2 slant column, we adopt the raw OMI datasets [42]; for albedo, we use the Moderate Resolution Imaging Spectroradiometer (MODIS) Bidirectional Reflectance Distribution Function (BRDF) reflectance from MCD43D07-09 [43] (with spatial resolution of 0.05° × 0.05°, and calculation of 16-day average every 8 days).For terrain pressure, elevation data from the Global Land One-Kilometer Elevation (GLOBE) database (with a resolution of 1 km x 1 km) [44] is averaged within the OMI pixel bounds to yield the pixel elevation.In BEHR v3.0B and onwards, base and perturbation geopotential profiles from WRF are also necessary to calculate the surface pressure, based on the approach outlined in [34]. 1 Spatial resolution: 9 km; 2 scaled from 0 to 1; 3 spatial resolution: 13 × 24 km 2

Grid Formation and Decomposition
As shown in Table 1, the different datasets mentioned in Section 3.2 have different array structures.After completing BEHR (v3.0A-v3.0C)run, AMF and tropospheric VCD will be 2D arrays with different grid points in each dimension during every satellite-passing hour.For consistent projection and fair comparison of spatial plots, we create a rectangular grid within our domain of interest (from 108°E to 118°E and 18°N to 26°N), with a mesh size of 0.1°.In this grid, we decompose the area into 8181 small grids, each with an area of 0.1° × 0.1°.For all satellite datasets, we first project every available pixel onto this grid hour by hour, by weighting the contribution of each pixel to the grid cell based on amount of areas that overlaps, then take the temporal average within each prescribed period.We then construct spatial plots based on the longitude, latitude and information of concerned attribute.

Results
This section provides the retrieval results of tropospheric NO2 column density obtained from OMI-NASA retrieval (Section 2.1), BEHR retrieval (Section 2.2) and WRF-CMAQ simulation (Section 2.3).Corresponding AMF spatial distributions within Southern China for OMI-NASA and BEHR retrieval for all four months are presented as part of the supplementary material of this paper (Figures S1 and Figures S2 in Supplementary Materials).

Comparison of OMI-NASA, BEHR and WRF-CMAQ VCDs
The AMF obtained by OMI-NASA retrieval generally has higher values (ranging from 0.5 to 0.8) for most pixels when comparing with BEHR v3.0C retrieval, as shown in Figure S1 and S2 in Supplementary Materials.By Equation ( 9), we can obtain the tropospheric NO2 VCD using different means of satellite retrieval and simulation.
Figure 2 shows the monthly average tropospheric NO2 VCD based upon the OMI-NASA retrieval.The tropospheric VCD in winter (represented by Jan 2015 in Figure 2(a)) is the highest out of all four seasons (many pixels exceed 10 16 molecules/cm 2 ), especially cities lying on coastal region such as Dongguan (DG), Foshan (FS), Guangzhou (GZ), Shenzhen (SZ) and Zhongshan (ZS), followed by spring and autumn (represented by Apr and Oct 2015), ranging from 3 to 9 × 10 15 molecules/cm 2 .The average tropospheric VCD in summer is the lowest, at around 2 to 6 × 10 15 molecules/cm 2 .As mentioned in Section 2.2, BEHR makes use of higher spatial resolution a priori NO2 profiles for retrieval, in combination with meteorological variables within southern China from WRF simulations.Therefore, it is believed to capture differences between retrieved NO2 VCDs for neighboring pixels well, especially when a pixel representing a rural area is situated next to a pixel that represents an urban region.Within many provinces in China, for example in Guangdong, the land-use pattern can change abruptly, in particular, the spatial variability of Guangzhou (GZ) is illustrated in [45].Figure 3 shows the spatial distribution of BEHR v3.0C retrieved tropospheric NO2 column density within 4 retrieval periods.Further statistical analysis and correlation between v3.0B and v3.0C retrieval are provided and discussed in Section 5.2.The spatial pattern looks quite similar to OMI-NASA retrieval (in Figure 2), except that BEHR captures some relatively higher values of NO2 VCD in the east and north of our domain, for example Meizhou (MZ) and Shanwei (SW) in the east, and Shaoguan (SG) in the north of Guangzhou.In BEHR retrieval, these areas show significantly higher amounts of NO2 VCD compared with OMI-NASA retrieval.Interestingly, BEHR shows greater NO2 VCDs in April than October, whereas OMI-NASA shows the reverse pattern.4.2.BEHR v3.0C v.s.OMI-NASA NO2 VCD BEHR v3.0C was developed during the course of writing this manuscript to address deficiencies in earlier versions that were only tested for North American conditions.We compare the BEHR v3.0C tropospheric NO2 VCD with that of the OMI-NASA product by calculating the percentage difference between the two, where percentage difference = (BEHR VCD -OMI-NASA VCD) / OMI-NASA VCD × 100%.Corresponding spatial plots are illustrated in Figure 5. From these plots, we find that in cities located in coastal regions, such as Guangzhou and Shenzhen, BEHR v3.0C retrieval gives lower NO2 VCD in January and October 2015, while higher NO2 VCD is obtained in April and July 2015.However, the trend is opposite for Chaozhou and Shanwei.For cities that are located in inland areas like Hezhou, Qingyuan and Shaoguan, as well as Haikou (near Hainan Province), BEHR retrieval gives higher NO2 VCD than the OMI-NASA product in all four periods.As shown in Figure S1 and S2 in Supplementary Materials, BEHR v3.0C generally has lower AMF comparing with OMI-NASA, possibly due to lack of lightning NOx information in WRF-CMAQ run, which produces the a priori NO2 profile for conducting BEHR run.Nevertheless, comparing with WRF-CMAQ simulation in Figure 4, we can conclude that BEHR retrieval is able to better capture the spatial variability in southern China, for example the changes of NO2 VCD or NOx emission within neighboring cities in Guangdong province.Moreover, having a priori inputs and satellite products of higher spatial resolution generally lead to higher VCD in cities, but the situation is more variable in China due to its highly heterogeneous land use.Many cities in southern China contain both rural and urban areas, and it is believed that urban areas with higher transportation flow rate and more industrial activities will emit higher amount of NO2, while rural areas with less commercial activities produce lower NO2 emissions.Complex mixture of urban and rural areas somehow explain why BEHR VCDs are not uniformly greater than NASA VCDs in cities.

Correction of BEHR v3.0B and Improvement to v3.0C
In the process of retrieval, we identified a flaw in BEHR v3.0B, then we developed v3.0C to repair this flaw.As described below, the issue appears to be associated with tropical conditions that were not tested during the development of v3.0B and doesn't affect conclusions related to application of v3.0B to North America [21,23].To understand and optimize the approach to high resolution retrievals, we compare different retrieval methods.In Section 5.1, we describe the abnormally large VCDs encountered in BEHR v3.0B retrieval within July 2015, and provide reasons for causing such discrepancies and ways of dealing with such issue.In Section 5.2, we focus on comparing BEHR v3.0B and BEHR v3.0C retrieval, and provide comments regarding the potential effects after tropopause pressure formulation is altered.

Causes of abnormally large VCDs in BEHR v3.0B
In previous literature based on OMI-NASA satellite retrieval [22,46], the reported tropospheric NO2 VCD within southern China (see Section 3.1) is within the range of 10 15 to 10 16 molecules/cm 2 .After updating AMF based on BEHR algorithm, all output files of BEHR v3.0A, v3.0B and v3.0C give reasonable numerical values of tropospheric NO2 VCD throughout our interested domain in January, April, and October 2015, respectively.However, some pixels of several dates in July 2015 have unreasonably high numerical values of tropospheric NO2 VCD (>10 17 molecules/cm 2 ).Table 2 provides details of occurrences of these unexpectedly high values after regridding is conducted in BEHR v3.0B retrieval.These abnormally large VCDs occur in pixels where both  <  and the cloud radiance fraction () is large.In BEHR v3.0B and v3.0C, if  <  , the cloudy-sky AMF is set to zero, representing a case where there is no above-cloud NO2.When the cloud fraction is also of large quantity, this will make the overall AMF (AMF ) very small by Equation (7) .Since the VCDs are inversely proportional to the AMF by Equation ( 9), this results in extremely large VCDs.Fundamentally, this error occurs when most of the NO2 column are obscured by clouds and the AMF calculation must attempt to reconstruct the total tropospheric NO2 VCD with limited observational information.Such error is more pronounced in the southern China domain than the US domain, mainly because in southern China domain, the tropopause pressure calculated from WRF temperature profiles in BEHR v3.0B is biased high by 180 hPa comparing with the tropopause calculated by NASA, whereas the difference is only ~50 hPa in the US domain.This high bias increases the frequency of pixels with  <  .This bias in the tropopause was traced to the coarse spacing of the WRF temperature profiles in the upper troposphere for the southern China domain, and can be corrected in BEHR v3.0B retrieval by removing any pixels for which  <  , however BEHR v3.0C replaces the WRF-derived tropopause pressure with the NASA-calculated one.Figure 6 shows the spatial plots of average tropospheric NO2 VCD in July 2015 before (Figure 6(a)) and after pixel filtering (Figure 6(b)) is adopted.Figure 6(c) shows the BEHR v3.0C retrieval result.From the figures, we conclude that both of the aforementioned approaches successfully remove pixels with anomalously large VCDs.

BEHR v3.0B v.s. BEHR v3.0C
Figure 7 shows the correlation plot between BEHR v3.0B and BEHR v3.0C for four different months in 2015.For Figure 7(c), we conduct filtering (Section 5.1) and eliminate 56 problematic pixels in BEHR v3.0B dataset and obtain the correlation plot for July 2015.Table 3 presents the statistical parameters that assess the relationship between NO2 VCD for the two BEHR retrievals.Readers can refer to S3 of the Supplementary Materials for details of these statistical parameters.From the relative position of data points with respect to y = x (blue line) shown in Figure 7, together with the slope of best-fit line for the four months, we conclude that BEHR v3.0B retrieval generally gives a higher amount of NO2 VCD than BEHR v3.0C.This happens in all 4 retrieval periods in 2015.For July 2015, even after filtering out problematic pixels in BEHR v3.0B, around 85% of data points lie below the blue line (as shown in Figure 7(c)), which indicates that tropopause pressure difference between satellite products is a key factor towards tropospheric NO2 VCD retrieval.Since the thermal tropopause calculation overestimates the tropopause pressure because of coarse spacing of temperature profiles in the UT, NASA tropopause is adopted in BEHR v3.0C to better capture the accurate height of tropical tropopause within southern China.This means that the NO2 VCD profile obtained by v3.0C retrieval is supposed to represent the actual situation better.Thus, we conclude that there is slight over-estimation in most circumstances when the v3.0B retrieval method is used.
Moreover, the R-value in winter and spring is slightly higher than that in summer and autumn; thus, the corresponding RMSE values in these seasons are lower.One possible reason is the changing dispersion rate of pollutants in different seasons.  1 1213.30.9958 2.55 ×10 14   Apr 2015 m = 0.926 c = -3.623×10 13  1809.90.9983 1.16 ×10 14 July 2015 1 (Before filtering) m = 0.123 c = 1.025 ×10 16   (After filtering) m = 0.7519 c = 5.029×10 12    Peer-reviewed version available at Remote Sens. 2018, 10, 1789; doi:10.3390/rs10111789 1 Filtering of pixels is conducted in July 2015 to eliminate problematic pixels, and 8125 valid pairs of data points remain for analysis

Discussion
In this section, we compare the satellite retrieval algorithms with available raw measurement datasets and obtain corresponding uncertainties and errors for each method, thus providing a source to assess the applicability of individual mechanisms in large-scale satellite retrieval all over the world.

MAX-DOAS Validation in Guangzhou
For validation purposes, we compare satellite retrieval results for all algorithms with MAX-DOAS tropospheric measurements in Guangzhou.Within our retrieved domain, Guangzhou is the largest city.It has enormous commercial and industrial activities, and thus high NOx emission [47], so obtaining time series of both ground NO2 concentrations and tropospheric NO2 VCD is extremely important for imposing future control policies.Here we adopt the tropospheric measurements conducted in [48] for comparison, for which a MAX-DOAS system (composed of CCD spectrograph) was installed on the roof of a nine-story building in Guangzhou Institute of Geochemistry of the Chinese Academy of Sciences (GIG-CAS) from March 2015 to March 2016 (Geographical location: latitude: 23°8'54''N, longitude: 113°21'32''E) [48].A spectrograph with spectral resolution of 0.35 nm captured radiation of wavelength 300-450 nm, as well as measuring sky radiance at different elevation angles [48,49].The tropospheric NO2 VCD measurement figures are based on homogeneous spatial distribution of NO2 within the lower troposphere along every azimuthal direction.Within our retrieval period, comparisons can only be made for April and July 2015, due to damage to the optical system [48].We first consider all valid daily measurements that are + or -0.5 hours of the satellite retrieval period, for example if the satellite retrieval period on a specific day is 06:00 and 07:00 UTC, then all MAX-DOAS measurements conducted from 05:30 -07:30 UTC will be taken into account.Based on all these available measurements, daily average is calculated, which provides us precise MAX-DOAS tropospheric NO2 VCD measurements for comparison purposes.Next, for our satellite-retrieved datasets, upon projection onto our preset grid (Section 3.3), we refer to the grid that represents Guangzhou (latitude: 23.1°N, longitude: 113.3°N) and obtain corresponding OMI-NASA and BEHR satellite-retrieved tropospheric NO2 VCDs for each day.Figures 8 and 9 show the comparison of MAX-DOAS, OMI-NASA, BEHR (v3.0C) retrieved NO2 VCD for April and July 2015, together with their corresponding error bound.Here we adopt ±15% of tropospheric NO2 column density values as its error bound.As shown in Figures 8 and 9, the OMI-NASA product underestimates the tropospheric NO2 VCD in Guangzhou throughout the whole retrieval period, although its time trend is in line with MAX-DOAS measurements.In contrast, both the magnitude and temporal pattern of the BEHR VCDs is similar to that observed by the MAX-DOAS instrument.This is true for BEHR v3.0A and v3.0B as well (not shown).Within the satellite passing hours of some dates (e.g. 10 April, 22 April, 5 July, 18 July, 20 July, 30 July), BEHR v3.0A retrieval gives relatively higher NO2 VCD than BEHR v3.0B, but the situation is reversed or no different for most other dates.In Guangzhou, BEHR v3.0B and v3.0C retrieval results are basically the same, where BEHR v3.0B gives slightly higher NO2 VCD than v3.0C in most of the dates, especially in the main peaks.

Numerical Uncertainties of each satellite retrieval method
We investigate the correlation between each satellite retrieval with MAX-DOAS measurements in Figure 10.The 30 data points represent 12 dates in April 2015 and 18 dates in July 2015 (as shown in the data points of Figure 8 and 9).In each of these plots, the blue line represents y = x.From Figure 10(a), we can see that there will be a large underestimation if we use the OMI-NASA product to retrieve NO2 VCD within China.Based upon Figure 10(b)-(d), we realize that versions of the BEHR algorithm may over-or under-estimate the VCDs comparing with the MAX-DOAS measurements, but do not systematically underestimate the NO2 VCD.However, due to updating of AMF, NO2 VCDs of some dates have become over-estimated in BEHR v3.0A retrieval.One possible reason is that we assume a fixed tropopause pressure, which is not true in general.By using variable tropopause top pressure, there is a clear improvement as shown in Figure 10(c) (BEHR v3.0B) and Figure 10(d) (BEHR v3.0C),where the R-values of the data pairs are 0.9338 and 0.9839 respectively, which are significantly higher than that in Figure 10(a) and 10(b), which have R-values of 0.7644 and 0.8468, respectively.The improvement can also be seen in the decrease in RMSE, from 8.961 × 10 15 molecules/cm 2 in the OMI-NASA retrieval to 3.992 ×10 15 molecules/cm 2 in BEHR v3.0B retrieval, then a further decrease to 2.083 ×10 15 molecules/cm 2 for BEHR v3.0C retrieval.The improvement in RMSE between version 3.0B and 3.0C is especially impressive and beyond what we expected to find.The small RMSE suggests that even if there is an absolute bias due to omitting lightning from the model used to calculate the prior, the bias is likely uniform across large swaths of the domain.
As noticed in Figure 10(c) and (d), most valid data points are closer to the y = x line for BEHR v3.0B and v3.0C retrieval, demonstrating that the BEHR NO2 VCDs are more accurate than the NASA VCDs in Guangzhou.This also suggests that the BEHR VCDs will better represent the NO2 concentration in other parts of China as well.Based on the data values, we compute the percentage difference between each satellite retrieved NO2 VCD and MAX-DOAS measurement value by Equations ( 11) and (12), respectively.
Here the absolute sign is to make the percentage difference non-negative for comparison purposes.A bar chart that shows the distribution of percentage differences is shown in Figure 11, together with the assumption that MAX-DOAS measurements give accurate results, we can conclude that for OMI-NASA product, degree of accuracy in southern China is 30% to 90%, while ~86.7% of dates for BEHR v3.0C are within 20% of MAX-DOAS VCD.
Future work will expand this product to a domain covering more of East Asia.It will be important to validate the retrieved VCDs in other cities throughout the larger domain to ensure that the agreement we see here is also representative of the accuracy of the BEHR retrieval in the larger domain.In particular, with more advanced measurement technologies for setting up MAX-DOAS systems within different cities in China, we can include more valid data points in our statistical analysis.This will enable more thorough assessments of different satellite retrieval methods, and estimations of spatial and temporal uncertainties within different cities around the world, and within different months or seasons throughout a year.Peer-reviewed version available at Remote Sens. 2018, 10, 1789; doi:10.3390/rs10111789

Conclusions
We adapt the BEHR approach to update AMF distribution within southern China and conduct NO2 retrievals.In particular, the tropospheric NO2 VCDs in four selected months within 2015 are retrieved.We take into account of the meteorological outputs from WRF model, and identify improvements to the BEHR algorithm based on considering atmospheric conditions and characteristics of southern China.We develop an improved algorithm v3.0C and successfully show that this version is highly correlated with tropospheric MAX-DOAS observations.Moreover, BEHR v3.0C has relatively small RMSE compared to the OMI-NASA retrieval, BEHR v3.0A and v3.0B, based on our MAX-DOAS datasets in Guangzhou.Such statistical improvement will be further verified when we have more tropospheric NO2 column measurements available from the surface.Future work will seek to expand this retrieval to larger domains in Asia.Obtaining tropospheric NO2 distributions of sufficiently high spatial resolution worldwide is crucial to estimate spatial distribution of ground-level NOx emission and NO2 concentrations through data assimilation techniques.Our methods advance the BEHR approach to obtaining an accurate and precise worldwide NO2 map with high spatial resolution in the foreseeable future.Nevertheless, most of the evaluation in this paper is within one location and comparing variations in time.The strength of BEHR retrieval is its variations in space at a single time, therefore one of our future goals is to conduct more comprehensive evaluation of the spatial pattern with high density surface observations.

Figure 1 .
Figure 1.This study's region of interest (the main part of southern China)

Figure 2 .
Figure 2. OMI-NASA monthly average tropospheric NO2 VCD in southern China.(a) Jan 2015; (b) Apr 2015; (c) Jul 2015; (d) Oct 2015.The units of the figures are molecules/cm 2 , ranging from 0 to 1.4 × 10 16 molecules/cm 2 .Here we only consider pixels that have Cloud Fraction ≤ 0.2, tropospheric AMF > 10 -6 and satisfy Quality Flags constraints stated in Section 2.2.Cities mentioned in Section 4.1 are labeled in the spatial plots.

Figure 4
Figure 4 provides the spatial distribution of tropospheric NO2 VCD obtained by WRF-CMAQ simulation, based upon the formulation discussed in Section 2.3.

Figure 7 . 1 . 3 .
Figure 7. BEHR v3.0C vs. BEHR v3.0B.(a) Jan 2015; (b) Apr 2015; (c) Jul 2015 (based on filtered datasets in BEHR v3.0B);(d) Oct 2015.The blue lines represent the line y = x, while the data points represent the pixel value at each small grid based upon the two retrievals; Both x and y axis are tropospheric NO2 VCD, units are molecules/cm 2 ; the color bar indicates linearly spaced points from 0 to 1.Table3.Comparison of statistical parameters for BEHR v3.0C (y-axis) and v3.0B (x-axis) retrieval (Least Square Fitting by assuming linear regression holds)

Figure 8 .
Figure 8. Tropospheric NO2 VCD in April 2015 obtained through MAX-DOAS measurements (Purple Triangles), and satellite retrieval: OMI-NASA (Red Squares), BEHR v3.0C (Green Circles) for 12 dates for which all datasets have available information.The error bound (±15%) indicates the uncertainty estimates of each measurement or retrieval result.

Figure 9 .
Figure 9. Tropospheric NO2 VCD in July 2015 obtained through MAX-DOAS measurements (Purple Triangles), and satellite retrieval: OMI-NASA (Red Squares), BEHR v3.0C (Green Circles) for 12 dates for which all datasets have available information.The error bound (±15%) indicates the uncertainty estimates of each measurement or retrieval result.

Table 1 .
Attributes of datasets used for satellite data retrieval

Table 2 .
Dates and pixels that return unusually high tropospheric VCDs after running BEHR v3.0B

Table 4 .
Comparison of statistical parameters for satellite retrieval algorithms (y-axis) and MAX-DOAS tropospheric measurements (x-axis) (Least Square Fitting by assuming linear regression holds)