## Abstract

A linear inversion algorithm to derive changes of surface skin temperature and atmospheric temperature and specific humidity vertical profiles using spatially and temporally averaged spectral radiance differences is developed. The algorithm uses spectral radiative kernels, which is the top-of-atmosphere spectral radiance change caused by perturbations of skin temperature and air temperature and specific humidity in the atmosphere, and is an improved version of the algorithm used in earlier studies. Two improvements are the inclusion of the residual and cloud spectral kernels in the form of eigenvectors of principal components. Three and six eigenvectors are used for, respectively, the residual and cloud spectral kernels. An underlying assumption is that the spectral shape of the principal components is constant and their magnitude varies temporally and spatially. The algorithm is tested using synthetic spectral radiances with the spectral range of the Atmospheric Infrared Sounder averaged over 16 days and over a 10° × 10° grid box. Changes of skin temperature, air temperature, and specific humidity vertical profiles are derived from the difference of nadir-view all-sky spectral radiances. The root-mean-square difference of retrieved and true skin temperature differences is 0.59 K. The median of absolute errors in the air temperature change is less than 0.5 K above 925 hPa. The median of absolute errors in the relative specific humidity changes is less than 10% above 825 hPa.

## 1. Introduction

Changes of atmospheric temperature and water vapor profiles, as well as surface and cloud properties (i.e., fraction, height, optical thickness, phase, particle size) at a global scale, as a response to increasing the carbon dioxide concentration are difficult to observe because the magnitude of natural variability is larger than the change of these properties occurring over a time scale equivalent to the lifetime of satellite instruments. In addition, even though these properties are currently observed from satellite-based instruments, the uncertainty in the calibration of these instruments usually makes it difficult to detect such changes within the lifetime of the instruments or the duration of the missions (Wielicki et al. 2017). Furthermore, uncertainty characterization in the retrieval using instantaneous cloud-cleared spectral radiances with a Bayesian approach for long-term and large-scale applications is complex (Smith and Barnet 2019). Another way to detect changes is to use spectral information caused by property changes. A possibility of detecting subtle atmospheric changes by spectral radiance change, known as spectral fingerprinting, is discussed by Goody et al. (1998), Leroy (1998), and Anderson et al. (2004). Generally, the spectral fingerprinting uses highly averaged spectral radiances and spectral signatures to separate the contribution by various atmospheric and surface properties. We refer this approach to the average-then-retrieve approach as oppose to retrieve-then-average approaches that derive atmospheric and surface properties using instantaneous spectral radiances. While it does not retrieve instantaneous properties, the average-then-retrieve approach offers several advantages compared to conventional retrieve-then-average approaches (e.g., Susskind et al. 2003, 2014; Chahine et al. 2006; Weisz et al. 2013; Irion et al. 2018; DeSouza-Machado et al. 2018; Shao and Smith 2019) in deriving property changes with time.

One of advantages of the average-then-retrieve approach is that once observed radiances are averaged, both instrument noise and natural variability are reduced, the larger the temporal and spatial averaging scales, the larger the reduction. The primal task in detecting the trend of atmospheric and cloud properties using averaged spectral radiances is separating cloud and atmospheric property changes contributing to the highly spatially and temporally averaged observed radiance change. The feasibility of such retrievals has been investigated in earlier studies using longwave spectral radiances by Leroy et al. (2008), Huang et al. (2010a), and Kato et al. (2011), shortwave spectral radiances by Feldman et al. (2011), Jin et al. (2011), and Roberts et al. (2011), and longwave spectral radiances combined with radio occultation data by Huang et al. (2010b). Pan et al. (2017) applied the technique to detect stratospheric temperature changes. Note that both retrieve-then-average and average-then-retrieve approaches are affected by instrument differences discussed by Smith et al. (2015), when multiple instruments are used over the time period.

In this study, we extend the work by Kato et al. (2014) and use computed (synthetic) spectral radiances to build a theoretical base of the average-then-retrieve approach. The cause of retrieval errors using averaged spectral radiances is different from the cause of retrieval errors using instantaneous radiances. The retrieval error depends on how well spectral difference are linearized and the size of the residual nonlinear term. Huang et al. (2010a) treat the nonlinear term as an uncertainty and include the uncertainty as a part of the covariance matrix. In this study, we treat the nonlinear residual term explicitly in the spectral radiative kernel matrix. To develop and formulate a robust algorithm to treat the residual term explicitly, we use 16-day mean spectral radiances in this study. In the flowing, section 2 discusses how synthetic spectral radiances and spectral kernels are generated. Instead of using monthly mean spectral radiances that were used in Kato et al. (2014), we use a shorter time scale to average spectral radiances in order to develop and test the method to account for the residual term caused by the linearization. As shown in section 3, the residual term is significant when 16-day mean spectral radiances are used. Section 3 also discusses the method to mitigate the impact of the residual term in the retrieval. Section 4 presents results and section 5 discusses and summarizes the error in retrieved skin temperature, air temperature and specific humidity profile changes.

