Interactive comment on “ Single Footprint Retrievals for AIRS using a Fast TwoSlab Cloud-Representation Model and All-Sky Radiative Transfer Algorithm

The manuscript ÂńSingle Footprint Retrievals for AIRS using a Fast TwoSlab CloudRepresentation Model and All Sky Radiative Transfer AlgorithmÂż of DeSouzaMachado et al., submitted to Atmopsheric Measurement Techniques presents a study in two parts. The first one compares two infrared cloudy radiative transfer models (RTM), a fast and complex cloudy RTM and a fast and simplified cloudy RTM. Based on their results they used in the second part the fast and simplified cloudy RTM model to improve water vapor and temperature profile retrievals from AIRS cloudy observations. The results presented here are very encouraging and comparison with other


Introduction
Since the early 2000's, a number of high spectral resolution, low noise, very stable new generation hyperspectral infrared (IR) sounders have been deployed on board Earth orbiting satellites, providing daily global Top of the Atmosphere (TOA) radiance spectra.In principle these TOA radiances can be inverted to estimate atmospheric temperature and humidity profiles, minor gas concentration, surface temperature, and some clouds parameters.
IR sounders have rather large nadir footprints of ~15 km diameter, consequently far less than 10% of scenes are cloud free.Retrievals using eigenvalue regression methods have been used with these all-sky (cloud and clear) radiances (see for example Weisz et al. (2007)) but these methods have no reliable error estimates for individual scenes.Existing NASA and NOAA operational retrieval systems for IR sounders use cloud-clearing (Susskind et al., 1998(Susskind et al., , 2003) ) and (Gambacorta, 2013).In addition, Numerical Weather Prediction (NWP) centers generally only assimilate IR sounder radiances that have been deemed clear.
Presently the NASA-AIRS operational soundings are performed using cloud-cleared radiances coupled with a clear-sky RTA (Susskind et al., 1998(Susskind et al., , 2003)).Cloud cleared radiances (CCRs) are synthesized using the differences in cloud amounts in a (typically) 3-by-3 set of adjacent Fields of View (FOVs) to produce a single effective estimate of the clear-sky radiance.This process increases the retrieval yield (to well above 10%) and provides some error estimates but simultaneously reduces the spatial resolution by a factor of three.Publicly available products from 1D variational retrievals include: 1. Atmospheric Infrared Sounder (AIRS): NASA, using cloud-clearing from 3x3 set of footprints (Susskind et al., 1998(Susskind et al., , 2003) ) 2. Cross Track Infrared Sounder (CrIS): The NOAA Unique Combined Atmospheric Processing System (NuCaps), also using cloud-clearing from 3x3 set of footprints (Gambacorta, 2013) 3. Infrared Atmospheric Sounding Interferometer (IASI): NOAA NuCaps, using cloud-clearing from 2x2 set of footprints (Gambacorta, 2013); 4. Infrared Atmospheric Sounding Interferometer (IASI): EUMETSAT, 2-step single-footprint retrievals: piecewise regression for all scenes nominally exploiting IASI in synergy with AMSU+MHS (IASI-only is fallback) followed by a physical retrieval using the Optimal Estimation Method (OEM) (Rodgers, 2000;Steck, 2001) on clear-scenes only (IASI-only) (August et al., 2012;EUMETSAT, 2016) The retrieval approaches mentioned above use various combinations of training to NWP forecasts from the European Center for Medium Range Weather Forecasting (ECMWF) either by regression (IASI EUMETSAT) or with neural nets (AIRS NASA) or use climatology (CrIS NOAA).All utilize co-located microwave soundings when possible.The development of a formal error estimate computation in the NuCaps algorithm is underway (personal communication, Gambacorta).The CCR approaches lead to complicated quality control issues, since cloud-clearing can fail, and the decisions made in assigning quality flags to the retrievals are not trivial.The cloud-clearing process is especially problematic (Zhou et al., 2005) when the cloud fields are homogeneous and cloud-clearing becomes unstable and inaccurate, which introduce errors into retrieved products.This is not necessarily a problem for weather-forecasting oriented applications, since the retrieval Quality Assurance (QA) can accurately determine when cloud-clearing failed.The extensive QA in the AIRS retrieval system deems as many as ~20% of observations The PCRTM MRO implementation (Chen et al., 2013) uses the full vertical cloud profiles in the ECMWF model data.When a 50 sub-column MRO is added to the RTA to represent the cloud sub-grid variability, the radiance computation slows by 10X as compared to a 5 sub-column MRO.The appendix show that hyperspectral infrared radiances typically contain ~2-4 degrees of freedom of cloud information, which could be parametrized by the cloud amount, fraction, and cloud top and bottom pressures.Our approach exploits this to reduce the cloud representation complexity from N-level model cloud fields for cloud ice water content and cloud liquid water content (and cloud cover) into two randomly overlapping slabs within the radiative transfer layers, greatly reducing the computational burden.The speed of the scattering calculations are then comparable to those under clear-sky conditions, and we show below that the radiances are as accurate as those from the MRO scheme.
The SARTA Two-Slab approach is then applied to single-footprint retrievals for an AIRS granule and compared to the existing NASA AIRS Level 2 retrievals.As noted above, a key issue is the proper initialization of the cloud parameters in our RTA.Model fields from ECMWF are used here to initialize the thermodynamic and scattering cloud fields.Although NWP models do a reasonably good job at estimating cloud parameters, it is very unlikely that the positions of the model clouds are correct at scales near the sounder spatial resolution, especially given the time mis-match between available forecast models and the observations (±1.5 hours).Hence, our cloud parameters are chosen using the closest matches between simulated and observed window region radiances, restricting the choices to model grid points close to the observation.This approach is key to the success of these single-footprint retrievals.
There are recent papers detailing hyperspectral Optimal Estimation based retrievals in the presence of clouds, see for example (Wu et al., 2017;Irion et al., 2017).Our approach is slightly different as it uses easily available NWP fields for initialization, and a simple cloud representation which allows for well defined jacobians to retrieve thermodynamic profiles and two cloud decks, leading to high yields.
The paper is organized as follows.The AIRS instrument and the use of the ECMWF model are summarized first, followed by a detailed description of the RTA models and the cloud representation schemes.We then examine the computed radiance differences for both clear-sky and all-sky for these two RTAs and discuss radiance differences arising from perturbations to the TwoSlab cloud representation schemes.Finally we outline a method to reduce the impact of the spatial/ temporal mismatch of observed versus modeled clouds, and use this together with the TwoSlab cloud representation to perform single footprint (cloudy) scene retrievals with an a-priori from the NWP model fields.

