Retrieval of Vertical Profile of Cirrus Cloud Effective Particle Size Using Reflected Line Spectra in 1.38 μm Band

This paper presents a new retrieval method for inferring the vertical profile of cirrus cloud effective particle size by using solar reflected line spectra in the 1.38‐μm band. The retrieval method is based on the maximum‐photon penetration principle coupled with the constrained linear inversion. This approach takes advantage of the vertical stratification of cirrus cloud effective particle size as well as absorption lines of water vapor of different intensity, which contain rich information on the vertical structure of cloud particle size. Reflected radiances at different wavenumbers provide the effective‐size information at different heights within cirrus associated with photon different penetration depths. Assuming a vertical profile of effective size monotonically decreasing toward cloud top and using results based on “exact” radiative transfer computations, we perform retrieval of the effective size for a number of model cirrus to check for algorithm accuracy. The retrieved profile of effective size is close to the model profile for cirrus optical depth less than about eight with an uncertainty range of 2.2–4.2 μm. In addition, we further carry out a sensitivity study involving the retrieved effective size in connection with different water vapor profiles and demonstrate that the difference from the model is only several percent except for the cloud top if an appropriate wavenumber set is selected. The results from this study suggest that the present method can be applied to realistic remote sensing of the vertical profile of cirrus cloud particle size.


Introduction
Cirrus cloud particle size and ice water path are fundamental parameters that determine the relative strength of the solar albedo (reflecting sunlight) and infrared greenhouse (trapping thermal infrared radiation) effects, essential for the discussion of cirrus and climate (Liou & Yang, 2016). Satellite remote sensing of cirrus optical depth (or ice water path) and effective particle size D e has been performed by many researchers (e.g., Jiang et al., 2019;Ou et al., 2013;Zhao et al., 2018). However, the inferred cloud effective particle size under the conventional assumption of homogeneous cloud is underestimated for cirrus clouds observed by the Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance/radiance datasets (Wang et al., 2009). Balloon-borne ice crystal replicator observations by Miloshevich and Heymsfield (1997) and Heymsfield and Iaquinta (2000) have shown a decrease of cloud particle size from cloud base to cloud top. Observations on water clouds show an opposite vertical stratification on cloud particle size, i.e., a monotonic increase of cloud particle size from cloud base to cloud top. If the vertical inhomogeneity is neglected, i.e., the conventional assumption of homogeneous cloud is applied, retrieved cloud particle sizes are overestimated (King & Vaughan, 2012) due to the stratification opposite to the cirrus cloud.
In this study, we utilize different penetration depth of photons within a cloud at different wavenumbers in the 1.38-μm band, which contains numerous water vapor absorption lines of various intensities and vertical stratification of cirrus cloud particle sizes, to infer vertical size profiles in cirrus clouds. For this purpose, we use the maximum-photon penetration principle to determine weighting functions, i.e., relative contributions from different layers within the cloud to the bidirectional reflectance at the top of the atmosphere. Reflectance measurements can be obtained from satellites such as Orbiting Carbon Observatory 2 (OCO-2) (Liou & Yang, 2016) and MODIS. We can then perform retrieval of the vertical profile of cirrus cloud effective particle size using the constrained linear inversion method (Liou, 2002) to solve the Fredholm equation of the first kind. Using the combined maximum-photon penetration principle and constrained linear inversion approach, we can resolve the vertical inhomogeneity of cirrus cloud particle size and therefore obtain more accurate cloud effective particle size averaged over the vertical direction than the conventional method. The organization of this paper is as follows: in section 2, we present cirrus ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

RESEARCH ARTICLE
10.1029/2020EA001119 Y. Gu and Y. Takano contribute equally to this work Key Points: • Development of a retrieval method for the vertical profile of cirrus cloud particle size using reflected radiances in the 1.38-μm band • Maximum-photon penetration principle and constrained linear inversion approaches are combined to determine ice crystal in cirrus clouds • The present approach can be applied to water vapor profiles with optical depth less than 8 by selecting a number of wavenumbers cloud models and methods to invert bidirectional reflectance to a vertical profile of cirrus cloud particle size. Reflectances from model cirrus clouds are used to validate the proposed method for different water vapor profiles in section 3. Conclusions are given in section 4.