## 2. Computation of synthetic radiances

Nadir spectral radiances are computed with a radiative transfer code developed by Liu et al. (2006). Temperature and specific humidity profiles are taken from the Modern-Era Retrospective Analysis for Research and Applications (MERRA; Rienecker et al. 2011). The spatial resolution of temperature and specific humidity profiles is 0.66° × 0.50° with 72 vertical levels in the atmosphere and the temporal resolution is hourly. Cloud fields are derived from *Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations* (*CALIPSO*; Winker et al. 2010), *CloudSat* (Stephens et al. 2008), and Moderate Resolution Imaging Spectroradiometer (MODIS; Minnis et al. 2008a,b) observations. Vertical cloud profiles are generated by the method described in Kato et al. (2010) and included in the *CALIPSO*–*CloudSat*–CERES–MODIS (CCCM) data product (Kato et al. 2011). Surface emissivities are from those derived from Infrared Atmospheric Sounding Interferometer (IASI) by Zhou et al. (2013).

Synthetic nadir-view spectral radiances are computed every near nadir view with the footprint size of approximately 20 km, following the Aqua orbit from 2007 through 2010. Geolocation of the center of footprints overlapping with the ground track of *CALIPSO* and *CloudSat* is included in the CCCM data product. The geolocation and time are used to sample temperature and humidity profiles from MERRA to extract matched profiles for every footprint. The wavenumber range considered in this study is from 649 through 1613 cm^{−1} and from 2181 through 2665 cm^{−1}, which is equivalent to the wavenumber range observed by Atmospheric Infrared Sounder (AIRS). Because cloud properties are taken from *CALIPSO* and *CloudSat* observations, those are viewed from slightly off nadir view (approximately 15° viewing zenith angle). Three cloud types per footprint and a clear-sky condition are used for all-sky spectral radiance computations. Clouds derived from *CALIPSO* and *CloudSat* are separated by three cloud types depending on the cloud-top pressure. Three cloud types are, high-, mid-, and low-level clouds. Clouds with the top pressure smaller than 400 hPa and greater than 700 hPa are classified, respectively, as high- and low-level clouds. Clouds with the top pressure in between are midlevel clouds. We compute the spectral radiance for every footprint using collocated temperature, humidity and cloud profiles depending on time and geolocation of the footprint. Resulting synthetic spectral radiances are averaged over 10° × 10° grids and over 16 days. For the purpose of development and testing the algorithm, we split all data into two time periods: the first half is from 2007 through 2008 and the second half is from 2009 through 2010. The first half is used to develop an optimal fingerprinting algorithm and the second half is used to evaluate the algorithm.

Nadir-view spectral radiances are also computed by perturbing surface, atmosphere, and cloud properties by the method described in Kato et al. (2014). Perturbed amounts are shown in Table 1 of Kato et al. (2014). We subsequently derive spectral radiative kernels by subtracting radiances computed with unperturbed properties from those computed with perturbed properties. Note that we only consider spectral difference due to water vapor variability and ignore the variability caused by other gases in this study (i.e., concentrations of carbon dioxide and other trace gases are fixed).

## 3. Spectral radiative kernels and linear inversion method

When surface, atmospheric, and cloud properties are perturbed from the mean value, the instantaneous nadir-view spectral radiance change *δ*$R$ can be expressed as

or

where $\delta R=R\u2061(t)\u2212R\xaf$, $R$(*t*) is the spectral radiance at time *t* corresponding properties of **a**(*t*), $R\xaf$ is the mean spectral radiance corresponding the mean properties $a\xaf$, $K$ is the kernel matrix, $e$ is the residual term caused by the error in $K$ or linearization, i.e., $e=\u2061[R\u2061(t)\u2212R\xaf]\u2212K\u2061[a\u2061(t)\u2212a\xaf]$, *i* and *j* are indices, respectively, for wavenumber and property, and overbar indicates the average over a time period. The *m* × *n* kernel matrix $K$ contains spectral radiance change at *m* wavenumbers caused by perturbations in *n* atmospheric properties (spectral kernel) and Δ**a** contains scaling factors to the perturbations from corresponding mean properties. Many atmospheric properties can influence the TOA longwave spectral radiance, including surface skin temperature, air temperature and specific humidity vertical profiles, cloud fraction, cloud optical depth, and cloud-top height.