The AIRS instrument and data
The Atmospheric Infrared Sounder (AIRS) on board NASA's polar orbiting EOS/Aqua platform has 2378 channels, covering the Thermal Infrared (TIR) spectral range (roughly 649-1613 cm -1 ) and shortwave infrared (2181-2665 cm -1 ).The full widths at half maximum satisfy ν/δν ∼ 1200.The (spectral dependent) noise is typically ≤ 0.2K.The instrument, operational since September 2002 is expected to continue operating until the early 2020's.In the comparisons presented later, we use about 1500 stable channels, avoiding channels that have deteriorated over time.AIRS has a 13.5 km nadir footprint from a ~705 km orbit, and scans about ±49.5 degrees from nadir.Radiances from AIRS have been shown to be very stable and accurate (Aumann et al., 2006).

The ECMWF model fields
The core ECMWF 0-10 day forecasts are produced using the Integrated Forecasting System (IFS) (Uppala et al., 2005;Dee et al., 2011).The topmost ECMWF level is 0.01 mb, with terrain following σ-levels from the surface to 0.01 mb.The level spacing is finest in the boundary layer and coarsest at the top.At each model grid point, the cloud fields include the Cloud Water/Ice Water Content profiles (CIWC(z),CLWC(z) in g/g), Cloud Cover profile (CC(z)) and Total Cloud Cover (tcc).
Here we use 91 level ECMWF model fields, at a horizontal resolution of 0.25 • (about 14 km at the equator, approximately the same size as the 13.5 km AIRS nadir footprint).AIRS is on a 1:30 pm equator ascending overpass orbit while ECMWF analysis/forecast are output at 3 hour intervals (8 per day) starting at 00.00 GMT.The closest in time forecast/analysis output is used to provide gridded fields, which are matched (using nearest grid point) to the AIRS L1B observations.This means the profile versus observed (latitude,longitudes) match ups are within 0.25 ±0.05 of each other, while the time differences are uniformly distributed within ±1.5 hours.
The topmost AIRS RTA pressure layer boundary is 0.005 mb, so US Standard temperature, water vapor and ozone fields (McClatchey et al., 1972) are appended above the 0.01 mb boundary.Standard profiles are also used for the remaining atmospheric gases, with carbon dioxide (CO 2 ) and methane (CH 4 ) concentrations set to 385 ppmv and 1.8 ppmv at the surface.
We use two different RTAs, described below, to simulate AIRS infrared radiances that differ primarily in the scattering radiative transfer.Both RTAs use the same AIRS 100 pressure layer scheme (Strow et al., 2003); layer thicknesses range from 0.25 km at the surface, 0.75 km at the upper tropopause, and about 4 km at 0.005 mb (about 80 km) which is the TOA for the model.

SARTA
The clear sky version (with gray cloud capability) of SARTA is used for the NASA AIRS Level 2 retrievals.Layer optical depths are generated using pre-computed predictor coefficients (Aumann and Pagano, 2002;Strow et al., 2003Strow et al., , 2006)).SARTA is trained using optical depths from the pseudo line-by-line (LBL) kCompressed Atmospheric Radiative Transfer Algorithm (kCARTA) package (De Souza-Machado et al., 2002).SARTA has been validated during dedicated AIRS campaigns (Strow et al., 2006).(Chou et al., 1999) algorithm.The PCLSAM algorithm recasts the extinction, single scattering albedo and asymmetry factor due to clouds and aerosols, into an effective absorption optical depth, and is used in other infrared fast models, see for example Matricardi (2005); Vidot et al. (2015); Liuzzi et al. (2016).For each SARTA AIRS layer that contains a cloud/aerosol, the total optical depth is then the sum of the atmospheric gas optical depth plus the cloud/aerosol effective optical depth.Fast, efficient clear sky radiative transfer can then be used to compute the TOA radiance, and to compute finite difference jacobians.Cirrus cloud scattering parameters come from Baum et al. (2007Baum et al. ( , 2011)), while water cloud scattering parameters come from Mie coefficients integrated over a gamma distribution of variance 0.1.

PCRTM
We bench-marked the SARTA TwoSlab model versus the radiance simulator based on the Principal Component Radiative Transfer Model (PCRTM) (Liu et al., 2006(Liu et al., , 2009) ) with full accounting of the cloud sub-grid variability (Chen et al., 2013).
PCRTM is a fast model that computes atmospheric optical depths based on the Line-by-Line Radiative Transfer Model (LBLRTM) (Liu et al., 2006).PCRTM uses an adding cloudy sky scheme (Liu et al., 2009), based on reflectance and transmittance trained using 32-stream Discrete Ordinates Radiative Transfer Program for a Multi-Layered Plane-Parallel Medium (DISORT) (Stamnes et al., 1988) to simulate the effects of ice and water clouds.Again, ice scattering coefficients come from Baum et al. (2011) while the refractive indices for water come from Segelstein (1981).Unlike conventional channel-based radiative transfer models which compute the radiance of each channel separately, the PCRTM calculates the scores (i.e. the coefficients) of pre-computed principal components (PCs) in the spectral domain, with the instrument spectral response function taken into account.The PC scores contain essential information about the radiances, and can be calculated by performing monochromatic radiative transfer calculations at a small number of frequencies.The spectral radiances are then computed by multiplying the PC scores with pre-computed PCs.With this approach, the PCRTM achieves both high accuracy and extremely fast computational and high storage efficiency (Liu et al., 2006;Chen et al., 2013).Chen et al. (2013) showed that the rootmean-square differences between the PCRTM and LBLRTM (a widely used line-by-line radiative transfer benchmark model (Clough et al., 2005)) are 0.67 K for the clear-sky case and 0.78 K for the overcast case, for a wavenumber range spanning 0-2000 cm -1 .Chen et al. (2013) implemented a radiance simulator using the PCRTM and taking cloud variability into account in the same way the ISSCP simulator does (Klein and Jakob, 1999).

