Radiative transfer acceleration based on the principal component analysis and lookup table of corrections: optimization and application to UV ozone profile retrievals
- 1Center for Astrophysics, Harvard & Smithsonian, Cambridge, MA, USA
- 2RT Solutions Inc., Cambridge, MA, USA
- 3Department of Atmospheric and Oceanic Science, University of Maryland, College Park, Maryland, USA
- acurrently at: Institute of Environmental Studies, Pusan National University, Busan, Korea
Correspondence: Juseon Bak (firstname.lastname@example.org)
In this work, we apply a principal component analysis (PCA)-based approach combined with lookup tables (LUTs) of corrections to accelerate the Vector Linearized Discrete Ordinate Radiative Transfer (VLIDORT) model used in the retrieval of ozone profiles from backscattered ultraviolet (UV) measurements by the Ozone Monitoring Instrument (OMI). The spectral binning scheme, which determines the accuracy and efficiency of the PCA-RT performance, is thoroughly optimized over the spectral range 265 to 360 nm with the assumption of a Rayleigh-scattering atmosphere above a Lambertian surface. The high level of accuracy (∼ 0.03 %) is achieved from fast-PCA calculations of full radiances. In this approach, computationally expensive full multiple scattering (MS) calculations are limited to a small set of PCA-derived optical states, while fast single scattering and two-stream MS calculations are performed, for every spectral point. The number of calls to the full MS model is only 51 in the application to OMI ozone profile retrievals with the fitting window of 270–330 nm where the RT model should be called at fine intervals (∼ 0.03 nm with ∼ 2000 wavelengths) to simulate OMI measurements (spectral resolution: 0.4–0.6 nm). LUT corrections are implemented to accelerate the online RT model due to the reduction of the number of streams (discrete ordinates) from 8 to 4, while improving the accuracy at the level attainable from simulations using a vector model with 12 streams and 72 layers. Overall, we speed up our OMI retrieval by a factor of 3.3 over the previous version, which has already been significantly sped up over line-by-line calculations due to various RT approximations. Improved treatments for RT approximation errors using LUT corrections improve spectral fitting (2 %–5 %) and hence retrieval errors, especially for tropospheric ozone by up to ∼ 10 %; the remaining errors due to the forward model errors are within 5 % in the troposphere and 3 % in the stratosphere.
Optimal-estimation-based inversions have become standard for the retrieval of atmospheric ozone profiles from atmospheric chemistry UV and visible (UV–Vis) backscatter instruments. This inversion model requires iterative simulations of not only radiances, but also of Jacobians with respect to atmospheric and surface variables, until the simulated radiances are sufficiently matched with the measured radiances. These ozone profile algorithms face a computational challenge for use in global processing of high spatial–temporal resolution satellite measurements, due to online radiative transfer (RT) computations at many spectral points from 270 to 330 nm; it is computationally very expensive to perform full multiple-scattering (MS) simulations with the polarized RT model. To reduce the computational cost, a scalar RT model can be applied together with a polarization correction scheme based on a lookup table (LUT) (Kroon et al., 2011; Miles et al., 2015). Another approach is to carry out online vector calculations at a few wavelengths (Liu et al., 2010) together with other approximations (e.g., low-stream, coarse vertical layering, Lambertian reflectance for surface and cloud, no aerosol treatment). However, the computational speed is still insufficient to process 1 d of measurements from the Aura Ozone Monitoring Instrument (OMI) within 24 h (30 cross-track pixels × 1644 along-track pixels × 14 orbits) with reasonable computational resources. Consequently, only 20 % of the available OMI pixels are processed to generate the operational ozone profile (OMO3PR) product (Kroon et al., 2011), and the spatial resolution is degraded by a factor of 4 to produce the research ozone profile (OMPROFOZ) product (Liu et al., 2010). With the advent of sophisticated inversion techniques and superior spaceborne remote sensing instruments, computational budgets have increased rapidly in recent years. Joint retrievals combining UV and thermal infrared (∼ 9.6 µm) have been investigated to better distinguish between upper- and lower-tropospheric ozone abundances from multiple instruments, e.g., OMI and TES (Fu et al., 2013), OMI and AIRS (Fu et al., 2018), or GOME-2 and IASI (Cuesta et al., 2013). The geostationary Tropospheric Emissions: Monitoring of Pollution (TEMPO) instrument, scheduled for launch in 2022, is specially designed for joint retrievals combining UV and visible (540–740 nm) radiances to enhance the performance of retrievals for ground-level ozone (Zoogman et al., 2017). Moreover, the temporal and spatial resolutions of upcoming geostationary satellite instruments are being improved, leading to a tremendous increase in the data volume to be processed; for example, daily measurements of TEMPO (with ∼ 2000 N/S cross-track pixels 1200 E/W mirror steps 8 times a day) are ∼ 30 times greater in volume than those of OMI. Therefore, accelerating RT simulations is one of the highest priority tasks to assure operational capability. For speed-up, LUTs have often been used in trace gas retrieval algorithms to serve as proxies for RT modeling or to perform corrections to online RT approximations. In recent years, applying neural network techniques and principal component analysis (PCA) to RT computational performance has received quite a lot of attention (e.g., Natraj et al., 2005; Spurr et al., 2013, 2016; Liu et al., 2016; Yang et al., 2016; Loyola et al., 2018; Nanda et al., 2019; Liu et al., 2020).
The goal of this paper is to improve both computational efficiency and accuracy of RT simulations in the OMI ozone profile algorithm (Liu et al., 2010) by combining a fast-PCA-based RT model with two kinds of correction techniques. The application of PCA to RT simulations was first proposed by Natraj et al. (2005) by demonstrating a computational improvement of intensity simulation in the O2A band by a factor of 10 and with ∼ 0.3 % accuracy compared to full line-by-line (LBL) calculations. This scheme has been deployed to the UV backscatter, thermal emission, and crossover regimes and has been extended for the derivation of analytic Jacobians, for vector RT applications, and for bidirectional surface reflectances (Kopparla et al., 2016, 2017; Natraj et al., 2010; Somkuti et al., 2017; Spurr et al., 2013). The RT performance enhancement arises from a reduction in the number of expensive full multiple scattering (MS) calculations; the PCA scheme uses spectral binning of the wavelengths into several bins based on the similarity of their optical properties and the projection to every spectral point of these full MS calculations which are executed for a small number of PCA-derived optical states. In addition to the adaption of a PCA-based RT model for our ozone profile retrieval, we have adopted the undersampling correction from our previous implementation (Kim et al., 2013; Bak et al., 2019); this enables us to use fewer wavelengths for further speed-up without much loss of accuracy. Furthermore, we have developed a LUT-based correction to accelerate online RT simulations by starting with a lower-accuracy configuration (scalar RT with no polarization, 4 streams, 24 layers) and then correcting the accuracy to the level attainable by means of a computationally more expensive configuration (vector RT, 12 streams, 72 layers). The stream value refers to the number of discrete ordinates in the full polar space; thus, for example, the term “12 streams” indicates the use of six upwelling and six downwelling polar cosine discrete ordinate directions. In previous work, PCA-based RT calculations were assessed mostly against LBL calculations, independently from the inverse model. Therefore, the PCA performance is likely to be overestimated in terms of operational capability, because operational algorithms have their own speed-up strategies with many approximations; this is the case for our ozone profile algorithm. As mentioned above, the PCA-based RT model is employed in this work to make forward-model simulations of OMI measurements for the retrieval of ozone profiles. Therefore, we evaluate the operational capability of our retrieval algorithm in terms of the retrieval efficiency as well as the accuracy, and we assess these relative to the current operational implementation.
This paper is structured as follows. Section 2 describes the current forward model scheme and evaluates the approximations made in RT calculations, with the determination of the configuration parameters for accurate simulations. The updated forward model scheme is introduced for the PCA-based RT model in Sect. 3.1, and the two kinds of correction schemes to use fewer spectral samples and a less accurate RT configuration are detailed in Sect. 3.2. The evaluation is performed in Sect. 4, and then we summarize and discuss the results in Sect. 5.
We first describe the current v1 SAO OMI ozone profile algorithm that was implemented in OMI Science Investigator-led Processing Systems (SIPS) to generate the research OMPROFOZ ozone profile product, publicly available at the Aura Validation Data Center (AVDC, https://avdc.gsfc.nasa.gov/index.php?site=1620829979&id=74, last access: 11 March 2021). It employs the OMI UV channel that is divided into UV1 (270–310 nm) and UV2 (310–380 nm). The spatial resolution of UV1 is degraded by a factor of 2 in order to increase the signal-to-noise ratio (SNR) in this spectral region. The full width at half maximum (FWHM) of the instrument spectral response function (ISRF) is ∼ 0.63 nm for UV1 and ∼ 0.42 nm for UV2, with spectral intervals of 0.33 and 0.14 nm, respectively. The total number of OMI wavelengths used in our spectral fitting for ozone profiles is 229, from 270–308 nm (UV1) and 312–330 nm (UV2). The RT model needs to simulate sun-normalized radiances as well as their derivatives with respect to the ozone profile elements and surface albedo. The Vector Linearized Discrete Ordinate Radiative Transfer (VLIDORT) model v2.4 (Spurr et al., 2008) was employed as a forward model in the v1 OMI ozone profile algorithm (Liu et al., 2010) implemented at SIPS. We have updated VLIDORT to the latest version v2.8 for this study as well as in the PCA-VLIDORT described in Sect. 3. Note that there is little difference between v2.4 and v2.8 in terms of simulation accuracy. The RT simulation is iteratively performed to ingest the atmospheric and surface variables adjusted through the physical fitting between measured and simulated spectra and simultaneously the statistical fitting between the state vector and the a priori vector. The retrieval is optimized within typically 2–3 iterations (up to 10 is permitted). The vertical grids of the retrieved ozone profiles in 24 layers are initially spaced in log (pressure) at (in atm, 1 atm = 1013.25 hPa) for 0 (surface), 23 (∼ 55 km), and with the top of atmosphere set for P24 (∼ 65 km). Each layer is thus approximately 2.5 km thick, except for the top layer (∼ 10 km). A number of RT approximations have already been applied in the current forward model to speed up the processing. In the remainder of this section, the current forward model scheme is described, with its flowchart depicted in the left panel of Fig. 1. An error analysis is performed for optimizing the RT model configuration to maximize the simulation accuracy.
In the first step, we select 93 wavelengths with variable sampling intervals, 1.0 nm below 295 nm, 0.4 nm from 295–310 nm, and 0.6 nm above 310 nm. The number of these wavelengths is smaller than the OMI native pixels (229 from 270–330 nm) by more than a factor of 2. The online radiative transfer model is run to generate the full radiance spectrum (single + multiple scattering) at these wavelengths in the scalar mode, with eight streams and a Rayleigh atmosphere divided into 25 layers – a grid that is similar to that for the retrieval, except for the top layer (∼ 55 to 65 km), which is further divided into two layers. In step 2, the scalar calculations done in step 1 are corrected using the online vector calculation at 14 wavelengths (visually shown with the vertical lines in Fig. 3b). In step 3 individual calculations are interpolated into 0.05 nm intervals with the undersampling correction, and the result is finally interpolated/convolved into OMI native grids in step 4.
Figure 2a shows the reference spectrum where Gaussian smoothing to 0.4 nm is applied to LBL calculations at the sampling rate (0.01 nm) of the ozone cross sections (Brion et al., 1993), which is used to evaluate the approximation errors related to undersampling. Figure 2b illustrates that LBL calculations are required to be performed at intervals of 0.03 nm or better. The undersampling correction applied in step 3 allows relaxation of the sampling rate without loss of the accuracy. This correction is based on the adjustment of the radiance due to the difference of the optical depth profiles between fine (λh) and coarse (λc) spectral grids as follows:
where is the weighting function with respect to the optical depth profiles and for gas absorption and Rayleigh scattering, (the number of atmospheric layers). However, as shown in Fig. 2c, the sampling rates (1.0, 0.4, 0.6 nm) used in the v1 forward model are too coarse to be corrected and hence are decided to be 0.3 and 0.1 nm in the v1 forward model. Figure 3 shows the errors due to RT approximations. As we mentioned above, the v1 forward model performs scalar simulations for all wavelengths, causing errors up to ∼ 10 % compared to vector simulations (Fig. 3a). And then the vector simulations are additionally performed at 14 wavelengths for adjusting the vector vs. scalar differences. However, As shown in Fig. 3b, second-order errors (∼ 0.2 %) remain due to neglecting the dependence of polarization effects on the fine structures of ozone absorption. Using eight streams causes errors of ∼ 0.05 % above 320 nm (Fig. 3c), whereas using 24 layers causes 1 % errors at shorter UV wavelengths (Fig. 3c). Moreover, to improve the v2 simulations we decide to set up 12 streams and 72 layers as well as more wavelengths in the polarization correction.
The right panel of Fig. 1 illustrates the flowchart of the updated forward model scheme (v2) which employs the PCA-based RT model to perform online scalar simulations using four streams and a 24-layer atmosphere for RT performance enhancement (step 1) and two kinds of correction schemes for accounting for approximation errors (steps 2 and 3). Section 3.1.1 gives an overview on how the PCA tool is combined with the VLIDORT version 2.8 model; full theoretical details may be found in Spurr et al. (2016) and Kopparla et al. (2017). Here, our paper gives details on how the PCA-based RT configuration is optimized for the application to UV ozone profile retrievals for maximizing the speed-up in Sect. 3.1.2. Section 3.2 specifies step 2, wherein the LUT-based correction is applied to approximation errors due to the use of a scalar model, a smaller number of streams, and coarser-resolution vertical grid. In step 3 the undersampling correction is adopted from the v1 implementation, but the Rayleigh scattering term of the Eq. (1) is neglected for the speed up with trivial loss of accuracy.
3.1 PCA-based RT model
3.1.1 General PCA procedure
The PCA-based RT process begins with a grouping of spectral points into several bins; atmospheric profile optical properties within each bin are similar. PCA is a mathematical transformation that converts a correlated mean-subtracted dataset into a series of principal components (PCs). To enhance RT performance, PCA is used to compress a binned set of correlated optical profile data into a small set of atmospheric profiles which capture the vast majority of the data variance within the bin. The layer extinction optical thickness Δn,i and the single scattering albedos ωn,i are generally subjected to PCA, where n and i are indices for atmospheric layers () and spectral points (), respectively. For each bin, the optical profiles ln Δn,i and ln ωn,i are composed of 2NL×NS matrix G in log space (, ). The mean-removed 2NL×2NL covariance matrix Y is then
where 〈〉 denotes spectral averaging over all grid points in a bin. This covariance matrix Y is decomposed into eigenvalues ρk and unit eigenvectors Xk through solution of the eigenvalue problem YXk=ρkXk, where the index k runs from 1 to 2NL. The scaled eigenvectors of the covariance matrix are defined as the empirical orthogonal function (EOF), , where the EOFs are ranked in descending order starting with those having the largest eigenvalues. The principal components (PCs) are the projections of the original data onto the eigenvectors, . The original dataset can then be expanded in terms of the mean value and a sum over all EOFs. As inputs to the RT simulation, the PCA-defined optical states are defined as F0=exp [〈G〉] and , corresponding respectively to the mean value and to positive and negative perturbations from the mean value by an amount equal to the magnitude of kth EOF. Therefore, Δn,i and ωn,i (i=1…NS) are expressed as follows:
For those optical quantities not included in the PCA reduction but still required in the RT simulations, the spectral mean values for the bin are assumed, as long as they have smooth monotonic spectral dependency or else are constant over the bin range. In our application, the phase functions and phase matrices for Rayleigh scattering are derived from bin-average values of the depolarization factor. Surface Lambertian albedos are constant in the RT simulation, but the calculated radiance is later adjusted to account for first-order wavelength dependency using surface albedo weighting functions. For larger bins, it is possible to include the depolarization ratio or the Lambertian albedo as additional elements in the optical dataset subject to PCA; this was investigated in another context by Somkuti et al. (2017).
In the PCA-based RT package, three independent RT models are combined in order to generate the full-scattering intensity field (IFULL) at each spectral point λi in a single bin as follows:
Two fast-RT models, the first order (FO) and 2STREAM (2S), are used to generate an accurate single scatter (SS) field (IFO) and an approximate multiple scatter (MS) field (I2S), respectively, for every spectral point. The scalar 2S model computes the radiation field with two discrete ordinates only. To derive the correction factors C(λi), we first compute (logarithmic) ratios of the full-scatter and 2S-based intensity fields calculated with PCA-derived optical states F0 and :
Intensity ratios at the original spectral points J(λi) are then obtained using a second-order central-difference expansion based on the PCA principal components Pki:
The correction factors C(λi)=exp [J(λi)] are then applied to the approximate simulation (I2S(λi)+IFO(λi)) according to Eq. (4) above. More details can be found in the literature (Natraj et al., 2005, 2010; Spurr et al., 2013, 2016; Kopparla et al., 2017).
So far, we have discussed generation of total intensity field, using values IFO(λi) and I2S(λi) from full-spectrum FO and 2S model calculation, as well as PCA-derived values IVLD(F), I2S(F), and IFO(F) based on PCA-derived optical states . The above procedure works with VLIDORT operating in scalar or vector mode; however, the 2S model is purely scalar and cannot be used if we want to establish PCA-RT approximations to the Q and U components of the Stokes vector with polarization present. Instead, we rely on just the VLIDORT and FO models and develop a PCA-RT scheme based on the differences between the values calculated by VLIDORT and FO for monochromatic and PCA-derived calculations, with an additive correction factor instead of the logarithmic ratios in Eq. (6) above. This was first introduced in Natraj et al. (2010) and is discussed in detail in Spurr et al. (2016).
Of greater importance for us is the need to derive PCA-RT approximations to profile Jacobians (weighting functions of the total intensity with respect to ozone profile optical depths). A PCA-RT Jacobians scheme was developed by Spurr et al. (2013) for total column Jacobians in connection with the retrieval of total ozone; this scheme involved formal differentiation of the entire PCA-RT system as outlined above for the intensity field. This is satisfactory for bulk property Jacobians, but for profile Jacobians it is easier to write (Efremenko et al., 2014; Spurr et al., 2016)
Here, , with similar definitions for the FO and VLIDORT partial derivatives with respect to a parameter ξ. The Jacobian correction factor D(ξ)(λi)=exp [L(ξ)(λi)] is determined using the same central-difference expansion as that in Equation (6), but with quantities
in place of J0 and in Eq. (5).
3.1.2 The binning scheme
The major performance saving is achieved by limiting full-MS VLIDORT calculations to those based on the reduced set of PCA-derived optical states F0 and . A general binning scheme has been developed over the shortwave region from 0.29 to 3.0 µm (Kopparla et al., 2016), whereby the entire region is divided into 33 specially chosen sub-windows encompassing the major trace gas absorption signatures; in each such sub-window there are 11 bins for grouping optical properties and up to four EOFs for each PCA bin treatment; with this scheme, radiance accuracies of 0.1 % can be achieved throughout the region. However, the binning scheme should be tuned to the specific application to get additional computational saving, and here we investigate the optimal set for spectral binning and the number of EOFs in the Hartley and Huggins ozone bands (265–360 nm).
Optical properties within each bin must be strongly correlated to reduce the number of EOFs required to attain a given accuracy. According to Kopparla et al. (2016), the UV region is divided at 340 nm, beyond which O2–O2 absorption must be considered. In our application, the spectral region 340–360 nm is further divided at 350 nm: in the first sub-window, ozone absorption is much stronger than O2–O2, while for the second (350–360 nm), O2–O2 absorption becomes dominant. The binning criteria are generally determined by similarities in total optical depth of gas absorption profiles τij as defined below:
where NL and Ng denote the number of atmospheric layers and atmospheric trace gases.
To evaluate the PCA approximation, the “exact-RT” model is performed, where accurate full-MS VLIDORT calculations are expensively performed at every wavelengths in addition to accurate SS calculations:
We first evaluate the impact of applying different binning steps and numbers of EOFs in Fig. 4 where the residuals (IPCA−IEXACT) are plotted as a function of Γg for the spectral window 265–340 nm at small and large SZAs, respectively. In this evaluation, the bins are equally spaced in Γg for the five steps from 0.20 to 1.0. For Γg < 1, where the extinction is strong enough that radiances are very small, the residuals are effectively reduced by having more bins rather than increasing the number of EOFs. In this optical range, using the first EOF is enough to capture the vast majority of the spectral variance, with the optimization of the binning step. However, the bins should be narrowly spaced with Γg intervals of at least 0.3–0.4 for those spectral grids for which Γg is less than −2. These spectral grids are correlated with the Hartley band above ∼ 300 nm, where radiance values rapidly increase due to decreasing ozone absorption, but the spectral variations are almost unstructured. The rest of our spectral region corresponds to the Huggins band above 310 nm, where spectral variations are distinctly influenced by local maxima and minima of ozone absorption. In this spectral region, PCA approximation errors can be greatly reduced by increasing the number of EOFs. However, it is interesting to note that the PCA approximation is not further improved by using four EOFs instead of three (not shown here). Figure 4 also illustrates the dependence of the PCA performance on SZA in the spectral range below 340 nm: for example, when two EOFs are applied with the binning step 0.4, errors are within ±0.02 % at smaller SZA but increase up to ±0.03 % at larger SZA. Therefore, as listed in Table 1, two sets of binning criteria are determined to keep the accuracy within 0.05 % for any viewing geometry. Based on the experiments shown in Fig. 5, the binning criteria are determined for the other sub-windows listed in Table 1, namely 340–350 and 350–360 nm: the former is set with bins at intervals of 1 and using the first two EOFs, while the latter is divided into a single bin with the first four EOFs. Figure 6 illustrates the binning criteria thus determined, demonstrating that the PCA performance keeps accuracies within 0.03 % when various sets of SZAs, ozone profiles, and vertical layers are implemented.
3.2 LUT-based correction
Two sets of LUTs are created: for high-accuracy (LUTH: vector, 12 streams, 72 layers) and low-accuracy (LUTL: scalar, 4 streams, 24 layers) configurations. The online PCA-VLIDORT model is configured to run in the LUTL mode. The correction spectrum is straightforwardly calculated as the ratio of the LUT-based spectrum (LUTH LUTL), but the radiance correction term is additionally adjusted to account for the different gas optical depth profiles used in online and LUT simulations. The RT results are corrected for each wavelength as follows:
where the subscripts “on” and “LUT” stand for online and LUT-based calculations, respectively; As and τn represent the surface albedos and gas absorption optical depths (n is the layer index). To construct LUTs, RT calculations are performed using the VLIDORT version 2.8 model for sets of geometrical configurations (θo,θ; solar zenith angle, viewing zenith angle), surface pressures for 22 climatological ozone profiles, and 92 wavelengths (265–345 nm) as listed in Table 2. The azimuth dependence is treated exactly using the 0–2 Fourier intensity components in a Rayleigh scattering atmosphere in conjunction with the associated cosine-azimuth expansion of the full intensity; see the discussion below. The 22 ozone profiles are constructed from the GOME ozone profile product (Liu et al., 2005), where the ozone profile shapes vary according to three latitude regimes and with the total column ozone amounts at 50 DU intervals. The 92 wavelengths are regularly sampled at 5 nm intervals below 295 nm and at 1.0 nm intervals up to 310 nm in the Hartley band, while they are irregularly sampled at the local minima and maxima of the ozone absorption structures in the Huggins band. The results of these RT calculations are separated into two components: the path radiance Iatm and the surface reflectance term Isfc according to Chandrasekhar (1960), so that the following relationship is employed to recover the full radiance:
Iatm represents the purely atmospheric contribution to the radiance in the presence of a dark surface (zero albedo), and in a Rayleigh scattering atmosphere this is given as a Fourier expansion in the cosine of the relative azimuth angle.
However, it is more convenient to write this in the form
In the LUTs, the three coefficients (I0, Z1, and Z2) are stored instead of Iatm. Note that the use of terms aq1 and aq2 is taken from Dave (1964); most of the angular variability in components I1 and I2 is captured analytically with these functions. In other words, Z1 and Z2 are angularly smooth and well-behaved (non-singular) functions, which helps improve angular interpolation accuracy with fewer points in the angular grids. The surface term is
In the LUTs, we store the transmission term T(θo,θ), which is the product of the atmosphere downwelling flux transmittance for a solar source with the upwelling transmittance from a surface illuminated isotropically from below and the geometry-independent term s∗ which is the spherical albedo from such a surface. This is the so-called “planetary problem” calculation (Chandrasekhar, 1960), and the code to obtain T and s∗ is now implemented in VLIDORT version 2.8 (Spurr and Christi, 2019). One of the key features of the VLIDORT code is its ability to generate simultaneously (along with the Stokes vector radiation field) any set of Jacobians with respect to atmospheric and surface optical properties. VLIDORT also contains an analytical linearization of the planetary problem. Indeed, in our Rayleigh-based application, we require Jacobians with respect to the albedo As and the ozone profile elements τ. First, for the albedo weighting function we have straightforward differentiation from Eq. (15) as follows:
For the optical depth derivative, is calculated from
All partial derivatives in this expression are returned automatically by VLIDORT. For a given ozone profile, wavelength, and surface pressure, the number of the LUT values specified in Table 3 is 770 (, nVar=8: ), which is much smaller than that of a LUT with dependence on eight relative azimuth angles and five surface albedo values (11 520 , nVar=3: ). LUT-based simulated radiances are evaluated against online simulations: the LUT interpolation errors are mostly less than 0.2 %–0.3 % (not shown here), except for extreme path length scenarios (e.g., ∼ 1 % at θo=87.0∘) as shown in Fig. 7a, b; however, the interpolation errors are quite similar to each other for LUTH and LUTL. Therefore, those errors are canceled out when performing corrections using these two LUTs, and thereby the overall error after LUT correction is much smaller than ∼ 0.05 % (Fig. 7c). Note that the accuracy is completely maintained with respect to both φ−φo and As, while the size of a LUT is reduced by a factor of 15. However, LUT corrections still contain ozone profile shape errors due to the use of 22 representative total ozone-dependent ozone profiles in the LUT. Figure 8 shows an example of the correction spectrum as a function of SZA, showing that polarization errors are mostly dominant, except at the high SZAs above 310 nm, where errors due to use of a low number of streams become significant, and for wavelengths below 300 nm, where the use of the coarse vertical layering scheme becomes the main source of uncertainty.
∗ Total ozone-based ozone profiles for three latitude regimes. The grid values represent the amount of total ozone (DU).
a Fourier coefficients of path radiance with respect to relative azimuth angle (AZA). b Total transmission of the atmosphere. c Spherical albedo of the atmosphere. d Total gas absorption optical depth profile. e n: number of atmospheric layers.
The PCA-RT model developed as described in this paper is implemented as the forward model component of an iterative optimal-estimation-based inversion (Rodgers, 2000) for retrieving the ozone profile from OMI measurements. In previous studies, the PCA-RT performance was evaluated against a suite of exact monochromatic baselines of fully accurate VLIDORT simulations. However, such exact RT calculations cannot be applied in the operational data processing system, especially when thousands of spectral points are involved; in other words, the operational capability of the PCA-RT approach has been overestimated in previous studies. Therefore, we evaluate the RT model developed against the existing forward model where many RT approximations are applied to meet the computational budget in the operational system.
a Undersampled (US) spectral intervals (nm) used to define wavelengths at which RT is actually executed. “0.3|0.1” represents the intervals divided at 305 nm, while for “1.0|0.4|0.6” those are divided at 295 and 310 nm. b High-resolution (HR) spectral intervals (nm) used to define wavelengths where undersampled simulations are interpolated before spectral convolution. c The number of discrete ordinates in the full polar space. d The number of atmospheric layers. e RT model is run in the vector (scalar) mode if polarization is true (false). f Online correction is performed for polarization errors, based on Liu et al. (2010). LUT-based correction is performed for RT approximation errors due to neglecting polarization as well as using four streams and 24 layers, developed in this study.
Table 4 contains sets of configurations for seven forward models. OMI spectra are simulated at the undersampled (US) intervals specified in the first column of this table and then interpolated at high-resolution (HR) intervals (second column) with the undersampling correction before convolution with OMI slit functions. In the v1 forward model, the US spectral intervals were set at 1.0 nm/0.4 nm intervals below/above 295 and 0.6 nm above 310 nm, while the HR spectral interval was set at 0.05 nm. In the updated RT model, the spectral points are selected at 0.3 nm (0.1 nm) intervals below (above) 305 nm, and the HR interval is set as 0.03 nm, which enables us to achieve very high accuracy, better than 0.01 %, as shown in Fig. 2c. In the reference configuration (abbreviated to “Ref”), VLIDORT is run in vector mode with 12 streams and 72 atmospheric layers so that the RT approximation errors are significantly reduced. The VLIDORT-based forward model is run with five sets of configurations (abbreviated to VLD in Table 4) to quantify the impact of RT approximations on ozone retrievals. Figure 9 compares the mean biases of the retrieved ozone profiles between VLD/PCA and Ref for three SZA regimes. VLD0 represents the v1 forward model configuration, demonstrating that the ozone retrieval errors due to the entire forward model errors range from ∼ 3.5 % for the large SZA regime to ∼ 5.5 % for the small SZA regime at the lower atmospheric layers but ∼ 2 % at the upper layers. The configuration VLD1 assesses the impact of undersampling errors on the retrievals, causing negative biases of up to 2.0 % below ∼ 20 km. Compared to the use of 12 streams, using 8 streams causes negligible impacts on ozone retrievals (VLD2) as the corresponding RT model approximation errors are negligible, except for extreme viewing geometries where the ozone retrieval errors are overwhelmed by instrumental measurement errors (a few %) rather than the forward model errors of ∼ 0.05 % as shown in Fig. 3c. The VLD3-based RT calculation is applied to ozone retrievals for evaluating online polarization correction, showing that the corresponding errors in tropospheric ozone retrievals are estimated as ±2 % at small SZAs. The evaluation for VLD4 demonstrates that the use of coarse atmospheric layering causes the largest errors (∼ 4.5 % in the troposphere, ∼ 1.5 % in the stratosphere). PCA0 represents the v2 forward model configuration while PCA1 is done with the highest accurate configuration except for PCA approximation. Retrieval errors due to PCA approximation are negligible except for the bottom few layers at smaller solar zenith angles (up to ∼ 1.5 %). Differences between PCA0 and PCA1 represent the ozone retrieval errors due to LUT errors, mostly related to the profile shape errors between LUT and online calculations. In Fig. 10, the comparison between VLD (v1 PROFOZ) and PCA (v2 PROFOZ) is performed for individual ozone profile retrievals. The large systematic errors of ∼ 5 %–15 % due to v1 forward model errors are greatly eliminated below 30 km. In addition, the variabilities of individual differences are significantly eliminated over the entire layers at high solar zenith angles. However, there are still some remaining retrieval errors up to −5 % in the troposphere and 3 % in the stratosphere due to v2 forward model simulation errors. Figure 11 further evaluates the v2 implementation. First of all, the comparison of the runtime (Fig. 11a) demonstrates that v2 is faster by a factor of 3.3 on average. Some spectral fit residuals are eliminated in the UV 1 band over the middle area of the swath (low latitudes), where the SZAs are relatively small, by up to ∼ 2 %; the corresponding improvements are found in the stratospheric column ozone. The amount that the stratospheric column ozone deviated from the reference is reduced by ∼ 0.2 % with the v2 implementation. On the other hand, the tropospheric column ozone retrievals show improvements for most cases, whereas the fit residuals of the UV2 band are slightly worse in the low latitudes. Note that the smaller fitting residuals could not directly lead to better ozone retrievals likely due to the presence of systematic measurement errors.
We have extended the PCA-based fast-RT method to improve computational challenges for OE-based SAO OMI ozone profile retrievals requiring iterative calculations of the radiance and its Jacobian derivatives. The PCA-RT model is designed to perform MS calculations for a few EOF-derived optical states which are developed from spectrally binned sets of inherent optical properties that possess some redundancy. In this study, the binning scheme is carefully turned for the UV ozone fitting window from 265 to 360 nm in such a way as to choose the number of EOFs to be as small as possible for each bin rather than always using the first four EOFs for all bins selected in previous studies. The spectral windows are divided into three sub-windows: (1) 265–340 nm, (2) 340–350 nm, and (3) 350–360 nm. Then, optical profiles are grouped into bins according to criteria based on the total gas optical depth, as specified in Table 1. We demonstrated that the PCA approximation errors for our application are within 0.03 % for any viewing geometry, optical depth profile, and vertical layering.
The existing (v1) forward model calculations are evaluated to determine the optimal configuration for the v2 forward model. RT approximation errors exist due to the use of 24 quite coarse vertical layers (2.5 km thick), which can cause radiance simulation errors of up to ∼ 1 % below 320 nm, and this leads to ozone retrieval errors of 2 %–4.5 % in the troposphere and 1.5 % in the stratosphere. Eight-stream calculations can result in radiance residuals of ∼ 0.05 % or less except at extreme viewing geometries, which causes trivial errors on ozone retrievals compared to other error factors. In spite of accounting for polarization errors using vector and scalar differences at 14 wavelengths, the retrieval accuracies are systematically worse by up to ∼ 2 % due to neglecting second-order polarization errors which are strongly correlated with ozone absorption features. We found that 72 atmospheric layers (∼ 0.7 km thick) and 12 streams should be used at least to improve the simulations comparable to those with 99 atmospheric layers and 32 streams. To reduce the impact of undersampling errors, we improve simulation intervals such as 0.3 nm below 305 nm and 0.1 nm above and thereby reduce the biases of the ozone retrievals by ∼ 1.5 % compared to the undersampled intervals used in the v1 simulation. Applying the PCA-RT approach allows us to reduce the number of MS calculations from the high-resolution optical dataset to 51 sets of EOF-derived optical states, but the performance savings are not enough to improve over previous RT approximations. To improve both efficiency and accuracy, we have developed a LUT-based correction for eliminating the RT approximation errors arising from the vector vs. scalar, 12 vs. 4 streams, and 72 vs. 24 layers. In conclusion, the updated PCA-based RT model combined with LUT corrections makes ozone profile retrievals faster than the v1 forward model by a factor of 3.3 on average. Fitting accuracies are significantly improved in the UV1 band by 2 % and comparable in the UV2 band, while the ozone profile retrievals are significantly improved, especially in the troposphere, by ∼ up to 10 %. However, there are still some remaining retrieval errors of up to −5 % in troposphere and 3 % in stratosphere due to the LUT correction errors and PCA approximation errors in the v2 implementation. The updated forward model is in preparation for reprocessing all OMI measurements (2004–current) for the next version of the PROFOZ product.
OMI Level 1b radiance datasets are available at https://doi.org/10.5067/Aura/OMI/DATA1002 (Dobber, 2007) (last access: 11 March 2021). The LUT database is attainable upon request.
JB and XL designed the research; RS provided oversight and guidance for using both VLIDORT and PCA-based VLIDORT; KY developed the LUT creation and interpolation scheme; XL contributed to analyzing ozone profile retrievals with different forward model approaches; JB conducted the research and wrote the paper; CRN, CCM, GGA, and KC contributed to the analysis and writing; CCM and GGA contributed to managing the computational resources.
The authors declare that they have no conflict of interest.
We acknowledge the OMI science team for providing their satellite data. Research at the Smithsonian Astrophysical Observatory is funded by the NASA Aura science team program (NNX17AI82G). Research at Pusan National University is funded by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education(2020R1A6A1A03044834). Both calculations and simulations are done on the Smithsonian Institution High-Performance Cluster (SI/HPC) computer system (https://doi.org/10.25572/SIHPC).
This research has been supported by the NASA Aura science team program (grant no. NNX17AI82G) and the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (grant no. 2020R1A6A1A03044834).
This paper was edited by Cheng Liu and reviewed by two anonymous referees.
Bak, J., Liu, X., Sun, K., Chance, K., and Kim, J.-H.: Linearization of the effect of slit function changes for improving Ozone Monitoring Instrument ozone profile retrievals, Atmos. Meas. Tech., 12, 3777–3788, https://doi.org/10.5194/amt-12-3777-2019, 2019.
Brion, J., Chakir, A., Daumont, D., and Malicet, J.: High-resolution laboratory absorption cross section of O3, Temperature effect, Chem. Phys. Lett., 213, 610–612, 1993.
Chandrasekhar, S.: Radiative Transfer, Dover Publications, Mineola, New York, USA, 1960.
Cuesta, J., Eremenko, M., Liu, X., Dufour, G., Cai, Z., Höpfner, M., von Clarmann, T., Sellitto, P., Foret, G., Gaubert, B., Beekmann, M., Orphal, J., Chance, K., Spurr, R. J. D., and Flaud, J.-M.: Satellite observation of lowermost tropospheric ozone by multispectral synergism of IASI thermal infrared and GOME-2 ultraviolet measurements over Europe, Atmos. Chem. Phys., 13, 9675–9693, https://doi.org/10.5194/acp-13-9675-2013, 2013.
Dave, J. V.: Meaning of Successive Iteration of the Auxiliary Equation in the Theory of Radiative Transfer, Astrophys. J., 140, 1292–1303, 1964.
Dobber, M.: OMI/Aura Level 1B UV Global Geolocated Earthshine Radiances 1-orbit L2 Swath 13 × 24 km V003, Greenbelt, MD, USA, Goddard Earth Sciences Data and Information Services Center (GES DISC), https://doi.org/10.5067/Aura/OMI/DATA1002, 2007.
Efremenko, D. S., Loyola, D. G., Spurr, R. J. D., and Doicu, A.: Acceleration of Radiative Transfer Model Calculations for the Retrieval of Trace Gases under Cloudy Conditions, J. Quant. Spectrosc. Ra., 135, 58–65, 2014.
Fu, D., Worden, J. R., Liu, X., Kulawik, S. S., Bowman, K. W., and Natraj, V.: Characterization of ozone profiles derived from Aura TES and OMI radiances, Atmos. Chem. Phys., 13, 3445–3462, https://doi.org/10.5194/acp-13-3445-2013, 2013.
Fu, D., Kulawik, S. S., Miyazaki, K., Bowman, K. W., Worden, J. R., Eldering, A., Livesey, N. J., Teixeira, J., Irion, F. W., Herman, R. L., Osterman, G. B., Liu, X., Levelt, P. F., Thompson, A. M., and Luo, M.: Retrievals of tropospheric ozone profiles from the synergism of AIRS and OMI: methodology and validation, Atmos. Meas. Tech., 11, 5587–5605, https://doi.org/10.5194/amt-11-5587-2018, 2018.
Kim, P. S., Jacob, D. J., Liu, X., Warner, J. X., Yang, K., Chance, K., Thouret, V., and Nedelec, P.: Global ozone–CO correlations from OMI and AIRS: constraints on tropospheric ozone sources, Atmos. Chem. Phys., 13, 9321–9335, https://doi.org/10.5194/acp-13-9321-2013, 2013.
Kopparla, P., Natraj, V., Spurr, R. J. D., Shia, R. L., Crisp, D., and Yung, Y. L.: A fast and accurate PCA based radiative transfer model: Extension to the broadband shortwave region, J. Quant. Spectrosc. Ra., 173, 65–71, 2016.
Kopparla, P., Natraj, V., Limpasuvan, D., Spurr, R. J. D., Crisp, D., Shia, R. L., Somkuti, P., and Yung, Y. L.: PCA-based radiative transfer: Improvements to aerosol scheme, vertical layering and spectral binning, J. Quant. Spectrosc. Ra., 198, 104–111, 2017.
Kroon, M., de Haan, J. F., Veefkind, J. P., Froidevaux, L., Wang, R., Kivi, R., and Hakkarainen, J. J.: Validation of operational ozone profiles from the Ozone Monitoring Instrument, J. Geophys. Res., 116, D18305, https://doi.org/10.1029/2010JD015100, 2011.
Liu, C., Yao, B., Natraj, V., Kopparla, P., Weng, F., Le, T., Shia, R., and Yung, Y. L.: A spectral Data Compression (SDCOMP) Radiative Transfer Model for High-Spectral-Resolution Radiation Simulations, J. Atmos. Sci., 77, 2055–2066, https://doi.org/10.1175/JAS-D-19-0238.1, 2020.
Liu, X., Chance, K., Sioris, C. E., Spurr, R. J. D., Kurosu, T. P., Martin, R. V., and Newchurch, M. J.: Ozone profile and tropospheric ozone retrievals from Global Ozone Monitoring Experiment: algorithm description and validation, J. Geophys. Res., 110, D20307, https://doi.org/10.1029/2005JD006240, 2005.
Liu, X., Bhartia, P. K., Chance, K., Spurr, R. J. D., and Kurosu, T. P.: Ozone profile retrievals from the Ozone Monitoring Instrument, Atmos. Chem. Phys., 10, 2521–2537, https://doi.org/10.5194/acp-10-2521-2010, 2010.
Liu, X., Yang, Q., Li, H., Jin, Z., Wu, W., Kizer, S., Zhou, D. K., and Yang, P.; Development of a fast and accurate pcrtm radiative transfer model in the solar spectral region, Appl. Optics, 55, 8236–8247, 2016.
Loyola, D. G., Gimeno García, S., Lutz, R., Argyrouli, A., Romahn, F., Spurr, R. J. D., Pedergnana, M., Doicu, A., Molina García, V., and Schüssler, O.: The operational cloud retrieval algorithms from TROPOMI on board Sentinel-5 Precursor, Atmos. Meas. Tech., 11, 409–427, https://doi.org/10.5194/amt-11-409-2018, 2018.
Miles, G. M., Siddans, R., Kerridge, B. J., Latter, B. G., and Richards, N. A. D.: Tropospheric ozone and ozone profiles retrieved from GOME-2 and their validation, Atmos. Meas. Tech., 8, 385–398, https://doi.org/10.5194/amt-8-385-2015, 2015.
Nanda, S., de Graaf, M., Veefkind, J. P., ter Linden, M., Sneep, M., de Haan, J., and Levelt, P. F.: A neural network radiative transfer model approach applied to the Tropospheric Monitoring Instrument aerosol height algorithm, Atmos. Meas. Tech., 12, 6619–6634, https://doi.org/10.5194/amt-12-6619-2019, 2019.
Natraj, V., Jiang, X., Shia, R.-L., Huang, X., Margolis, J. S., and Yung, Y. L.: Application of principal component analysis to high spectral resolution radiative transfer: A case study of the O2−A band, J. Quant. Spectrosc. Ra., 95, 539–556, https://doi.org/10.1016/j.jqsrt.2004.12.024, 2005.
Natraj, V., Shia, R. L., and Yung, Y. L.;. On the use of principal component analysis to speed up radiative transfer calculations, J. Quant. Spectrosc. Ra., 111, 810–816, https://doi.org/10.1016/j.jqsrt.2009.11.004, 2010.
Rodgers, C. D.: Inverse methods for atmospheric sounding – Theory and practice, in: Series on atmospheric oceanic and planetary physics, vol. 2, edited by: Rodgers, C. D., World Scientific Publishing Co. Pte. Ltd., Oxford, UK, 256, https://doi.org/10.1142/9789812813718, 2000.
Somkuti, P., Boesch, H., Vijay, N., and Kopparla, P.: Application of a PCA-based fast radiative transfer model to XCO2 retrievals in the shortwave infrared, J. Geophys. Res.-Atmos., 122, 10477–10496, https://doi.org/10.1002/2017JD027013, 2017.
Spurr, R. J. D.: Linearized pseudo-spherical scalar and vector discrete ordinate radiative transfer models for use in remote sensing retrieval problems, in: Light Scattering Reviews, edited by: Kokhanovsky, A., Springer, New York, USA, 2008.
Spurr, R. J. D. and Christi, M.: The LIDORT and VLIDORT Linearized Scalar and Vector Discrete Ordinate Radiative Transfer Models: An Update for the last 10 years, in: Light Scattering Reviews, vol. 12, edited by: Kokhanovsky, A., Springer, New York, USA, 2019.
Spurr, R. J. D., Natraj, V., Lerot, C., Van Roozendael, M., and Loyola, D.: Linearization of the principal component analysis method for radiative transfer acceleration: Application to retrieval algorithms and sensitivity studies, J. Quant. Spectrosc. Ra., 125, 1–17, https://doi.org/10.1016/j.jqsrt.2013.04.002, 2013.
Spurr, R. J. D., Natraj, V., Kopparla, P., and Christi, M.: Application of Principal Component Analysis (PCA) to Performance Enhancement of Hyperspectral Radiative Transfer Computations, in: Principal Component Analysis: Methods, Applications and Technology, chapter 3, Nova Science Publishers, Hauppauge, NY, USA, available at: https://novapublishers.com/shop/principal-component-analysis-methods-applications-and-technology/ (last access: 11 March 2021, ISBN 978-1-53610-889-7, 2016.
Yang, Q., Liu, X., Wu, W., Kizer, S., and Baize, R. R.: Fast and accurate hybrid stream PCRTM-SOLAR radiative transfer model for reflected solar spectrum simulation in the cloudy atmosphere, Opt. Express, 24, A1514–A1527, 2016.
Zoogman, P., Liu, X., Suleiman, R. M., Pennington, W. F., Flittner, D. E., Al-Saadi, J. A., Hilton, B. B., Nicks, D. K., Newchurch, M. J., Carr, J. L., Janz, S. J., Andraschko, M. R., Arola, A., Baker, B. D., Canova, B. P., Chan Miller, C., Cohen, R. C., Davis, J. E., Dussault, M. E., Edwards, D. P., Fishman, J., Ghulam, A., González Abad, G., Grutter, M., Herman, J. R., Houck, J., Jacob, D. J., Joiner, J., Kerridge, B. J., Kim, J., Krotkov, N. A., Lamsal, L., Li, C., Lindfors, A., Martin, R. V., McElroy, C. T., McLinden, C., Natraj, V., Neil, D. O., Nowlan, C. R., O'Sullivan, E. J., Palmer, P. I., Pierce, R. B., Pippin, M. R., Saiz-Lopez, A., Spurr, R. J. D., Szykman, J. J., Torres, O., Veefkindz, J. P., Veihelmann, B., Wang, H., Wang, J., and Chance, K.: Tropospheric Emissions: Monitoring of Pollution (TEMPO), J. Quant. Spectrosc. Ra., 186, 17–39, https://doi.org/10.1016/j.jqsrt.2016.05.008, 2017.