When Δ$R$ is computed as the difference of spectral radiances from two time periods *t*_{1} and *t*_{2}, then

where $a\u2061(t1)\u2212a\xaf$ and $a\u2061(t2)\u2212a\xaf$ are the deviation of properties for time period *t*_{1} and *t*_{2} from mean properties defined as

*n*_{1} and *n*_{2} are numbers of observations, respectively, over time period 1 and time period 2 and overbar indicates temporal and spatial averaging. We rewrite Eq. (3) as

where

and

In Eq. (8), *δ*$R$ is defined as $\delta R\u2061(t)=R\u2061(t)\u2212\u2061[1/\u2061(n1+n2)]\u2061(\u2211n1R+\u2211n2R)$. The residual term $e\xaf$ in Eq. (5) is originated from two sources. The first term $e\xaf1$ is from the residual associated with linearizing instantaneous radiance change and the second term $e\xaf2$ is caused by using averaged spectral kernels. Both terms are affected by contributions from property changes that are not included in spectral kernels. Subgrid-scale spatial variabilities and temporal variabilities within the averaging period are responsible for $e\xaf2$. For example, the sensitivity of radiances to a perturbation of lower-atmosphere properties depends on upper-atmospheric properties. The mean upper-atmospheric properties multiplied by the mean lower-atmospheric property change does not necessarily give the same sensitivity computed by averaging instantaneous radiance changes (Kato et al. 2014). The mean spectral kernel defined by Eq. (6) is the average over two periods. However, we use $\u2061(1/n1)\u2211n1K\u2061(t1)$ to simplify the algorithm so that the difference $[1/\u2061(n1+n2)]\u2061[\u2211n1K\u2061(t1)+\u2211n2K\u2061(t2)]\Delta a\xaf\u2212\u2061[\u2061(1/n1)\u2211n1K\u2061(t1)]\u2061[a\u2061(t1)\u2212a\xaf\xaf]$ is included in $e$_{2}, i.e., $e\xaf2={\u2061(1/n2)\u2211n2K\u2061(t2)\u2061[a\u2061(t2)\u2212a\xaf]\u2212(1/n1)\u2211n1K\u2061(t1)\u2061[a\u2061(t1)\u2212a\xaf]}\u2212\u2061[\u2061(1/n1)\u2211n1K\u2061(t1)]\u2061[a\u2061(t1)\u2212a\xaf\xaf]$. Wavenumber ranges from 649 through 1613 cm^{−1} and from 2181 through 2665 cm^{−1} is used for the retrieval.

Top panels of Fig. 1 show the mean spectral radiance change $\Delta R\xaf$ computed with Eq. (3) by the red line and the change computed with $\Delta R\xaf=K\xaf\Delta a\xaf$ (Hereinafter, we call $K\xaf\Delta a\xaf$ the $K$-based or kernel-based mean spectral radiance changes) by the green line for clear sky conditions in four selected 10° × 10° regions. The difference between two lines is $e\xaf$ in Eq. (5). Plots of $e\xaf$ for all-sky conditions is shown in bottom panels. The figure shows that $e\xaf$ for all-sky conditions is much larger than $e\xaf$ for clear-sky conditions. Cloud profiles are highly variable and thus greatly contribute to spatial and temporal variabilities, which is the reason for a larger $e\xaf$ for all-sky conditions than $e\xaf$ for clear-sky conditions. Figure 2 shows the global distribution of the spectrally integrated residual $e\xaf$ in kernel-based mean spectral radiance change $K\xaf\Delta a\xaf$ compared with spectral radiance change Δ$R$ on all 10° × 10° grid boxes for clear-sky (top left) and all-sky (top right) conditions. The error is expressed as the ratio of root-mean-square (RMS) of $ei\xaf$ summed over the spectral range relative to the RMS of the summation of $\Delta Ri\xaf$. The top-right plot of Fig. 2 indicates that the magnitude of $e\xaf$ can exceeds 50% of the RMS of summed spectral radiance differences.

### a. Spectral kernel matrix