Cloud Model Field Conversion
Here we describe the TwoSlab cloud representation and the MRO cloud models.The latter is used exclusively with PCRTM, and the former with SARTA except in Section 4.2 when both PCRTM and SARTA use the TwoSlab model for inter-RTA comparison purposes.The MRO model has been previously documented and is briefly summarized at the end of this section.

TwoSlab conversion
Our cloud representation scheme replaces the N-level NWP cloud vertical profiles by one or two randomly overlapping finite width slabs clouds.The NWP cloud liquid water and cloud ice water content (CLWC,CIWC) profiles (in g/g) are integrated to obtain the column loading of the clouds (in g/m 2 ), and also to determine the slab cloud top/bottom pressures.The CIWC,CLWC profiles, cloud cover profiles and total cloud cover are used to determine the slab cloud fractions.Effective particle sizes then need to be assigned to the clouds.
Infrared sensors cannot see through optically thick clouds and are mostly sensitive to the emission from cloud upper boundary while emission throughout the cloud can contribute to the outgoing radiance for less optically thick clouds.The TwoSlab model is very flexible when placing the slabs, for example (a) at the weighted mean or centroid (C) of the cloud ice or cloud liquid profile or (b) near the most prominent cloud profile peak (P), which is best for optically thin clouds.
In practice, the cloud content profile CXW C(z) (where X = I, W for ice or water cloud) is smoothed before construction of the two cloud slabs (ice and water).The NWP cloud profiles usually result in the code identifying one water and one ice cloud, though the CLWC/CIWC pairs could produce two liquid or two ice slab clouds, or just a single slab cloud.
Figure 1 shows two examples of slab cloud outputs.The left panel is the simpler case where the NWP cloud profile (in blue) is singly peaked.The right panel shows a case where the profiles are much more complex, and the cyan (water cloud) slab is placed higher in the atmosphere closer to the peak of the weighting function.The integrated cloud amount g/m 2 is proportional to the width of the cyan/magenta bars.
Assuming one ice and one water cloud slab are produced, the cloud fractions are constrained as follows: where TCC is the NWP total cloud cover (TCC).From this we compute the individual clouds fractions and their overlap using the following criteria : 1.If there is only one cloud present, its cloud fraction is set to TCC.
2. For two ice or two water clouds, the fraction for the first cloud c 1 is set as T CC × R where R is a random number between 0 and 1.Then check to see if TCC is less than a random number; if so set c overlap to be c 1 × RR where RR is also random, else the overlap fraction is set to 0. c 2 then follows from Eq. 1.
3. If there is one ice and one water cloud, the cloud fractions are set according to c water , c ice (described below); after that c overlap is set using Eq. 1.
For the third case, the water cloud slab fraction comes from weighting the NWP cloud cover (CC) profile using the cloud liquid water content profile c water = CLWC(z) × CC(z)/ CLWC(z).The ice cloud fraction c ice is similarly determined.
Water particle effective diameters are 20 µm plus a uniformly distributed random offset.The ice effective particle size are estimated from a temperature based parametrization by (Ou and Liou, 1995), where the NWP temperature profile is used to associate the ice cloud slab top pressure with a cloud top temperature.

Radiance computation
The channel all-sky radiance r i (ν) is computed using four weighted radiance streams where f clr is clear fraction, cx i , i = 1, 2 is the exclusive cloud type i fraction and c overlap is the cloud overlap between the two cloud types; the exclusive cloud fraction being related to the cloud fraction via the relationship cx i = c i − c overlap .The model currently exclusively uses ice or water clouds when computing the radiances r (1) i (ν) associated with the cloud types; r clr i (ν) is the clear sky contribution while r (12) i (ν) is the radiance contribution from the cloud overlap.Since the atmospheric gas optical depth computation dominates the run time, computing four radiance streams is not a speed penalty, and the overall average run time per profile is about double that for a single clear sky radiance computation.

Maximum Random Overlap Conversion
The MRO cloud processing for the PCRTM model is described in Chen et al. ( 2013), and will only be briefly summarized here.
MRO converts the NWP water and ozone levels profiles to 100 layer profiles.For each layer, the cloud ice water content and cloud liquid water content mixing ratios are converted to a cloud optical depth.The optical depths at each layer are summed.
Layers above 440 mb are considered ice clouds, and layers in the lower atmosphere are assigned to water clouds (Rossow andSchiffer, 1983, 1991).The effective water diameter is set at 20 µm while the effective ice diameter is again temperature dependent, based on the parametrization in Ou and Liou (1995).The cloud cover profile cc(z) is used to generate 50 subcolumns using MRO (Chen et al., 2013) for which one radiance is computed per sub-column; the final radiance is an average over these sub-columns.