Models and Methods
Single-scattering properties of cirrus cloud particles are computed by the geometric ray tracing method (e.g., Takano & Liou, 1989a) assuming polydisperse solid columnar hexagonal ice crystals. Updated complex refractive index of ice complied by Warren and Brandt (2008) is used for the computation. We use 8 size distributions, Cold Ci (23.9 μm), T = −60°C Ci (30.4 μm), Cs (41.5 μm), T = −20°C Ci (57.9 μm), Nov 1 Ci (75.1 μm), Nov 2 Ci (93.0 μm), Oct 25 Ci (110.4 μm), and Ci Uncinus (123.6 μm), where numbers in parenthesis denote cloud effective particle size (Fu & Liou, 1993). We interpolate cloud particle effective size of inhomogeneous cirrus using bidirectional reflectances for these eight size distributions. For application to the vertical structure of cirrus, a monotonic decrease of cloud effective particle size D e from the cloud base upward shown in Figure S1a is used as a first approximation, mimicking the balloon-borne observations by Miloshevich and Heymsfield (1997) and Heymsfield and Iaquinta (2000). The optical depth of each subcloud layer Δτ c is given by τ c /8. Extinction coefficients for eight subcloud layers are expressed in terms of β ek = τ c /(8Δz k ), where τ c is the cloud optical depth, and Δz k is the geometric thickness of subcloud layer which is given by the average vertical distance of a diamond and the next diamonds above and below, as shown in Figure S1a. Figure S1b shows the discretized effective size as a function of relative cloud optical depth τ/τ c , which is used in the retrieval. Figure S2a shows the bidirectional reflectance R i ν , serving as the measurement input for a vertically inhomogeneous, plane-parallel cirrus cloudy atmosphere defining the cloud structure illustrated in Figure S1b with a cloud optical depth τ c of 1. They are computed by the "exact" doubling-adding method for radiative transfer (Takano & Liou, 1989b) coupled with the correlated k-distribution method for water vapor absorption (Liou, 2002). Throughout the paper, the cirrus cloud top is placed at z t = 9 km with a thickness Δz c = z tz b of 1 km embedded in a midlatitude summer atmosphere (McClatchey et al., 1972). The solar and emergent zenith angles θ 0 and θ are, respectively, 60°and 0°, while the relative azimuth angle Δϕ = ϕ -ϕ 0 is 0°. The surface albedo α s is assumed to be 0.1. These conditions are used throughout the paper, unless stated otherwise. Using the trial-and-error method in order for magnitudes of bidirectional reflectance to scatter over a wide range in Figure S2a, 16 wavenumbers are selected with a line resolution of 1 cm −1 from 7100-7400 cm −1 where the non-zero surface reflection is neglected for the retrieval of D e due to strong water vapor absorption, as shown in Figure 5.34a of Liou and Yang (2016). Figure S2b displays cloud effective particle size D * e ν ð Þ corresponding to R i ν in Figure S2a, with interpolation using a lookup table consisting of ( R h νk , D ek ), where k = 1, …, 8, R h νk denote the nadir bidirectional reflectances corresponding to homogeneous cloudy atmospheres with eight effective sizes D ek . Deviations of D * e ν ð Þ from a constant value are associated with the vertical inhomogeneity of cirrus cloud particle sizes. The wide range of variations in D * e ν ð Þillustrates that the results displayed in Figure S2(b) whose spectral range and wavelengths selected contain rich information on cirrus particle size vertical profile. The interpolated D * e ν ð Þ simulating the effective size retrieved using a conventional method can be expressed as: where D e (τ) is a vertical effective size profile shown in Figure S1(b), the term w ν (τ, τ c )represents a weighting function, i.e., relative contributions from different layers within clouds to the bidirectional reflectance at the top of the atmosphere. Equation 1 is the well-known Fredholm equation of the first kind in which w ν (τ, τ c ) is a kernel such that D e (τ) is a function to be inverted from a set of D * e ν ð Þ . Using the maximum-photon penetration principle, the weighting function can be written as follows (2) when a differential layer Δτ is added to the base of a sublayer of the cloud.
Þrepresents only reflected photons that penetrate to a maximum depth between τ and τ + Δτ. Therefore w v (τ, τ c ) is normalized fraction of reflected photons that reach the bottom differential layer Δτ. Further, we discuss the constrained linear inversion for the retrieval of vertical profile of effective size D e (τ) such that equation 1 is discretized as follows: where j (=1, …, M) is the wavenumber index, k (=1, …, N) denote the index of cloud optical depth, and A jk = w jk Δτ k . Alternatively, using matrix-vector form, equation 3 can be expressed by where A −1 is the inverse of a matrix A. A direct inversion of equation 4 may produce negative values of D e , i.e., the solution is unstable. For this reason, we adopt the constrained linear inversion method (Liou, 2002): where γ is an arbitrary smoothing coefficient and A T denotes the transpose of A. We select H in equation 5 such that the sum of squares of the first differences of the solution is minimized as follows: The coefficient γ in equation 5 is used to determine χ 2 defined in the form: where D comp e À Á k is the effective sizes computed from equation 5, and D model e À Á k are the model effective sizes defined in Figure S1(b).  Figure S2a. When water vapor absorption is weaker, weighting functions, e.g., at 7171.5 cm −1 (orange ○), increase into deeper cloud layers for τ~0.3, induced by a larger single scattering albedo due to multiple scattering effects. When water vapor absorption is strong, weighting functions, e.g., at 7226.5 cm −1 (red ×) increase into shallower cloud parts (down to τ~0.2) due to a smaller single scattering albedo. When water vapor absorption is weak, we see more information at lower cloud levels, because bidirectional reflectance R i ν approaches to a constant value slower when cloud optical depth τ becomes larger. In the case of strong water vapor absorption, the reverse is shown. Figure 1(b) displayed weighting functions for τ c = 4. In comparison to τ c = 1, relative intervals between weighting functions at different wavenumbers are narrower especially at lower cloud levels. This is partly because of the effect of cloud scattering which has little dependence on wavenumber related to the narrow range of wavenumber dominates the effect of water vapor absorption. Also, when τ c is larger,