One of the key requirements to achieve a small retrieval error is that the sum of individual contributions approximately agrees with spectral radiance change. Results by Kato et al. (2014) and Thorsen et al. (2018) suggest that the agreement is relatively good when a longer time scale (e.g., monthly or longer) and a larger spatial scale (e.g., zonal) are used because changes are generally small. However, a significant difference exists when a shorter time scale and smaller spatial scale such as used in this study (16 days and 10° × 10° grid box) is used. Using a shorter time scale and a smaller spatial scale, we demonstrate the method to treat the difference between the sum of individual contribution and spectral radiance change in order to reduce the retrieval error in this study.

To reduce the retrieval error, we rewrite $e\xaf$ as $\u2061(K\xafe\Delta a\xafe+e\xaf*)$, where $K\xafe$ is spectral kernel representing the spectral feature of the residual $e\xaf$ and $\Delta a\xafe$ is a scaling factor. With this assumption, Eq. (3) is rewritten as

Explicitly including $K\xafe$ in the spectral kernel matrix is different from the algorithm used by Huang et al. (2010a) and Kato et al. (2014). Because the residual is large when clouds are present, as shown by Fig. 1, the variability of spectral radiance due to clouds that are not canceled out in the averaging is in part responsible for $K\xafe$. In Eq. (9) we also form the cloud spectral kernel $K\xafc$ that is the spectral radiance change due to cloud property changes. We then add both $K\xafe$ and $K\xafc$ to the original spectral kernel matrix $K\xaf$ to form a new spectral kernel matrix $K\xaf*$. Skin temperature, atmospheric temperature and specific humidity spectral kernels are included in $K\xafT,q$. We use the logarithmic derivative for specific humidity spectral kernels so that the vertical profile of the relative change of specific humidity is retrieved. The scaling factors to $K\xafc$ and $K\xafe$ are retrieved in addition to the perturbations in surface skin temperature and atmospheric temperature and specific humidity profiles. Both cloud and residual spectral kernels are formed in a similar way. In the next section, we describe the method to form the cloud and residual spectral kernels.

#### 1) Construction of spectral radiative kernel $Kc$ for cloud properties

Cloud profiles are highly variable spatially and temporally. In addition, the spectral shape of nadir-view radiance change caused by cloud property variability is not necessarily constant because the spectral shape depends on, for example, the water vapor amount above the cloud-top height. When the spectral shape of a fixed amount of perturbation for a given variable varies, but a constant spectral shape is used for the retrieval, it causes a problem in retrieving cloud property change in inverting Eq. (9). To express cloud variability in the form of spectral kernels and scaling factors, we derive cloud spectral kernels using cloud radiative effect, which is defined as all-sky minus clear-sky nadir-view spectral radiances, in order to isolate the spectral radiance change due to cloud property change from temperature and specific humidity changes.

The top plot of Fig. 3 shows all ensembles of mean cloud radiative effect in 16-day periods in 2007–08 on one selected grid box (20°–30°S, 70°–80°E). To show the spectral shape more clearly, we normalize them by dividing these cloud radiative effects by corresponding broadband radiances so that spectrally integrated value of the normalized spectral radiance equals to one. The normalization is used to show the similarity of the spectral shape in Fig. 3 only and it is not used in the actual algorithm. The bottom plot of Fig. 3 shows that spectral shapes of the normalized cloud radiative effect are similar, which facilitates the use of the cloud radiative effect in the linear inversion technique to retrieve cloud property changes.

To express the variability of the spectral cloud radiative effect, we apply empirical orthogonal function (EOF) decomposition onto the cloud radiative effect matrix $R$_{c} over each 10° × 10° grid box. The element of $R$_{c}, denoted by $Rcki$, is cloud radiative effect at the *i*th wavenumber averaged from the *k*th 16-day period in 2007–08. First 3 principal components are retained and cloud radiative effect is decomposed as $\u2211i\u2061(pc)i\u2061(ec)iT$, where pc_{i} is the *i*th dominant principal component and ec_{i} is corresponding eigenvector after the EOF decomposition. The first 3 principal components explain, on average, 99% of the variabilities of the cloud radiative effect for at least 91% of grid boxes. We subsequently use ec_{1}, ec_{2}, and ec_{3} as spectral radiative kernels in $K\xafc$ for cloud property changes, and projection of cloud radiative effect change onto $K\xafc$ is the corresponding cloud property changes $\Delta a\xafc$. Therefore, $K\xafc\Delta a\xafc$ is the response of kernel-based cloud radiative effect change caused by the changes in cloud properties. Individual cloud property variability such as cloud fraction and cloud optical thickness variability, is not separated using cloud kernel $K$_{c}. Elements of Δ**a**_{c}, therefore, no longer represent the physical property of clouds but they are a mixture of multiple cloud properties.