Clear-Sky Comparisons
An earlier inter-comparison of the SARTA and PCRTM clear sky models is presented in Saunders et al. (2007).In this subsection we assess the more recent spectroscopy embedded in the SARTA and PCRTM codes, using ECMWF thermodynamic Figure 2 shows the calculated BT biases between SARTA and PCRTM clear sky models along with AIRS noise levels.The top panel shows the mean differences, while the bottom shows the standard deviations.The mean bias between SARTA and PCRTM clear calculations is within AIRS noise levels at all channels, except in the methane region (1300 cm -1 ) and some channels in the water vapor 6.7 µm region.This is due to differing methane and water vapor spectroscopy and continuum models in these two RTAs.In addition, PCRTM uses a density weighted layer temperature that may introduce differences.The PCLSAM algorithm approximations in SARTA are more accurate for absorptive clouds that are more likely in the mid-IR.However, the DISORT-based scattering in PCRTM is more accurate if the cloud representation is correct.In general it would be reasonable to expect the differences to increase with optical depth and/or cloud fraction.In addition, in the TIR the single scattering albedo of water clouds is generally larger than that of ice, so we would also expect larger differences for water clouds.
To evaluate the SARTA (PCLSAM) versus PCRTM (DISORT) radiance differences, we used 1000 scenes spanning all climate types from AIRS on 2009/03/01.After matching the thermodynamic and cloud NWP fields to the observations, and subsequent conversion of the input cloud profiles to slab clouds, SARTA and PCRTM were run twice: (a) A clear sky run where no cloud effects are included, and (b) A all-sky run using the TwoSlab cloud representation derived from ECMWF.
In the TIR window region cloud forcing (difference between observed BT and surface temperature) can be as large as 100 K (for the Deep Convective Cloud (DCC) cases).For our sample set the mean (AIRS observation-SARTA TwoSlab RTA simulation) difference at 820 cm -1 is -1.8K compared to a mean cloud effect of 27.9 K. Similarly at 1231 cm -1 the mean difference is -2.4K compared to a mean cloud effect of 28.7 K.The corresponding biases for the PCRTM TwoSLab simulations are -2.4K at 820 cm -1 and -2.6 K at 1231 cm -1 .The effects of the clouds becomes less noticeable for channels sensing high in the atmosphere, such as the 650-700 cm -1 and 1400-1600 cm -1 regions.These comparisons shows the differences between the scattering RTAs is much smaller than the mean cloud effect.
Figure 3 shows the mean and standard deviations between SARTA TwoSlab/PCRTM TwoSlab using double differences, where the mean of the clear sky differences is removed from the mean of the all-sky differences.The double difference removes residual spectroscopic differences, allowing one to attribute the remaining differences to the RTA scattering algorithms.
Differences are seen in the window where cloud scattering is most significant, but overall they are less than 0.5K.A detailed examination showed that the differences for ice clouds was proportional to the ice cloud fraction, while there was comparatively more scatter in the differences between PCLSAM and DISORT for water clouds at any cloud fraction.
These comparisons indicate that our implementation of the PCLSAM model is a fast, yet simple and effective method to include scattering effects in the TIR, as has also been shown by Matricardi (2005); Vidot et al. (2015); Liuzzi et al. (2016).Plots of the 1231 cm -1 all-sky calculations using PCRTM/MRO are very similar to Figures 4 and 5, namely large areas of the tropical oceans having almost zero bias, and with noticeable mis-matches of cold cloud tops.

Window Channel PDFs
Here we explore the similarities and differences between the observations, SARTA/TwoSlab, and PCRTM/MRO by examining the radiance Probability Distribution Functions (PDF)s and the scene dependence of the mean BT differences, again for the 1231 cm -1 AIRS channel.
Cloud mis-match errors will contribute to a significant portion of the standard deviation between observations and calculations, as is evident in from Figures 4 and 5 1.In the tropics the observations have more cold (DCC) scenes, also seen in (Shahabadi et al., 2016).example (Shahabadi et al., 2016)).The differences between SARTA (P/C) and PCRTM-MRO are generally small compared to the observation versus model differences.
Note that for the highest clouds SARTA (C) is hotter than SARTA (P), which makes sense since these are very optically thick clouds and the TIR weighting function peaks at the cloud top.For warm scenes, the MRO simulations are again closer to SARTA (P) than SARTA (C) indicating that, as expected, placing clouds near the weighting function peak in the TwoSlab algorithm is preferable to the centroid.-4.2 ± 9.5 -3.5 ± 9.9 -5.4 ± 9.6 0.7 ± 4.4 -1.2 ± 3.5 -1.9 ± 3.5 N. Polar 6.9 ± 9.4 -2.9 ± 7.1 -1.2 ± 7.8 -0.3 ± 8.1 1.7 ± 4.5 2.6 ± 5.2 0.9 ± 5.2 S. Polar 9.9 ± 9.2 -4.9 ± 6.5 -2.8 ± 8.0 -3.5 ± 7.3 2.1 ± 5.0 1.4 ± 4.3 -0.7 ± 4.3 We note the PDF correlations between observations and all the calculations, and also among the calculations themselves is typically 0.9 and above, which reinforces the point that the NWP fields from ECMWF capture much of the atmospheric variability that is observed; however the mis-match between observed and model cloud tops lead to biases on the order of 2-4 K with standard deviations on the order of 10 K.
The implications of this are that by shifting the position of the model clouds, one could significantly mitigate the differences, and hence have a-priori micro-physical properties of the TwoSlab clouds, for a physical single footprint all-sky retrieval.This idea is exploited later in the section on applying the TwoSlab code for use in single pixel all-sky retrievals.We close this sub-section with Figure 7 where the BT1231 observations and calculations are plotted as pdfs for the tropics, mid-latitudes and polar regions.Note that we have limited the pdfs to contain data only over non-frozen oceans, using the ECMWF sea ice fraction model field.As in Figure 6 the black curves are observations, while SARTA/TwoSlab (P)eak and (C)entroid calculations are in (blue/cyan) and PCRTM/MRO calculations are in red.While the tropical PDFs are quite similar between observations and both RTAs, the mid-latitudes suggest ECMWF produces more clouds than observed.The polar PDF calculations are again quite similar, but with even large disagreements with observations.Accounting for this in detail is not the focus of the paper, but some insight was gained by comparing the cloud fields in Granule 001 obtained over the Southern Ocean/Antarctica (SOA) against those in Granule 137 which was filled with Marine Boundary Layer (MBL) clouds off the coast of Namibia.For the SOA granule the mean CIWC and CLWC were 1 × 10 −5 g/g and 0.3 × 10 −5 g/g with peaks centered at 825 ±50 mb and 875 ±50 mb respectively, and mean cloud fractions of 0.3 at about 800-900 mb, and less than 0.1 above that.
Conversely the mean CIWC amounts for the MBL granule were almost 20 times larger at roughly the same vertical position (while the CIWC amounts were 10 times lower and cloud fractions were about the same); however the computed radiances were about 5-10 K cooler than the surface temperatures, in rough agreement with the observations.The mean surface pressure and temperature for the SOA granule was 985 mb and 273 K respectively, while the mean atmospheric temperature at 850 mb was 262 K.The bottom-most panel of Figure 7 shows the observed peak was close to 260 K compared to the computed peak at about 270 K.All this evidence points to one of two possibilities about the polar over-ocean clouds in the NWP model inability to statistically reproduce the observations.They could either at too low an altitude, or at the right altitude but not optically thick enough or with low cloud fractions.