Computed Results and Discussion
Þ in equation 2 becomes smaller associated with the principles of invariance (Liou, 2002). For this reason, the relative intervals between weighting functions of different wavenumbers becomes smaller for τ c = 4.  Figure S6 shows the retrieved effective size for τ c = 8 as a function of optical depth (a) and height (b). Table 1 lists cloud effective particle sizes D e averaged over the vertical direction as a function of the cloud optical depth along with the model average. Differences from the model are within~1 μm and the relative error is less than 2%. We note that agreement is excellent with respect to the average effective size D e . Although effective sizes differ among 1.6, 2.2, and 3.7 μm MODIS channels, the present method has no such difference, even in the case for vertically inhomogeneous clouds.  We perform sensitivity experiments in order to investigate whether the proposed method can be applied to various water vapor profiles. Figure 4 shows the retrieved D e for different cloud heights but keeping other conditions unchanged. Retrieved values D e for cloud top heights z t of 10, 11, and 12 km show negligible differences. Figures S7 -S10 show the corresponding weighting functions for z t = 8, 10, 11, and 12 km. To evaluate R i ν τ c ð Þ in the denominator of equation 2, we assume T ¼ T b þ T t ð Þ =2 for representative temperatures of the entire cloud, where T b and T t are, respectively, cloud-base temperature and cloud-top temperature. T is larger than T t . As the temperature lapse rate is nearly constant at 7 ≤ z ≤ 12 km, the cloud temperature T becomes T -ΔT, when the cloud height is raised by a certain height. Thus, the initial T and T t becomes T -ΔT and T t -ΔT, respectively. On the basis of Clausius-Clapeyron equation (see, e.g., Figure 1.12 of Liou & Yang, 2016), the saturation water vapor pressure varies exponentially. Hence, we can obtain the following relation where ρ s denotes the saturation water vapor density over ice. Further, consider an upper level of the cloud (τ/ τ c < 0.5) corresponding to T ≈T t . Using equation 8, we compare variation of R i ν τ c ð Þ in the denominator with that of R i ν τ ð Þ in the numerator of equation 2, when the cloud height is raised and the cloud temperature is reduced by ΔT. R i ν τ c ð Þincreases more than R i ν τ ð Þ due to more reduction of water vapor amount, as expressed by equation 8. Thus the peak of weighting function becomes smaller when cloud moves higher, as shown in Figures S7 -S10, except for ν = 7117.5 cm −1 where water vapor absorption is exceptionally strong. Further, we consider a lower level of the cloud (τ/τ c > 0.5) corresponding to T ≈T b . R i ν τ c ð Þin the denominator is com- (2) becomes larger due to weaker water vapor absorption and the weighting function also becomes larger with an increase of cloud height, as shown in Figures S7 -S10. In the case of z t = 8 km and τ c = 4 in Figure 4(b) whose water vapor amount is the most abundant due to the lower height, the retrieved D e deviates from the model D e more than that in Figure 4(a) whose τ c is 1. This is because the effect of water vapor absorption is enhanced by multiple scattering. Figure 4 suggests that the deviation of retrieved D e from the model is related to water vapor amount. By selecting another set of wavenumbers where water vapor absorption is weaker, as shown in Figure 5, we perform another retrieval which shows better agreement between retrieved and model D e (Figure 6). In contrast to the paragraph [33] of King and Vaughan (2012) who take above cloud total column water vapor amount into account as the state vector, above cloud water vapor absorption does not affect retrieved cloud particle size D e , as mentioned in Text S1.   Figure 7 illustrates the retrieved D e for different cloud thickness Δz c (= 0.5, 1, 2, and 3 km) but keeping the other conditions unchanged. When the cloud top height is reduced, the cloud base height is raised. Even though the cloud top temperature is changed, the cloud middle temperature and approximate average Figure 6. Retrieved and model cirrus cloud effective particle size for the case of z t = 8 km but keeping the same cloud thickness Δz c = 1 km as a function of cloud optical depth for (a) τ c = 1 and (b) τ c = 4. But we use another wavenumber set shown in Figure 5 in addition to the wavenumber set shown in Figure S7.  Figure 7 shows the effect of water vapor amount with the same temperature, whereas Figure 4 shows the effect of water vapor amount but with different temperature. When the cloud thickness Δz c increases, the retrieved D e at an upper part of the cloud (τ/τ c < 0.5) becomes smaller, while the retrieved D e at a lower level of the cloud (τ/τ c > 0.5) becomes larger. Figures   Figure S11. S11 -S13 show the corresponding weighting functions for Δz c = 0.5, 2, and 3 km. When the cloud thickness Δz c increases, the water vapor amount in the cloud also increases so that R i ν τ c ð Þ in the denominator of equation 2 decreases. When the cloud thickness Δz c increases, the water vapor amount at an upper part of the cloud decreases due to higher cloud height so that R i ν τ ð Þin the numerator of equation 2 increases. As a result, the weighting function at an upper level of the cloud (τ/τ c < 0.5) becomes much larger, as shown in Figures S11 -S13. We then consider a lower level of cloud (τ/τ c > 0.5). R i ν τ c ð Þ is comparable to R i ν τ ð Þ in equation (2) associated with τ≈ τ c . But the water vapor amount in a differential layer Δτ increases with an increase of cloud thickness Δz c , so that R i ν τ þ Δτ ð Þ becomes smaller related to the stronger water vapor absorption. Therefore, the weighting function also decreases with an increase of cloud thickness Δz c , as shown in Figures S11 -S13. To ameliorate the poor agreement between the retrieved and model D e in Figure 7(a), we select another set of wavenumbers, as shown in Figure 8. Intervals of weighting functions of different wavenumber in Figure 8 are wider than those in Figure S11 except for ν = 7117.5 cm −1 . For this reason, the weighting functions in Figure 8 are more orthogonal than those in Figure S11. Consequently, we find better agreement in Figure 9 using the new set of wavenumbers. The three sets of wavenumber selection are listed in Table S1. Figure 10 shows the retrieved D e for different atmospheric profiles but keeping other conditions the same as in Figure 2. TRP, MLS, SAS, MSW, and SAW denote, respectively, tropical, midlatitude summer, subarctic  summer, midlatitude winter, and subarctic winter model atmospheres (McClatchey et al., 1972). TRP is close to MLS; SAS is close to MLW. TRP, MLS, SAS, MLW, and SAW in Figure 10 correspond, respectively, to zt = 8, 9, 10, 11, and 12 km in Figure 4, since the cloud temperature decreases monotonically. Figures S14 -S17 show the corresponding weighting functions for TRP, SAS, MLW, and SAW model atmospheres. Also, Figures S14 -S17 are close to Figures S7 -S10, such that weighting function at an upper (lower) part of the cloud decreases (increases) in order of TRP, SAS, MLW, and SAW.