#### 2) Construction of spectral radiative kernel to reduce residual

After including spatially and temporally averaged spectral kernels of surface skin temperature, temperature profile, specific humidity profile, and cloud radiative effect in the spectral kernel matrix in Eq. (9), we determine the residual term $e\xaf$ using synthetic spectral radiances for the time period from 2007 through 2008. Once $e\xaf$ is derived, we apply EOF decomposition on $e\xaf$ to each grid box and retained first 6 principal components, similar to the EOF decomposition applied to the cloud kernel. The first 6 principal components explain 99.5% to 100% variability of the residual term for most of the 10° × 10° grid boxes. Corresponding the first 6 eigenvectors are used to form $K\xafe$. Using 6 eigenvectors allows changing the spectral shape of $e\xaf$ by changing the contribution from the 6 elements.

When $K\xafc$ and $K\xafe$ are formed, they are included in the spectral kernel matrix to form $K\xaf*$. The bottom-right plot of Fig. 2 shows the ratio of the RMS summation of $e\xaf*$relative to the RMS summation of $\Delta R\xaf$ for all-sky conditions (right). The ratio is less than 1% for most of 10° × 10° grid boxes, which is a significant reduction from the ratio shown in the top plots especially for all-sky conditions (top-right plot) of Fig. 2.

### b. Linear inversion

We derive relative changes of specific humidity and changes of surface skin temperature, air temperature profile, cloud, and residual terms $\Delta a\xaf*$ by inverting Eq. (9). In doing so, we add two constraints in the inversion. First constraint is the prior information,

Equation (10) defines $w\xaf\Delta a$ as the difference between the solution of the perturbation $\Delta a\xaf*$ in atmospheric and cloud properties and a priori information $\Delta a\xafp$. The second constraint is that scaling factors need to be smooth between neighboring layers,

where $G$ takes the difference of two elements of $\Delta a\xafp$ or $\Delta a*\xaf$ next to each other, which forces the change of atmospheric profiles to be smooth with height. Equation (11) defines $w\xafG$ as the difference between the smoothed solution of atmospheric and cloud perturbations and smoothed a priori information.

Following the method in Dubovik (2004), when these two constraints are included, the solution is

where $S$_{$R$}, $S$_{a}, and $S$_{$G$} are, respectively, temporal covariance matrices constructed from **e*****, $w$_{Δa}, and $w$_{$G$} as defined in Eq. (9), (10) and (11). Generally, only their diagonal components are retained when inversing $S$_{a} and $S$_{$G$}. Inversion of $S$_{$R$} is computed as $\u2211\mu =1k\lambda \mu \u22121s\mu Ts\mu $, in which *λ*_{μ} and **s**_{μ} are the first *k* eigenvalues and eigenvectors of $S$_{$R$}. We use *k* = 10 between 60°N and 60°S and *k* = 8 for the latitude greater than 60°. The result of inversion is not very sensitive to the choice of the *k* value. These covariance matrices are

where error matrices $e$*, $w$_{Δa}, and $w$_{$G$} are, respectively, *m* × *l*, *n* × *l*, and *n* × *l* matrices, and *l* is the number of sets of two time periods (i.e., 23 pairs of periods per year), overbar indicates average matrix $X\xaf=(1/l)IX$ and $I$ is a *l* × *l* matrix of which all elements are 1. Because we retrieve the difference of surface, atmosphere, and cloud properties averaged over two time periods, we use no property changes between two periods as a priori information, i.e., $\Delta a\xafp=0$. Therefore, the last two terms in the second parenthesis of Eq. (12) are zero. Temporal samples of $e$*** are collected from all differences between any two consecutive 16-day average of spectral radiances and atmospheric profiles from 2007 to 2008 over a 10° × 10° grid box. Covariance matrices for different grid boxes are formed separately. Note that our covariance matrices are different from the covariance matrix used in Huang et al. (2010a) who include the uncertainty due to natural variability, spectral kernel shape, and nonlinearity.

## 4. Results

The covariance matrices described by Eqs. (13), (14), and (15) are formed with the dataset from 2007 through 2008. We then apply the linear inversion of Eq. (12) to the synthetic spectral radiance differences from 2009 to 2010 in order to retrieve changes in surface skin temperature and temperature profiles and relative changes of specific humidity profiles between any two consecutive pairs of 16-day mean radiances averaged over a 10° × 10° grid box. The surface temperature and other atmospheric profiles used for computing synthetic spectra are treated as the truth to evaluate our retrieval results. The difference between true and retrieved values with one pair of 16-day mean radiances is referred as an error. Because the sign and magnitude of the error changes when different pairs of 16-day mean radiances are used, the error is a bias with unknown sign.