Spectral Comparisons
Finally we compare the observed and calculated spectra for ocean scenes where the ECMWF model sea surface temperature is very accurate and the emissivity is well known.As in Figure 7, we divide these observations into tropical, mid latitude and polar zones (using boundaries at ±30 and ±60), with over 200,000 observations per zone.Regions with ice contamination (according to ECMWF) were removed to avoid uncertainties in emissivity for the simulations.As with the PDF plots, the spatial mismatches mis-matches between modeled and actual clouds will still lead to significant differences between observations and calculations, as seen in from Figure 4.
These comparisons are summarized in Figure 8.In each panel the blue, green and red refer to the tropical, mid-latitude and polar regions; the solid lines (marked (S)) refer to the SARTA/TwoSlab calculations, while the lines with 'X' (marked (P)) As expected from studying the 1231 cm -1 PDFs, the largest differences between the observations and calculations are in the polar (red) region.The tropical biases were typically the smallest, though their standard deviation is the largest (as the models and calculations have to span warm surface temperatures all the way to cold DCC cloud tops).The largest spectral biases are in the window regions, which is to be expected as cloud effects are most readily seen in these regions.
The above plots are similar to all-sky monthly global averaged biases seen in (Shahabadi et al., 2016).They are put into context by considering the cloud forcing effect, which we simply define as the mean BT1231 difference between clear sky calculations and observations.In the tropics/mid-latitudes/polar regions they are 7.4 K, 10.2 K and 12.6 K respectively.Conversely the PCRTM/MRO versus SARTA/TwoSlab differences are a factor of 10 smaller, at -0.78 K, -0.78 K and 1.89 K respectively.Again our main emphasis here is to validate the accuracy of our simple SARTA TwoSlab algorithm relative to the more rigorous PCRTM-MRO RTA, not to evaluate the accuracy of the ECMWF model clouds.
6 Application of TwoSlab Code : All-sky Retrieval One full AIRS granule (6 minutes, 12150 spectra) worth of all-sky retrievals using the SARTA/TwoSlab code are provided here using the OEM method, which provides natural diagnostics such as DOF and AK that are extremely helpful in understanding the information content of our retrieval approach in the presence of highly variable clouds.This is a "proof of concept" retrieval which has been separately tested on several days worth of AIRS observations.While considerable efforts have been put into selection of various regularization matrices and channels, we anticipate additional fine-tuning in the future.

OEM Approach
We follow normal OEM notation here, where the observation vector y(ν) (brightness temperature) is modeled by the radiative forward model operator F and (ν) is the combined instrument and forward model noise, where ν is the wavenumber, and x is the thermodynamic and cloud state.The solution is regularized with matrix R using a cost function J (Rodgers, 2000;Steck, 2001) given by The first two terms are the observation and background penalties respectively, while the last is a constraint to reduce the amount of humidity supersaturation (Phalippou, 1996;Deblonde and English, 2003).Our regularization matrix contains both empirical regularization (Tikonov) and a-priori covariance-based terms.J is minimized using a nonlinear Gauss-Newton iterative approach.(Rodgers, 2000;Steck, 2001) where K is the jacobian.At present the observations covariance matrix S −1 is diagonal, combining a linear sum of forward model error (≤ 0.2 K per channel) and AIRS channel dependent noise error.The retrieval currently uses no AIRS shortwave channels (past 2000 cm -1 ) but uses almost the same 500 longwave channels used in the AIRS L2 retrieval.

State vectors and OEM parameters
The state vector consists of surface and profile temperatures (Kelvin), the logarithm of the water vapor profile (molecules/cm 2 ), as well as a (logarithmic) multiplier for the ozone profile and two slab cloud amounts.The jacobian matrix K in Equation 5 is built using these parameters.For any FOV the iterations were halted after i max = 5 iterations, or if there was no improvement in the χ 2 after iteration i ≤ i max .
The regularization matrix R is block diagonal for the temperature and water vapor profiles, and is a linear combination of altitude dependent covariance matrices (with exponential decaying off-diagonal elements) and a L 1 type first-derivative Tihkonov smoother.Diagonal terms were added for the remaining state variables being retrieved.The humidity saturation penalty is of the form J sat = i=N i=1 J i , where J i is set to 0 if the relative humidity (RH) of the i -th layer is less than 100 %, else it is computed from r(log 10 RH(i) 100 ) 3 where r = 100.0.For this paper we start with smoothed climatological profile, and use 2 K uncertainties for the temperature and surface temperature, 60% uncertainty for water vapor profiles and 10% cloud loading uncertainty.We start with ECMWF surface temperatures since they are likely to be better than climatology.Land and ocean surface emissivity is set from a database (see Masuda et al. (1988); Zhou et al. (2011)).