Conclusions
In this study, we find the monotonic relation between effective size D ek and bidirectional reflectance R h νk , such that photons of different wavenumbers penetrate different depths of stratified cirrus cloud in terms of cloud particle size associated with differences in water vapor absorption. Using these characters, we retrieve the vertical profiles of cirrus cloud effective particle size. Specifically, we have developed a retrieval method for the vertical profile of cirrus cloud particle size using the maximum-photon penetration principle coupled with the constraint linear inversion approach to minimize sum of square of the first differences of the solution using reflected line spectra in the 1.38-μm band. We have performed the retrieval of effective cloud particle size as functions of cirrus cloud optical depth and height for optical depth of less than or equal to 8. Except for the cases involving the cloud top and bottom, the retrieved effective size agrees well with the "exact" model effective size. In addition, sensitivity studies for the retrieved size on cloud height, cloud thickness and atmospheric profile have been carried out to check the validity of present algorithm for various water vapor profiles. Average root mean square error for cloud effective particle size is 4.2 and 2.2 μm, respectively, for τ c = 1 and 4. Cirrus cloud optical depths are usually less than 6 (e.g., Kox et al., 2014). Thus, the present method can be applied to realistic remote sensing of the vertical profile of cirrus cloud particle size using the 0.76-μm O 2 A-band hyperspectral reflectance observed by the Orbiting Carbon Observatory 2 (OCO-2) (Liou & Yang, 2016), although there is still room to select better sets of wavenumbers and another retrieval method, e.g., the nonlinear iterative inversion.