### a. Surface skin temperature

Figure 4 compares the retrieved surface skin temperature changes with true changes. We obtain 29 808 retrievals, 46 retrievals on each 10° × 10° grid box. Among all true and retrieved values, respectively, 35 and 44 values are outside 5 times the standard deviation. Therefore, we exclude 44 cases from Fig. 4 in computing statistics. Retrieved skin temperature changes agree with the truth relatively well. The correlation coefficient between them is 0.98 and the RMS difference is 0.59 K. The global distribution of mean absolute errors in the retrieved surface skin temperature change on each grid box is shown in Fig. 5. Retrieval errors are quite small compared with true differences. Larger retrieval errors mostly exist over land, especially at high latitudes where skin temperature changes are larger (Fig. 5, bottom). As a consequence, the relative error of surface skin temperature change is small. The relative error is less than 25% of true values for 52% of all retrievals and the mean error of 11% of grid boxes are less than 25% of corresponding mean true skin temperature changes.

### b. Temperature and humidity profiles

Figure 6 shows correlation coefficients between true and retrieved temperature and relative specific humidity changes for two consecutive 16-day means over 10° × 10° grid boxes as a function of height. The box and error bar represent, respectively, the spread of 50% and total populations of the correlation coefficient distribution including all grid boxes. The vertical bar in the box indicates the median. The figure shows that the inversion algorithm can retrieve temperature and relative specific humidity changes relatively well above the boundary layer. The larger error in the stratosphere is due to a small absolute amount of water vapor in the layers and its small change (e.g., Huang et al. 2010a). Huang et al. (2010b) show that retrieval of stratospheric water vapor change can be improved when radio occultation data are combined. The correlation coefficient between true and retrieved temperature changes above 874 hPa is larger than 0.9 for most of the grid boxes. In addition, the correlation coefficient of relative specific humidity changes above 748 hPa is larger than 0.9 for almost all grid boxes. In the boundary layer, the correlation coefficient between the true and the retrieved values are significantly lower and vary depending on region. The poor retrieval accuracy in the boundary layer could be caused by weaker spectral signal from lower-tropospheric changes due to cloud cover. In addition, the similarity of spectral kernels of temperature, skin temperature, and low-level cloud fraction changes might also contribute to the larger error. Using spectral radiances averaged by selecting clear-sky and nonoverlapping clouds separately might be needed to reduce the error in boundary layer temperature and specific humidity.

Figure 7 shows the absolute error of temperature and relative specific humidity changes. Because a few large retrieved values distort the mean, Fig. 7 shows the full distribution with the error bar, 50% populations with the box, and median with vertical bar in the box. Similar to skin temperature results, outside 5 times the standard deviation cases are excluded. True and retrieved values that are outside 5 times the standard deviation are less than 1% of all retrievals at all heights. Similar to the result of correlation coefficients, absolute retrieval errors are small above the boundary layer. The median of absolute errors in the temperature change is less than 0.5 K except for the lowest three layers. The median of absolute errors in the relative specific humidity changes is less than 10% above 825 hPa (top plots of Fig. 7).

In addition, the relative mean error is computed in two different ways. The first relative mean is computed by averaging the absolute error divided by the absolute value of the truth

where *δx* and Δ*x* are, respectively, the retrieval error and true value and *n* is the number of retrievals. The second relative error is computed by the mean of absolute errors divided by the mean of absolute values of true change

When there are small values of Δ*x* associated with large *δx*, $r\xaf1$ may not represent the mean. Similarly, when there are large values of Δ*x* associated with small *δx*, $r\xaf2$ may not represent the mean. When all Δ*x* are nearly equal, the means computed by Eqs. (1) and (2) are nearly the same. The median of $r\xaf1$ is between 20% and 30% for temperature changes and 30%–50% for relative specific humidity changes (top plots of Fig. 8), while $r\xaf2$ is within 40% of mean of absolute values of true temperature change and within 50% of mean of absolute values of true relative specific humidity change above the boundary layer (bottom plots of Fig. 8). Large differences between $r\xaf1$ and $r\xaf2$ within the boundary layers shown in top and bottom plots of Fig. 8 indicate that true values are highly variable in the boundary layer. This is consistent with Fig. 7, showing that the upper end of the distribution of true values is large. These extreme cases are expected to disappear as changes between two periods become smaller when we use a longer temporal and larger spatial mean compared to those used in this study. Retrieval results are expected to improve because the change of the spectral radiance is expected to be more linear (i.e., the residual term is smaller) in a longer temporal and larger spatial mean. Testing this possibility is left for the future work.