Retrieval a-priori
The highly non-linear effect of clouds on infrared radiative transfer makes retrieval success highly dependent on an accurate linearization point for cloud parameters.This fact has made it difficult to create physically-gased (not statistical) robust singlefootprint infrared hyperspectral retrievals.Since typical infrared all-sky retrievals only have 2-4 degrees of freedom for cloud fields (see Appendix 10) it is essential that the linearization point and a-priori for clouds be as accurate as possible.Fortunately, NWP model forecasts such as from ECMWF provide reasonably accurate cloud fields derived using the best physics possible for an operational model.However, as we have shown earlier, perfect spatial placement of model clouds is not possible, and winds can move forecast clouds significantly during the ±1.5 hour time difference between observations and an ECMWF analysis/forecast file.
The ECMWF cloud fields are statistically quite accurate in the sense that they reproduce similar spatial distributions of window channel brightness temperatures as the AIRS observations.Our approach is to compare the observed window channel brightness temperatures to those we compute from the ECMWF model (and cloud) fields using our TwoSlab approximation.
We then find the closest spatial match between a particular window region observation and nearby simulated window channel radiance (in a least squares sense).The cloud fields from the ECMWF grid point with the closest matching radiance are then used as the linearization point for the retrieval for the AIRS scene.We only retrieve cloud loading after the linearization, while keeping particle effective size, cloud top and bottom pressure and cloud fractions unchanged.The temperature and water vapor profile linearization point and a-priori is not taken from ECMWF, instead a monthly climatology is used.For example, Figure 5 shows the overall spatial agreement between observed and simulated clouds over the Tropical West Pacific, the ECMWF cloud fields are often offset from the observations (and are a factor of two fewer).The 300 hPa ECMWF (u,v) wind fields for this granule suggests that these convective regions are moving westward at approximately 36 km/h.Since the forecast is approximately 0.8 hours previous to the AIRS overpass, the model clouds could move by up to 30 km after the model forecast, or roughly two AIRS FOVs.Thus, just the time delays between forecast and observations can contribute to the cloud mis-match.
The temperature and water vapor profile a-priori state for the tests reported here used using co-located 10-year monthly averaged AIRS L3 fields (March 2004(March -2013)).The water vapor a-priori uncertainty was set to 60% and the temperature uncertainty to 2 K while the surface temperature used a 0.3 K uncertainty.
We chose a climatology for the a-priori in order to more easily demonstrate the performance of the retrieval algorithm.If a single-footprint retrieval was used for production of a long time series of AIRS Level 2 products we believe a reanalysis (ECMWF, MERRA) would be a more suitable a-priori in that it would provide accurate profile estimates below thick clouds.
The OEM framework will naturally provide this capability since the degrees-of-freedom of the retrieval shrink rapidly as the cloud thickness increases.The red curve just shows that we can find and include nearby clouds in the ECMWF model that have similar brightness temperatures to what is observed.

Single Granule Case Study
This section focuses on AIRS granule 039 from March 11, 2011, a day scene over the Tropical Western Pacific containing many DCCs.Figure 9 maps the 1231 cm -1 BT observations.The three lines marked (A), (B) and (C) are at AIRS scan angles of roughly -23, 0 and +24 • .Figure 10 shows the observations (black), retrievals (red) and original ECMWF calculations (blue) for these three lines.The significant offsets in the original (ECMWF) clouds is apparent by comparing the blue and black curves.The red curve shows how well the final retrieval simulates the window temperatures, due to a combination of moving the ECMWF clouds fields to the appropriate pixels together with the adjustment of the cloud fields during the retrieval.The final calculations largely reproduce the observations.A more thorough examination of (longitude,latitude) versus BT1231 scatter plots after the retrieval show that the spatial patterns after the retrieval, including the cold DCC, correlate extremely well to the patterns seen in Figure 9.

Spectral Biases and DOF
A first step in testing retrieval performance is to examine the final retrieval residuals and their standard deviations shown in Figure 11.Also shown are the initial differences between all-sky radiances simulated with ECMWF (clouds in their original positions) and the AIRS observations (blue curve).Most of the bias relative to ECMWF is due to clouds and some significant water vapor differences.Our retrieval BT bias and standard deviation are given by the red curve.(The standard deviations from the original ECMWF co-located cloud fields have been multiplied by 0.2 to fit on the graph.)Most channel biases are in the 0.1-0.2Krange although larger biases are evident in the 650-700 cm -1 stratospheric region.Further work is needed to determine the cause of these deviations.We have included all retrievals here, including those with low degrees of freedom where the a-priori is weighted heavily.
The retrieval standard deviations are far smaller than those for ECMWF in the region of strong water lines past 1400 cm -1 .
In this region cloud cover is unimportant except for the deepest DCCs, so this reduction in the standard deviations indicates skill in the retrievals.This is discussed in more detail below.addition, although ECMWF is quite accurate globally, it is difficult to judge how well ECMWF represents truth for a single granule but it is likely far better than our a-priori climatology.
These statistics have been separated by the total number of DOFs in the retrieval.Figure 13 shows scenes with DOFs less than 6.5 (thick cloud, 76 FOVs), and while Figure 14] have DOFs greater than 9 (almost or completely clear, 8165 FOVs).
The small number of scenes with low DOFs also limit the utility of their statistics.The ECMWF and climatology profiles were multiplied by the retrieval averaging kernels for these comparisons.
-4 -2 0 2 4 T (K) The low DOF case (Figure 13) shows little difference between retrieved profile temperature and a-priori in the lower troposphere.This is not too surprising for a tropical granule.The cloud contamination that caused these low DOFs led the retrieval to stick to the a-priori in the troposphere.Small movements in the upper-troposphere and stratosphere did occur that gave similar disagreements with the a-priori and ECMWF.The same is true for water vapor for the low DOF case, any differences between the retrieval and either the a-priori or ECMWF are in the upper troposphere, where the cloud contamination is lower.
The high DOF cases in Figure 14 indicate that the retrieval moved slightly away from the a-priori towards ECMWF at almost all levels.Results are a bit more mixed for water vapor.The retrieval standard deviation from ECMWF is lower than with the a-priori in most of the troposphere.
A much more definitive diagnostic of retrieval performance is given in Figure 15, which is a curtain plot of relative humidity along the line denoted by "B" in Figure 9.Here we compare a series (from top to bottom) of our retrieval, ECMWF, AIRS L2 retrieval and a-priori humidity profiles, with the vertical axis being pressure and horizontal axis being latitude.The bottom panel is a plot of the BT 1231 cm -1 channel which is an excellent proxy for cloud top height and opacity.The a-priori and ECMWF water vapor structures have been multiplied by the averaging kernels.