### c. Global distribution of retrieval error in temperature and humidity profiles

The statistical evaluation of the agreement between the true and the retrieved changes discussed in the previous section demonstrates that our retrieval error for changes in temperature and specific humidity profiles is relatively small above the boundary layer. Although the agreement between the true and the retrieved changes degrades in the boundary layer, the distribution of errors is nearly symmetric around 0 (Fig. 9). Therefore, when a trend of temperature or specific humidity changes is derived based on many retrievals, the bias in the trend is expected to be small.

In addition to the statistical evaluation of retrieval errors from all grid boxes, we also examine the spatial distribution of retrieval errors. Figures 10 and 11 show the spatial distribution of mean absolute retrieval errors for, respectively, temperature and relative specific humidity changes in the upper, middle, and lower troposphere. We compute the standard deviation of 46 retrieval errors for each 10° × 10° grid box. We then exclude retrieved errors outside 5 times the standard deviation in plotting Figs. 10 and 11. The true temperature change increases from tropics to the polar regions throughout the troposphere (the bottom plot of Fig. 10). The magnitude of the temperature change in higher latitudes is greater in the Northern Hemisphere than in the Southern Hemisphere, especially in the lower troposphere, because of a larger land area. As a consequence, the retrieval error is small over the tropical ocean while a larger error occurs in higher latitudes. In addition, a larger error occurs over land in the midlatitude. The number of regions with the mean absolute error in the temperature change less than the corresponding absolute mean true temperature change for the upper, middle, and lower troposphere is, respectively, 47%, 60%, and 56% of the total number of regions of 648.

While the spatial pattern of relative specific humidity change is less dependent on latitude in the upper troposphere, the change in the lower troposphere increases from the tropics to higher latitudes. Similar to the spatial pattern of the error in the temperature change, a large retrieval error tends to occur over regions where large true changes are present. This implies the error tends to be larger when the residual term caused by the linearization is larger. Figure 11 shows, however, that the error in retrieved relative specific humidity change in the boundary layer over midlatitudes and polar regions is significantly larger. The regional absolute mean retrieval error is within the corresponding regional absolute mean true relative specific humidity change for 69% and 56% of 648 grid boxes at 346 and 547 hPa, while the regional mean retrieval error for most grid boxes exceeds the corresponding regional mean true change at 825 hPa. This in part is due to a few large errors occurring at 825 hPa for most grid boxes. Presence of large errors in retrieving the relative specific humidity change is apparent in the top-right plot of Fig. 7.

## 5. Discussion and summary

We developed a linear inversion algorithm to derive changes of surface skin temperature and atmospheric temperature profile and relative changes of specific humidity profile using temporally and spatially averaged all-sky spectral radiances. The algorithm is similar to the algorithm used in Kato et al. (2014). Two noteworthy differences from the earlier algorithm are inclusion of residual and cloud spectral kernels in the eigenvector of principal component form. Six eigenvectors are used for the residual spectral kernels. An underlying assumption is that the spectral shape of the principal components is constant and their magnitude varies temporally and spatially. Similarly, three eigenvectors of principal components are used for the cloud spectral kernels. Although the direct link between physical cloud properties and principal components is not established, inclusion of cloud spectral kernels is needed to derive temperature and relative specific humidity profile changes using all-sky spectral radiance differences. As Fig. 3 indicates, the spectral shape of cloud radiative effect is relatively constant, although the spectral shape of cloud fraction, optical thickness, and cloud-top height might depend on background conditions. The use of principal components mitigates variability of the spectral shape of individual cloud kernels.

We investigated the error in retrieved skin temperature, temperature and specific humidity profiles using synthetic nadir-view spectral radiances averaged over 16 days and over a 10° × 10° grid box. The RMS difference of retrieved and true skin temperature differences is 0.59 K. The median of absolute errors in the temperature change is less than 0.5 K except for the lowest three layers. The median of absolute error in the relative specific humidity changes is less than 10% above 825 hPa. Although the error is larger in the boundary layer, the distribution of the errors is nearly symmetric around 0. Because a large portion of spectrum used for the retrieval is sensitive to the upper troposphere, the information content of lower-troposphere variables in the mean spectral difference might be small. Therefore, selecting spectral ranges that are sensitive to the lower atmosphere to focus on the lower troposphere or the use of mean clear-sky spectral radiances might improve the retrieval results for the lower troposphere.