10
There is little water vapor structure in the a-priori (fourth panel) given that is a monthly average of many years.The retrieved water vapor (top panel) shows significant structure along this track, with some small instabilities in the regions of thick high The ability of the retrieval to catch some of the upper-troposphere variability near 150 hPa seen in ECMWF indicates good vertical resolution as well.Note that the retrieval does not use any information from ECMWF, except for the clouds.We also comment that the UTLS humidity from the AIRS L2 is significantly lower than either ECMWF or our retrieval, with the blanked out areas indicating where the surface AIRS L2 QA flags were bad.
The higher tropopause RH that was initiated by climatology remained unaffected by the retrieval; this could be alleviated by an improved first guess of the thermodynamic state, as well as choosing WV channels that peak very high in the atmosphere.In the future we plan to use reanalysis as a-priori, which will be adjusted by the retrieval when there are low to medium optically thick clouds.The use of a fixed shape for the ozone profile is also a limitation of the present retrieval that will be removed in later versions.

Cloud Parameter Changes
Comparisons between the initial EMCWF TwoSlab cloud parameters (found by matching window BTs to nearby ECMWF scenes) and the retrieved cloud parameters have a number of understandable differences.It is well known that NWP models do not produce as many deep convective clouds as observed so it is understandable that the mean ice cloud fraction changed from less than 0.5 to higher values ranging from 0.6 to 0.9 (these can be quite thin ice clouds).The water cloud fractions increased slightly from when less than 0.3, and generally decreased for the higher factions.In addition, the frequency of high ice cloud tops (less than 250 mb) increased while the rest decreased; for water clouds the largest increase in frequency of occurrence was seen between 500-700 mb.
A quick validation of our ice cloud optical depths is achieved by comparing of AIRS L2 ice optical depths versus our retrieved ice cloud loading (in g/m2).While the thermodynamic AIRS profiles are generated at 45 km (3x3 AIRS footprint) resolution, the ice cloud AIRS L2 retrievals are generated for single footprints (Kahn et al., 2014) after the L2 thermodynamic retrievals, in a separate step that keeps all other retrieval variables constant.Figure 16 is a comparison of AIRS L2 ice optical depths versus our retrieved ice cloud loading (in g/m2), clearly showing a proportional relationship.This is very encouraging in that we are also retrieving the full thermodynamic state and water cloud parameters simultaneously.
The retrieval used cloud heights derived from matching to nearby ECMWF cloud fields.

Conclusions
A fast infrared radiative transfer algorithm with the ability to handle two scattering layers (from clouds, aerosols, volcanic dust) has been described and compared to a more sophisticated, and often slower approach (Maximum Random Overlap).
Our ultimate goal is to perform single-footprint retrievals with hyperspectral IR sounder radiances.In particular we wish to handle the very common case of two cloud layers (water, ice) in order to provide accurate, higher spatial resolution retrievals of temperature and water vapor (and other minor gases).This approach uses the observed radiances in the retrieval, rather than 32 However, if the a-priori cloud parameters are not sufficiently accurate, it can be very difficult for the retrieval to converge quickly, or at all.Our approach uses NWP model fields (here ECMWF) to initialize the cloud model fields.Four sub columns (at most) are needed to compute a radiance for one scene, which is a small speed penalty in fast RT models, where most of the time is spent in computing the atmospheric optical depths.The TwoSlab model can be an order of magnitude faster than typical implementations of MRO and has nearly the same accuracy, both in terms of mean spectral radiances and radiance PDFs.The spectral bias between all-sky AIRS observations and calculations are dominated by spatial location mis-matches between actual and forecast clouds.Both approaches used the ECMWF cloud fields, and in general both differed from observations similarly.For example, PCRTM/MRO is slightly more accurate in the tropics than SARTA TwoSlab, while the opposite is true in polar regions.However, the comparisons of RTA simulations to observations are both limited by the accuracy of the NWP model fields and especially by small spatial mis-matches between NWP and observed clouds.The larger errors of both RTA approaches in the polar regions indicate that ECMWF clouds have too low optical depths.
We demonstrated the feasibility of the SARTA TwoSlab approach by performing single-footprint retrievals using an AIRS tropical granule.Our approach to the retrieval cloud initialization and a-priori, which is key to the successful results shown here, was to use NWP cloud fields in the region of the footprint of interest based on matching simulated and observed radiances.
These matched cloud fields are then converted from N-layer NWP cloud fields to the two-layer SARTA Two-Slab.Together these steps provide a method for robust, and fast, retrievals from single-footprint all-sky hyperspectral spectra.
A major advantage of single-footprint retrieval using the OEM approach are the retrieval quality diagnostics that are provided within the OEM framework.We demonstrated that the retrieval DOFs are reduced in the presence of thicker clouds.However, we were able to reproduce much of the water vapor variability in ECMWF (assumed to be relatively accurate, partly because it agrees with our retrievals) when using a climatology for the water vapor and temperature a-priori.
Existing AIRS L2 retrievals fail in scenes with thick clouds and where the 3-by-3 set of radiances used for cloud-clearing are too homogeneous (which is not always in the case of thick clouds as seen in Figure 12).This is particularly troublesome for long-term climate studies in that the AIRS L2 sampling is incomplete and may alias certain climate variables, especially for water vapor where microwave retrievals (that are part of the AIRS L2 system) have more limited value.
The retrieval approach examined here may be able to address sampling limitations of the existing AIRS retrievals by (1) using single-footprint retrievals that are not affected by cloud-clearing failures for highly homogeneous scenes, and (2) using a-priori information in a statistically correct way under condition of thick clouds.A possible approach is to use a reanalysis for the a-priori rather than climatology in order to insert the best possible information in cases where the DOFs of the retrieval are very low.
This work does not represent a rigorous analysis of the accuracy of our retrieval approach, but only a proof-of-principle that the technique appears viable.In particular, the temperature retrievals are not stressed in a tropical environment, although our results suggest significant skill for water vapor.The retrieval tests shown here were mostly all over ocean where the surface Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-261Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 28 July 2017 c Author(s) 2017.CC BY 4.0 License.We extended SARTA to handle clouds and aerosols, based on Parametrization of Clouds for Long-Wave Scattering in Atmospheric Models (PCLSAM) Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-261Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 28 July 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 1 .
Figure 1.Example of cloud vertical profiles, reduced to one or two slabs.The red and blue curves come from the NWP model, while the cyan and magenta are the resulting locations (and loadings) for the slabs.
profiles and surface parameters to compare clear sky radiances computed from the models.We use 1600 randomly chosen night time scenes observed by AIRS on 2009/03/01 for an inter-RTA clear sky simulation comparison.The locations span all climate zones over ocean and land, as well as all AIRS scan angles.Night time scenes are used to avoid non-local thermodynamic equilibrium (De Souza-Machado et al., 2007) and solar surface reflectivity during the daytime in the 4 µm shortwave region.Both of these effects are handled differently by SARTA and PCRTM and are not relevant to this paper.Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-261Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 28 July 2017 c Author(s) 2017.CC BY 4.0 License.