While spectral radiances averaged over 16 days and over a 10° × 10° grid box are used in this study, the temporal averaging length and spatial averaging scale are not fixed and can be larger. The algorithm developed in this study can be used for spectral radiances averaged over longer temporal and larger spatial scales. As the temporal averaging length and spatial averaging scale increase, the magnitude of the residual term is expected to decrease and all changes are linearly combined because variability due to clouds decreases with increasing spatial scale faster than the variability due to temperature and humidity (Brindley et al. 2015). Therefore, we expect that the retrieval error decreases as increasing spatial (e.g., 10° zonal) and temporal (e.g., monthly, seasonal, or annual) scales. However, to understand how the retrieval error changes with temporal averaging length and spatial averaging scale is left for a future study. In addition, the next step is to apply this algorithm to observed spectral radiance differences. In doing so, issues that are ignored in this study such as the uncertainty in modeling spectral radiances and uncertainty in surface spectral emissivities need to be addressed. When actual observations are used, the validation of temperature and humidity changes over a large area and a long time period needs to be done using multi-instruments and multidata products. All observed changes, temperature, humidity, radiation budget, clouds, and aerosols need to be consistent once they are converted to irradiance changes at the top of atmosphere. That is the sum of all contributions to the top-of-atmosphere broadband irradiance change needs to agree with observed broadband irradiance change.

The application of the algorithm is not limited to the retrieval of temperature and relative specific humidity changes. For example, the algorithm can be used to correct bias errors in temperature and specific humidity vertical profiles using the difference of computed and observed spectral radiances. When the algorithm discussed in this paper is used for the modeled and observed spectral radiance difference, inverted temperature and specific humidity differences are the bias error in temperature and specific humidity profiles that are used for modeling the spectral radiance. The algorithm used in this way helps in producing climate data records constrained by observations. For example, Kato et al. (2018) use difference of AIRS derived upper tropospheric temperature and specific humidity and those from reanalysis product in order to reduce the bias error in computed surface and atmospheric irradiances. The bias correction is applied to monthly mean gridded temperature and specific humidity profiles. The inversion algorithm developed in this study can be used to correct the bias error in temperature and specific humidity profiles from reanalysis, provided that the error in spectral emissivity and gaseous absorption cross section used for spectral radiance computations is either small or known.

## Acknowledgments

We thank Dr. Oleg Dubovik of Laboratoire d’Optique Atmosphérique and Dr. Stephen Leroy of Atmospheric and Environmental Research Inc. for discussions and suggestions at the earlier stage of the work and two anonymous reviewers for providing constructive comments that improved the manuscript. This research was performed when Fang Pan was supported by an appointment to the NASA Postdoctoral Program (Fang Pan) at Langley Research Center, administered by Universities Space Research Association through a contract with NASA. Seiji Kato, Fred Rose, and Alexander Radkevich were supported by the NASA CERES project.

## REFERENCES

*J. Quant. Spectrosc. Radiat. Transfer*

*J. Climate*

*Bull. Amer. Meteor. Soc.*

*Atmos. Meas. Tech.*

*J. Geophys. Res.*

*Bull. Amer. Meteor. Soc.*

*J. Geophys. Res.*

*J. Climate*

*Atmos. Meas. Tech.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*J. Climate*

*J. Climate*

*J. Climate*

*J. Climate*

*J. Climate*

*Appl. Opt.*

*IEEE Trans. Geosci. Remote Sens.*

*Geophys. Res. Lett.*

*J. Climate*

*J. Climate*

*J. Geophys. Res.*

*J. Geophys. Res. Atmos.*

*Remote Sens.*

*J. Appl. Meteor. Climatol.*

*J. Geophys. Res.*

*IEEE Trans. Geosci. Remote Sens.*

*J. Appl. Remote Sens.*

*J. Climate*

*J. Geophys. Res. Atmos.*

*Bull. Amer. Meteor. Soc.*

*CALIPSO*mission

*Bull. Amer. Meteor. Soc.*

*IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.*

## Footnotes

^{a}

Current affiliation: JD Digits, Beijing, China.

Denotes content that is immediately available upon publication as open access.

Photopolarimetry in Remote Sensing, G. Videen, Y. Yatskiv, and M. Mishchenko, Springer, 65–106