5Figure 2 .
Figure 2. Spectral differences for clear sky calculations between SARTA and PCRTM for ocean night scenes.Also shown are the AIRS noise levels.Typical difference between PCRTM and SARTA are less than 0.25 K, except in the methane region (1300 cm -1 ) where SARTA does not use line mixing, and in some of the strong water vapor lines (1400+ cm -1 ).

5Figure 3 .Figure 4 .Figure 5 .
Figure 3. Biases and standard deviations between 1000 TwoSlab computations using SARTA and PCRTM.See the text for details on these double-difference signals.
. Conversely as mentioned earlier, the dynamic range spanning 200K to 330K is seen both in observations and calculations from both cloud representation models.The left panel of Figure 6 is a plot of the corresponding histograms or un-normalized probability distribution functions (PDFs); the bins are 1 K wide.The curves are the night time observations (black), SARTA/TwoSlab calculations (blue/cyan) and PCRTM/MRO calculations (red).The Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-261Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 28 July 2017 c Author(s) 2017.CC BY 4.0 License.blue and cyan SARTA TwoSlab calculations differ in the positioning of the slabs: the (P)eak of cloud weighting function or the (C)entroid of the cloud profile respectively, as described in Section 3.1.The calculated radiance histograms are very similar compared to their differences with the observations.If the 1231 cm -1 PDFs are subset for different geographical regions discrepancies between the computed and observed PDFs are easier to evaluate compared to the global PDFs.For example:

5 2 .Figure 6 .
Figure 6.Left: 1231 cm −1 channel brightness temperatures PDFs (probability distribution functions) for night scenes (global).Locations (1),(2) and (3) highlight large differences that are discussed in the text, while the horizontal lines and associated circles in Location (4) represent the BT1231 cm −1 means and standard deviations.Right: Mean biases as a function of scene temperature for various RTAs tested here.Typical standard deviations between the different RTAs is about 10K, while the typical standard deviations between observations and calculations vary from 30K (cold scenes) to 10 K (warm scenes).

Figure 7 .Figure 8 .
Figure 7. Brightness temperature PDFS for the night non-frozen ocean 1231 cm -1 channel for the data used in Figure 8 separated tropical, mid-latitude and polar regions.The whiskers correspond to the brightness temperature means and standard deviations.

Figure 12 Figure 11 .
Figure 12 maps retrieval DOFs for this granule.Comparisons to Figure 9 show the clear correlation between low window channel BT, due to clouds, and low retrieval DOFs.The DOFs range from 5 for thick DCC (ice cloud loadings of 100+ g/m2)and increase to larger than 10 for nearly clear scenes.These DOF values are more than a factor of two smaller than what is reported in Appendix B, where a diagonal-only error covariance matrix is used when assessing the N-level cloud information content.The red circles in the figure are the 60% of the observations in this granule where the L2 retrieval failed all the way down to the surface.6.2.2 Thermodynamic retrievalsTemperature and water vapor retrieval statistics are shown in Figures13 and 14.Statistical measures for a single granule, especially in the tropics, are of limited value for understanding retrieval accuracy, especially for temperature.However, they do help indicate nominal performance and provide a measure of the impact of clouds on the DOFs and retrieval accuracy.In

Figure 12 .
Figure 12.Retrieved degrees of freedom for G039 show evident dependence on observed BT1231 (which depends on cloud loading).The red circles denote locations of the AIRS FORs (3x3) where the retrieval quality down to the surface is 2, meaning missing or do not use; this happened for 60% of the L2 FORs.

Figure 14 .
Figure14.Same as Figure13except now uses cases for thin clouds, defined by DOFS to be > 10; 6220 profiles were used here (compared to 7296 if we ignore the AIRS L2 QA and only use high DOF retrievals, again with very similar plots being obtained).

Figure 15 .
Figure 15.Comparing Relative Humidity (colorbar) along track "B" in Figure 9. From top to bottom we have our retrieval, ECMWF, AIRS L2 retrieval and a-priori humidity profiles.The bottom panel is a plot of the observed BT 1231 cm -1 channel.The black and red circles in the top panel are the positions of ice and water clouds, with the circle size denoting the OD.29

Figure 16 .Figure 17 .
Figure 16.Retrieved ice cloud amount (in g/m2), compared to AIRS L2 ice cloud optical depths.The colorbar is the log(10) of the number of points in the bin.
Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-261Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 28 July 2017 c Author(s) 2017.CC BY 4.0 License.derived equivalent cloud-cleared radiances that are presently used for the NASA AIRS Level 2 products.The complexity of true cloud structures cannot be retrieved with hyperspectral IR radiances, and we have shown that only a maximum of 2-4 cloud parameters can be derived from a single scene, suggesting that only a simple RTA is needed.

Table 1
summarizes the 1231 cm -1 window channel global statistical comparisons including differences with observations and differences among the various RTAs.First note that the high standard deviations for observations minus computed radiances are dominated by spatial/time mis-matches, and are not necessarily indicative of model limitations.For the purposes of this paper, the most interesting result is that the SARTA (P) agrees better with observations in the mean, while SARTA (C) generally agrees better with MRO.The uncertainties in the model cloud fields are sufficiently large (especially given spatial/time mismatches) that these comparisons are not sufficient to indicate which scattering model is more accurate.Table1.Night land and ocean 1231 cm -1 biases for 2011/03/11 in K. "P" and "C" denote the cloud slabs placed at the (P)eak of the weighting function and the cloud profile (C)entroid respectively.