Retrieval of terahertz ice cloud properties from airborne measurements based on the irregularly shaped Voronoi ice scattering models

Currently, terahertz remote sensing technology is one of the best ways to detect the microphysical properties of ice clouds. Influenced by the representativeness of the ice crystal scattering (ICS) model, the existing terahertz ice cloud remote sensing inversion algorithms still have significant uncertainties. In this study, based on the Voronoi ICS model, we developed a terahertz remote sensing inversion algorithm of the ice water path (IWP) and median mass diameter (D_me) of ice clouds. This study utilized the single-scattering properties (extinction efficiency, single-scattering albedo, and asymmetry factor) of the Voronoi, sphere, and hexagonal column ICS models in the terahertz region. Combined with 14 408 groups of particle size distributions obtained from aircraft-based measurements, we developed the Voronoi, sphere, and column ICS schemes based on the Voronoi, sphere, and column ICS models. The three schemes were applied to the radiative transfer model to carry out the sensitivity analysis of the top-of-cloud (TOC) terahertz brightness temperature differences between cloudy and clear skies (BTDs) on the IWP and D_me. The sensitivity results showed that the TOC BTDs between 640 and 874 GHz are functions of the IWP, and the TOC BTDs of 380, 640, and 874 GHz are functions of the D_me. The Voronoi ICS scheme possesses stronger sensitivity to the D_me than the sphere and column ICS schemes. Based on the sensitivity results, we built a multi-channel look-up table for BTDs. The IWP and D_me were searched from the look-up table using an optimal estimation algorithm. We used 2000 BTD test data randomly generated by the RSTAR model to assess the algorithm's accuracy. Test results showed that the correlation coefficients of the retrieved IWP and D_me reached 0.99 and 0.98, respectively. As an application, we used the inversion algorithm to retrieve the ice cloud IWP and D_me based on the Compact Scanning Submillimeter-wave Imaging Radiometer (CoSSIR) airborne terahertz radiation measurements. Validation against the retrievals of the Bayesian algorithm reveals that the Voronoi ICS model performs better than the sphere and hexagonal column ICS models, with enhancement of the mean absolute errors of 5.0 % and 12.8 % for IWP and D_me, respectively. In summary, the results of this study confirmed the practicality and effectiveness of the Voronoi ICS model in the terahertz remote sensing inversion of ice cloud microphysical properties.

1 Introduction

Ice clouds account for about 20 %–30 % of the total global cloud mass (Liou, 1992; Rossow and Schiffer, 1991). They play an essential regulatory role in the global radiation balance and the water cycle (Hong et al., 2016; Stephens et al., 2012). Microphysical properties such as the ice water content, ice particle size, orientation, and shape are the main influencing factors of the scattering and radiative properties of ice clouds (Li et al., 2022; Yi et al., 2017; Zhao et al., 2018; Chen and Zhang, 2018) and, in turn, affect the radiation budget (Heymsfield et al., 2013, 2017; Rossow, 2014). The latest report of the 6th Intergovernmental Panel on Climate Change (IPCC6) (Forster et al., 2021) identifies cloud radiative properties and their feedback effects as the largest source of uncertainty in the overall climate feedback, with errors in ice clouds being one of the most significant factors (Zhang et al., 2021). Therefore, the accurate acquisition of microphysical properties of ice clouds is of great importance for studying global climate change and improving the accuracy of numerical model simulations (Baran, 2009, 2012). Remote sensing techniques are one of the most effective means of obtaining microphysical and radiative properties of ice clouds, mainly including ground-based (Cimini et al., 2006), aircraft-based (Fox et al., 2017), and satellite-based remote sensing observation technologies (Yang et al., 2015, 2018). Currently, large numbers of passive sensors (visible, infrared, and microwave detectors) have been developed, and related ice cloud retrieval algorithms have been reported in a substantial body of literature (Nakajima and King, 1990; Nakajima et al., 1991, 2019; Nakajima and Nakajima, 1995; Platnick et al., 2003, 2017; Fox et al., 2019; Brath et al., 2018). The visible and infrared wavelengths are sensitive to the visible optical depth and cloud top (Minnis et al., 1993a, b). The millimetre-wave ice cloud remote sensing technique is more suited to detect vertical cloud properties. Sensors such as the Millimeter-wave Imaging Radiometer (MIR) (Racette et al., 1992) and Special Sensor Microwave Water Vapor Sounder (SSM/T-2) have been used in several studies of ice water path (IWP) and particle size retrievals (Lin and Rossow, 1996; Liu and Curry, 1998, 1999). MIR channels at 89, 150, and 220 GHz have been used by Deeter and Evans (2000) and Liu and Curry (2000) to retrieve IWP and particle size in cirrus anvils over the tropical ocean. Compared to passive sensors, the cloud satellite radar (CloudSat), with an onboard millimetre-wavelength (94.05 GHz) radar, and the radar–lidar cloud product (DARDAR) (Ceccaldi et al., 2013) present new opportunities to infer the microphysical properties of ice clouds on a global scale. In practice, the current detectors and approaches are confined to a limited range of ice particle sizes (Cho et al., 2015). For example, visible and infrared sensors are only sensitive to particles smaller than 50 µm (Evans et al., 2005). Additionally, microwave detectors are mainly useful for particles larger than 500 µm (Fox, 2020; Fox et al., 2019). There is a pressing need to develop an effective frequency region for obtaining comprehensive information about ice particles. To bridge the gap between the technologies of visible–infrared and microwave measurement of ice clouds, the terahertz (THz) wavelength between the infrared and microwave regions has the potential advantage of complementing existing visible–infrared and microwave techniques (Gasiewski, 1992).

The terahertz wave is the sub-millimetre wave with frequencies in the range of 0.1–10 THz and corresponding wavelengths of 0.03–3 mm, comparable to the size of ice particles in ice clouds. Many theoretical studies (Evans and Stephens, 1995a, b; Evans et al., 1998) have shown that the passive terahertz wave has a higher detection capability and sensitivity to the IWP and particle size of ice clouds (Liu et al., 2020). Although terahertz waves were proposed for ice cloud remote sensing in the 1960s (Gao et al., 2016), the technology of terahertz radiometry of ice clouds lagged behind the theory (Evans et al., 2005). It is only within the last decade that terahertz radiometry has been applied to aircraft-based and satellite-based ice cloud remote sensing. With advances in terahertz detectors, researchers have successively developed aircraft-based terahertz radiometers (Gao et al., 2016), such as the Sub-millimeter Wave Cloud Ice Radiometer (SWCIR) (Evans et al., 2002), the Compact Scanning Submillimeter-wave Imaging Radiometer (CoSSIR) (Evans et al., 2012), and the International SubMillimeter Airborne Radiometer (ISMAR) (Fox et al., 2017). Also, research institutions have developed satellite-based terahertz radiometers, including the Superconducting Submillimeter-Wave Limb Emission Sounder (SMILES) (Inatani et al., 2000), Ice Cloud Imager (ICI) (Eriksson et al., 2020; Kangas et al., 2014), and IceCube (Wu et al., 2014), which were proposed in 2013 and are still under development.

Several methods have been reported to retrieve the ice clouds' IWP and particle size from terahertz brightness temperature. Inversion methods can be divided into physical and linear regression methods (Weng et al., 2019). The linear regression method is efficient and convenient, but it lacks a definite physical mechanism and is heavily dependent on the accuracy of the pre-calculated database. The physical method includes the Bayesian algorithm (Evans et al., 2005, 2002), the look-up table (LUT) method (Li et al., 2017, 2018; Liu et al., 2021), and the neural network algorithm (Jimenez et al., 2003). Evans et al. (1998) modelled terahertz brightness temperature using a polarized radiative transfer model based on eight ice particle shapes calculated by the discrete dipole approximation (DDA) method. The study found that the ice particle shape plays a vital role in modelling ice cloud scattering in the terahertz region. Evans et al. (2002) developed a Monte Carlo Bayesian integration (MCBI) algorithm to retrieve ice clouds' IWP and median mass diameter (D_me) from simulated SWCIR brightness temperatures. According to the validation results, the total median retrieval error for clouds with IWP of more than 5 g m⁻² is roughly 30 % for IWP and 15 % for D_me (Evans et al., 2002). Moreover, Evans et al. (2005) utilized the MCBI method to retrieve the IWP and D_me of ice clouds from the CoSSIR brightness temperatures (referred to as the CoSSIR-MCBI hereafter). During the retrieval procedure, one of five particle shapes (spherical snow, aggregates of frozen droplets, aggregates of hexagonal plates, and aggregates of plates and hexagonal columns) developed by the DDA method (Evans and Stephens, 1995a) was selected for each ice cloud retrieval. The CoSSIR-MCBI results are validated by the Cloud Radar System data and showed a good agreement of radar backscattering with errors smaller than 5 dB. Later, Jimenez et al. (2007) used the neural network method combined with the radiative transfer code to retrieve the IWP and D_me of ice clouds. The ice cloud microphysical input is based on one of three randomly oriented particle shapes (solid columns, hexagonal plates, and four-bullet rosettes) simulated with the DDA method (Evans and Stephens, 1995a; Evans et al., 1998) in the microwave region. Results showed overall median relative errors of around 20 % for IWP and 33 µm for D_me for a mid-latitude winter scenario and 17 % for IWP and 30 µm for D_me for a tropical scenario. Based on early studies and the background, Buehler et al. (2007) proposed a formal scientific requirement for a passive sub-millimetre-wave cloud ice mission. The requirements are that the low IWP should be less than 10 g m⁻², the high IWP should be less than 50 %, and the particle diameter should be less than 50 µm. Li et al. (2016) investigated the effects of five ice crystal scattering (ICS) models from the single-scattering property database of Hong et al. (2009) at frequencies ranging from 100 to 1000 GHz, namely aggregates, hollow columns, flat plates, rosettes, and spheres (Van de Hulst, 1957), on the transmission characteristics of terahertz radiation. Results showed that the ice particle shape is one of the dominant factors affecting terahertz radiation. Lately, Fox (2020) evaluated seven ice particle habits from the ARTS scattering database (Eriksson et al., 2018). They showed that the randomly oriented large column aggregate performs best when simulating brightness temperatures between 183 and 664 GHz.

In summary, the physical method is based on the radiative transfer principle and relies on the forward physical model simulation and effective ice crystal scattering (ICS) model. Different assumptions of ice cloud microphysical properties (shape, size, orientation, and particle size distribution of ice particles) in the forward physical model significantly affect the retrieval of the IWP and particle size of ice clouds. Therefore, a practical and representative ICS model is essential for the ice cloud remote sensing in the terahertz region.

Many aircraft observations have demonstrated that ice clouds comprise a complex and diverse range of non-spherical ice particles (Lawson et al., 2006, 2019; Liou, 1992). It is still challenging for one specific light-scattering method to precisely calculate the single-scattering properties for non-spherical particles with different size parameters (SZPs). The SZP is the ratio of the equivalent-volume sphere's circumference dimension (or π times particle maximum diameters) to the incident wavelength (Baran, 2012; Nakajima et al., 2009). So far, the light-scattering methods can be roughly divided into the approximation method (AM) and the numerical simulation method (NM). The AM method is based on ray-tracing techniques, including the geometrical optics method (GOM) (Macke et al., 1996; Takano and Liou, 1989), which is suitable for large particles. The NM includes the finite-difference time domain (FDTD) (Yang and Liou, 1996a; Yee, 1966) and the DDA (Draine and Flatau, 1994; Yurkin and Hoekstra, 2007) methods, which are appropriate for small particles. Combined with the advantages of the AM and NM methods, several improved geometrical optics approximation (GOA) methods, including the geometric optics integral equation (GOIE) method (Yang and Liou, 1996b), have been developed. Moreover, the boundary element method (Groth et al., 2015; Kleanthous et al., 2022) has been recently applied to complex ice particles. Based on the above-mentioned light-scattering calculation methods, a series of regularly shaped ICS models have been designed, including hexagonal columns, hexagonal plates, bullet rosettes, droxtals, and so on (Yang et al., 2000a, 2013; Fu et al., 1999; Takano and Liou, 1989). Since the regularly shaped ICS models are not fully representative of the scattering properties of natural ice particles, researchers have developed complex ICS models. Yang et al. (2013) developed the surface-roughened non-spherical ICS models and applied them to the production of MODIS C6 ice cloud products (Platnick et al., 2017). C.-Labonnote et al. (2000, 2001) and Doutriaux-Boucher et al. (2000) developed an ICS database for the inhomogeneous hexagonal monocrystal (IHM) containing embedded inclusions (air bubbles and aerosols). The IHM database was applied for the ice cloud retrievals from the French satellite Polarization and Directionality of the Earth's Reflectance (POLDER) measurements (Deschamps et al., 1994). Furthermore, Baran and Labonnote (2007) and Baran et al. (2014a) developed an ensemble ice particle model made of hexagonal column ice particles for use in the Met Office Unified Model Global Atmosphere 5.0 (GA5) configuration (Baran et al., 2014b). Unlike the above-mentioned regularly shaped ICS models, Letu et al. (2012) and Ishimoto et al. (2013) analysed ice particles of different shapes and sizes from many NASA aircraft observations and developed an irregularly shaped complex Voronoi ICS model. The single-scattering property database of the Voronoi ICS model in the visible and infrared regions uses a combined FDTD, GOIE, and GOM approach. The Voronoi ICS model has been adopted for generating official ice cloud products for the Second Generation gLobal Imager (SGLI)/Global Change Observation Mission – Climate (GCOM-C) (Letu et al., 2012, 2016; Nakajima et al., 2019); AHI/Himawari-8 (Letu et al., 2018); and the Multi-Spectral Imager (MSI)/Earth Cloud Aerosol and Radiation Explorer (EarthCARE) satellite programmes (Illingworth et al., 2015), which will be launched in 2023. Furthermore, Li et al. (2022) showed the effectiveness of the Voronoi ICS model in describing ice clouds' global visible and infrared radiative properties in the Community Integrated Earth System Model (CIESM). As a result, researchers have proven the effectiveness of the database of the Voronoi ICS model in the visible and infrared regions in the ice cloud satellite remote sensing and climate model applications (Letu et al., 2016). Recently, Letu et al. (2012) and Ishimoto et al. (2013) have successfully extended the spectral range of the single-scattering property database of the Voronoi ICS model to the terahertz wave region. The database of the Voronoi ICS model in the terahertz region was adopted by Baran et al. (2018) as standard data for the modelling and evaluation of the ISMAR (International SubMillimeter Airborne Radiometer), which the European Space Agency (ESA) and Met Office jointly developed (Kangas et al., 2014). The results showed good evaluation performance of the Voronoi ICS model. However, the effectiveness of the Voronoi ICS model in the terahertz region in a practical retrieval of ice cloud microphysical properties from terahertz radiation has yet to be discovered.

Motivated by the above-mentioned situations, this study aims to apply the Voronoi ICS model to the ice cloud remote sensing retrieval of IWP and particle size from aircraft-based terahertz radiation measurements. To achieve this goal, we use the Voronoi ICS model to create a parameterization scheme (referred to as the Voronoi ICS scheme hereafter) for the ice cloud scattering property in the terahertz region. The Voronoi ICS scheme is employed in the RSTAR radiative transfer model (Nakajima and Tanaka, 1986, 1988) to build the LUT for upward terahertz brightness temperature differences between cloudy and clear skies (BTDs) at the top of the ice cloud (TOC) between multi-channel frequencies. The LUT based on the Voronoi ICS model is compared with that of the sphere and hexagonal column ICS models. Finally, the CoSSIR-MCBI evaluates the retrieval results searched iteratively from the LUTs for the Voronoi, sphere, and hexagonal column models. This paper is organized as follows: Sect. 2 introduces the data and radiative transfer model used in this study. Section 3 describes the ice cloud parameterization scheme, retrieval algorithm, and validation with the CoSSIR-MCBI algorithm. Section 4 presents the results of the retrieved IWP and particle size; a comparison of the Voronoi, sphere, and hexagonal column ICS models; and the validation of the retrieval algorithm. Section 5 presents the conclusions of this study.

2 Data and model

2.1 Single-scattering property database for the Voronoi, sphere, and column ICS models

This study used the single-scattering property database of three ice crystal habits (Voronoi, sphere, and hexagonal column) in the terahertz radiative transfer forward simulation. For the Voronoi ICS model, the single-scattering property database contains 31 ice particle sizes ranging from 0.25 to 9300 µm and covers 20 terahertz channels with frequencies ranging from 10 to 874 GHz (see Table 1), corresponding to wavelengths from 0.03 to 3 cm. The hexagonal column ICS model (referred to as the column ICS model hereafter) is randomly oriented and was defined by Yang et al. (2000b). For the column ICS model, the aspect ratios $a/L$ (defined as the ratio of the semiwidth a of a particle to its length L) are defined as 0.35 and 3.48, respectively, when L is less than 100 µm and greater than or equal to 100 µm. The single-scattering property database of the column ICS model used in the study was developed by Hong (2007) and Hong et al. (2009) using the DDA method. For the column ICS model, the single-scattering property database includes 38 maximum particle dimensions ranging from 2 to 2000 µm and covers 21 terahertz channels with frequencies ranging from 90 to 874 GHz (see Table 1).

Table 1 The frequency channels and maximum particle dimensions of ice particles included in the single-scattering property database of the Voronoi and column ICS models.

The single-scattering properties – including the extinction efficiency, single-scattering albedo, and asymmetry factor of the Voronoi, sphere, and column ICS models in the terahertz region – are used to calculate the terahertz scattering properties of ice clouds. For the Voronoi ICS model, the single-scattering properties are derived from the single-scattering property database developed by Letu et al. (2012, 2016) and Ishimoto et al. (2012) using a combination of FDTD and DDA methods. The DDA method is used to calculate the single-scattering properties of moderate ice particles (SZP>30). The FDTD method is used for small ice particles (SZP<30). As the particle size increases, the shape of the Voronoi ICS model changes and becomes complicated. The geometrical characteristics of Voronoi ICS model shape with increasing ice crystal size have been shown in Fig. 3 and Table 1 in Ishimoto et al. (2012). The mass–dimension and area–dimension power law relationships of the Voronoi ICS model are defined by Ishimoto et al. (2012) and are described in Eqs. (1)–(3) as shown below.

$$\begin{array}{}\text{(1)}& m=\mathrm{0.00528}{D}^{\mathrm{2.1}}\end{array}$$

(in centimetre–gram–second system of units, cgs)

$$\begin{array}{}\text{(2)}& {\displaystyle }{A}_{\mathrm{r}}={\displaystyle \frac{\mathrm{4}G}{\mathit{\pi }{D}^{\mathrm{2}}}}\text{(3)}& {\displaystyle }{A}_{\mathrm{r}}=\mathrm{0.20}{D}^{-\mathrm{0.29}},\end{array}$$

where m is the mass, G is the cross-sectional area, A_r is the area ratio, and D is the maximum particle dimension of the Voronoi ICS model. The single-scattering property database of the sphere ICS model is developed based on the exact solution of the Lorentz–Mie theory (Van de Hulst, 1957). The sphere ICS database contains the same ice particle sizes and terahertz frequencies as the Voronoi ICS database. For the Voronoi and sphere ICS models, calculations of their single-scattering properties utilized the real and imaginary parts of ice from the newest library of the refractive index provided by Warren and Brandt (2008). The refractive indices of the Voronoi and sphere ICS models at frequencies of 10–874 GHz are computed at a temperature of 266 K, according to Warren and Brandt (2008). The refractive indices of the column ICS model at frequencies of 90–340 GHz are derived from Warren (1984) at a temperature of −30 ^∘C.

2.2 Airborne measurements and auxiliary data

The measured terahertz brightness temperature from CoSSIR/ER-2 aircraft during the TC4 mission on 17 and 19 July 2007 is utilized in this study (available at https://espoarchive.nasa.gov/archive/browse/tc4/ER2, last access: 26 December 2021). During the TC4 field campaign, the CoSSIR measured brightness temperatures in channels from 183.3 to 873.6 GHz (183.3 ± 1.0, 3.0, 6.6, 220, 380.2 ± 1.8, 3.3, 6.2, and 640 V and 874 GHz), all with matched beamwidths of about 4^∘ (Evans et al., 2012). According to the studies by Li et al. (2016) and Liu et al. (2021), the atmospheric windows are near 640 and 874 GHz, and the atmospheric absorption peaks are near 380 GHz. The leading absorbing gases are water vapour and ozone. Therefore, both the 640 and 874 GHz brightness temperature are affected by ice clouds, while the brightness temperature of 380 GHz is insensitive to ice cloud microphysical properties. Hence, the 380 GHz minus 640 GHz brightness temperature differences can highlight the brightness temperature depression caused by ice clouds. And the 640 GHz minus 874 GHz brightness temperature differences can reflect the difference in the scattering properties of differently shaped ice clouds. According to Li et al. (2016), the differences between 640 and 874 GHz also can offset the regional errors due to different latitudes and atmospheric profiles. In this study, we choose the centre frequencies of 380, 640, and 874 GHz.

The particle size distributions (PSDs) describe the relationship between ice particle size and particle number concentration and are essential in the calculation of the average scattering properties of ice clouds. In this study, we select 14 408 groups of PSDs from aircraft observation sampling data (available at https://www.ssec.wisc.edu/ice_models/microphysical_data.html, last access: 18 February 2022) (Heymsfield et al., 2013), which are based on 11 field flight observation experiments covering ice cloud observations in mid-latitude and tropical regions. The 11 field programmes span a wide range of locations (ranging from 12^∘ S to 70^∘ N and from 148^∘ W to 130^∘ E) and encompass the temperature range 0 to −86 ^∘C, with altitudes from near the surface to 18.7 km. This dataset is representative of the wide range of conditions where ice clouds are found in the troposphere and lower stratosphere on a near-global scale (Li et al., 2022; Heymsfield et al., 2013). For the fitting of PSDs for the Voronoi, sphere, and column ICS models, we adopted the gamma distribution following Heymsfield et al. (2013) in the form of Eq. (4):

$$\begin{array}{}\text{(4)}& n\left(L\right)=N_{\mathrm{0}}{L}^{\mathit{\mu }}{e}^{-\mathit{\lambda }L},\end{array}$$

where L is the maximum particle dimension, n(L) is the particle concentration per unit volume (e.g. cm⁻³), N₀ is the intercept, λ is the slope, and μ is the dispersion. The physical meaning of the PSDs is that n(L)×dL is the number of particles per unit area.

The total water vapour and ozone column data provided by the ERA5 reanalysis data are used as input data for the radiative transfer model to simulate clear-sky brightness temperature. The ERA5 reanalysis data are the fifth-generation reanalysis data product developed by the European Centre for Medium-Range Weather Forecasts (ECMWF) (Hans et al., 2019). The total water vapour and ozone column data used have a horizontal resolution of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ and a temporal resolution of 1 h. The retrieved results are validated against the IWP and D_me (Evans et al., 2005) from the CoSSIR-MCBI algorithm during the same period. The IWP and D_me are available at https://espoarchive.nasa.gov/archive/browse/tc4/ER2 (last access: 20 May 2022).

2.3 Radiative transfer model

The RSTAR radiative transfer model is a set of numerical radiative transfer models developed by Nakajima and Tanaka (1986, 1988) for the plane-parallel atmosphere. The calculated wavelengths can cover from 0.17 to 1000 µm, and the assumed plane-parallel atmosphere could be divided into 50 layers from sea level to a maximum altitude of 120 km. The RSTAR model contains six atmospheric profiles (tropical, mid-latitude summer, mid-latitude winter, high-latitude summer, high-latitude winter, and US standard atmosphere), including vertical profiles of temperature, pressure, water vapour, and ozone. Calculating gas absorption in RSTAR is based on the K-distribution method developed by Sekiguchi and Nakajima (2008), which considers important atmospheric radiative gases (H₂O, CO₂, O₃, N₂O, CO CH₄, etc.) and trace gases. The K-distribution method parameters were obtained from the HITRAN-2004 database. The RSTAR radiative transfer model assumes the simulated scene is composed of a homogeneous ice cloud layer.

3 Method

Figure 1 shows the general flowchart for the inversion of the IWP and D_me of ice clouds using the CoSSIR brightness temperature measurements. For convenience, the difference between the 640 GHz BTD and the 874 GHz BTD is simplified to BTD_2–3, and the difference between the 380 GHz BTD and the 640 GHz BTD is simplified to BTD_1–2. We refer to the difference between the BTD_1–2 and BTD_2–3 as BTD_1–3. Firstly, based on the single-scattering property database of the Voronoi, sphere, and column ICS models, the atmospheric radiative transfer model RSTAR is used to build LUTs of top-of-cloud (TOC) BTD_2–3 and BTD_1–3 for three ICS models, respectively. With the assumptions of an a priori value of IWP and D_me, the initial TOC BTD_2–3 and BTD_1–3 simulations can be searched from the LUTs for three ICS models. Secondly, we use the RSTAR to construct a clear-sky LUT of TOC brightness temperature (BT) for 380, 640, and 874 GHz, respectively. The inputs are different ozone and water vapour column values under clear-sky conditions (see Sect. 3.2 for details). Then, the ERA5 reanalysis data are used to estimate the clear-sky TOC BT based on the clear-sky LUT. The measured cloudy-sky TOC BTs of 380, 640, and 874 GHz are obtained from the CoSSIR measurements. Consequently, the TOC BTD_2–3 and BTD_1–3 measurements can be calculated using the measured cloudy-sky TOC BT minus the interpolated clear-sky TOC BT. Finally, the cost function will provide an estimate of the coherence between the measured and simulated TOC BTD_2–3 and BTD_1–3. The IWP and D_me are obtained by using the Gaussian–Newtonian non-linear optimization estimation method. The retrieved results are validated against the IWP and D_me retrieved from the CoSSIR-Evans algorithm. In the RSTAR forward simulation, the input parameter of the ice particle size is R_e, which is defined by Eq. (5):

$$\begin{array}{}\text{(5)}& R_{\mathrm{e}}={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{r}^{\mathrm{3}}n\left(r\right)\mathrm{d}r}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{r}^{\mathrm{2}}n\left(r\right)\mathrm{d}r}},\end{array}$$

where r is the equivalent-volume sphere's particle radius, and n(r) is the particle size distribution. According to Sieron et al. (2017), a mass-weighted size D_me would be a more appropriate parameter size than the area-weighted size for describing the PSDs. The D_me is given by Eq. (6):

$$\begin{array}{}\text{(6)}& D_{\text{me}}={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}Dm\left(D\right)n\left(D\right)\mathrm{d}D}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}m\left(D\right)n\left(D\right)\mathrm{d}D}},\end{array}$$

where D is the maximum particle dimension of ice particles, m is the particle mass, and n(D) is the particle size distribution. To validate the retrieved D_me using the D_me from the CoSSIR-Evans algorithm, the transformation from the R_e to the D_me is necessary due to the different definitions of R_e and D_me. Hence, we developed a conversion relationship between the R_e and D_me combined with different particle size distributions. Based on the integration of r and D over 14 408 PSDs, multiple groups of R_e and D_me are calculated according to Eqs. (5) and (6). Numerical fitting was used to build a relationship between the R_e (independent variable) and D_me (dependent variable) over 14 408 PSDs. The conversion formula between R_e and D_me is demonstrated as follows:

$$\begin{array}{}\text{(7)}& D_{\text{me}}=a+b{R}_{\mathrm{e}}+c{R}_{\mathrm{e}}^{\mathrm{2}},\end{array}$$

where a, b, and c are fitting coefficients obtained by numerical fitting and provided as input. Based on this relationship, we have unified all the R_e into D_me for comparability for the Voronoi, sphere, and column ICS models.

Figure 1 The overall flowchart of the retrieval of the IWP and D_me of ice clouds based on the Voronoi, column, and sphere ICS models.

3.1 Parameterization of ice cloud optical properties

To apply the Voronoi, sphere, and column ICS models to the RSTAR radiative transfer model for forwarding simulation, three parameterization schemes (referred to as the Voronoi, sphere, and column ICS schemes hereafter) for the scattering properties of ice clouds in the terahertz spectrum need to be constructed. The ice cloud optical properties depend on the single-scattering properties of the ICS and the PSDs. The effective ice cloud particle size measures the average ice cloud particle size for a given PSD. Different methods have been used to determine the effective ice cloud size, and according to Baum et al. (2005a, b), the effective particle diameter for a given PSD is determined by Eq. (8):

$$\begin{array}{}\text{(8)}& D_{\mathrm{e}}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\displaystyle \frac{{\int }_{L_{\text{min}}}^{L_{\text{max}}}V\left(L\right)n\left(L\right)\mathrm{d}L}{{\int }_{L_{\text{min}}}^{L_{\text{max}}}A\left(L\right)n\left(L\right)\mathrm{d}L}},\end{array}$$

where D_e is the effective particle diameter, and V and A are the volume and projected area of Voronoi and sphere models. Then, the spectral ice cloud optical properties (mass-averaged extinction coefficients, single-scattering albedo, asymmetry factor, and mass-averaged absorption coefficients) for the Voronoi and sphere ICS schemes are calculated for all PSDs given by Eqs. (9)–(12):

$$\begin{array}{}\text{(9)}& {\displaystyle }{K}_{\text{ext}}\left(\mathit{\lambda }\right)={\displaystyle \frac{{\int }_{L_{\text{min}}}^{L_{\text{max}}}{Q}_{\text{ext}}(\mathit{\lambda },L)A\left(L\right)n\left(L\right)\mathrm{d}L}{{\mathit{\rho }}_{\text{ice}}{\int }_{L_{\text{min}}}^{L_{\text{max}}}V\left(L\right)n\left(L\right)\mathrm{d}L}},\text{(10)}& {\displaystyle }\mathit{\varpi }\left(\mathit{\lambda }\right)={\displaystyle \frac{{\int }_{L_{\text{min}}}^{L_{\text{max}}}{Q}_{\text{sca}}(\mathit{\lambda },L)A\left(L\right)n\left(L\right)\mathrm{d}L}{{\int }_{L_{\text{min}}}^{L_{\text{max}}}{Q}_{\text{ext}}(\mathit{\lambda },L)A\left(L\right)n\left(L\right)\mathrm{d}L}},\text{(11)}& {\displaystyle }g\left(\mathit{\lambda }\right)={\displaystyle \frac{{\int }_{L_{\text{min}}}^{L_{\text{max}}}g(\mathit{\lambda },L){\mathit{\sigma }}_{\text{sca}}(\mathit{\lambda },L)n\left(L\right)\mathrm{d}L}{{\int }_{L_{\text{min}}}^{L_{\text{max}}}{\mathit{\sigma }}_{\text{sca}}(\mathit{\lambda },L)n\left(L\right)\mathrm{d}L}},\text{(12)}& {\displaystyle }{K}_{\text{abs}}\left(\mathit{\lambda }\right)={\displaystyle \frac{{\int }_{L_{\text{min}}}^{L_{\text{max}}}{Q}_{\text{abs}}(\mathit{\lambda },L)A\left(L\right)n\left(L\right)\mathrm{d}L}{{\mathit{\rho }}_{\text{ice}}{\int }_{L_{\text{min}}}^{L_{\text{max}}}V\left(L\right)n\left(L\right)\mathrm{d}L}},\end{array}$$

where K_ext(λ) represents spectral-mass-averaged extinction coefficients (m² g⁻¹), ϖ(λ) is spectral single-scattering albedo, g(λ) is a spectral asymmetry factor, and K_abs(λ) represents spectral-mass-averaged absorption coefficients (m² g⁻¹). Q_ext, g, Q_sca, and Q_abs are extinction efficiency, asymmetry factor, scattering efficiency, and absorption efficiency for Voronoi, sphere, and column models. Based on the parameterized scattering properties of ice clouds and the effective particle diameter of ice clouds, we developed a parameterization scheme for the scattering properties of ice clouds by establishing the spectral bulk-scattering properties of ice clouds as functions of the effective particle diameter of ice clouds. The least squares method is used to obtain first-order and third-order polynomial fitting equations according to Eqs. (13)–(16):

$$\begin{array}{}\text{(13)}& {\displaystyle }{K}_{\text{ext}}\left(\mathit{\lambda }\right)=a_{\mathrm{0}}+{\displaystyle \frac{a_{\mathrm{1}}}{D_{\mathrm{e}}}},\text{(14)}& {\displaystyle }\mathit{\varpi }\left(\mathit{\lambda }\right)=b_{\mathrm{0}}+b_{\mathrm{1}}{D}_{\mathrm{e}}+b_{\mathrm{2}}{D}_{\mathrm{e}}^{\mathrm{2}}+b_{\mathrm{3}}{D}_{\mathrm{e}}^{\mathrm{3}},\text{(15)}& {\displaystyle }g\left(\mathit{\lambda }\right)=c_{\mathrm{0}}+c_{\mathrm{1}}{D}_{\mathrm{e}}+c_{\mathrm{2}}{D}_{\mathrm{e}}^{\mathrm{2}}+c_{\mathrm{3}}{D}_{\mathrm{e}}^{\mathrm{3}},\text{(16)}& {\displaystyle }{K}_{\text{abs}}\left(\mathit{\lambda }\right)=d_{\mathrm{0}}+d_{\mathrm{1}}{D}_{\mathrm{e}}+d_{\mathrm{2}}{D}_{\mathrm{e}}^{\mathrm{2}}+d_{\mathrm{3}}{D}_{\mathrm{e}}^{\mathrm{3}},\end{array}$$

where a_i, b_j, c_j, and d_j ($i=\mathrm{0},\mathrm{1}$; $j=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}$) are fitting coefficients and are spectral functions of the terahertz wavelength. Values of the above coefficients for the Voronoi ICS scheme are listed in Appendix A (Tables A.1–A.4).

3.2 Look-up table

Before constructing the LUT, the sensitivity of different terahertz brightness temperatures to the IWP and D_me of ice clouds needs to be analysed to build a representative and efficient LUT. According to the sensitivity results (see Sect. 4.3), the BTDs between the cloudy and clear-sky conditions at 380, 640, and 874 GHz frequencies are simplified to BTD₁, BTD₂, and BTD₃, respectively, and the differences between the two BTDs are represented by a dash. Based on the Voronoi, sphere, and column ICS schemes, the RSTAR is used to construct three LUTs (Voronoi, sphere, and column LUTs) of BTD_2–3 and BTD_1–3 for cloudy- and clear-sky conditions. For the cloudy-sky LUT, the number of the IWP is set to 12, and the range of values is defined in log10 space from 0 to 3.36 in steps of 0.28. The number of R_e's is set to six, and the range of values is defined in log10 space from 1.6 to 3.1 in steps of 0.25. For the clear-sky LUT, the number of total ozone columns is set to seven, and the range of values is from 200 to 500 in steps of 50. To be consistent with aircraft observations of terahertz brightness temperatures, BTD_2–3 and BTD_1–3 are simulated at the mean altitude (20 km). The US standard atmospheric profile is used in the simulation, and the cloud top temperature is assumed to be the same as the atmospheric temperature at that level. The surface is assumed to be a black body (emissivity equals 1) with a temperature of 288.15 K.

3.3 Optimal estimation inversion method

Based on the terahertz BTD LUTs for the Voronoi, sphere, and column ICS schemes established in the previous section, the IWP and D_me are retrieved using an optimization method with Gaussian–Newtonian non-linear iterations. Based on the atmospheric radiative transfer transmission in the terahertz spectrum, if the background field under clear-sky conditions is subtracted from the cloudy-sky condition, the TOC BTDs are only functions of the cloud microphysical property parameters, which are given by Eqs. (17) and (18):

$$\begin{array}{}\text{(17)}& {\displaystyle }Y=F\left(X\right)+\mathit{ϵ},\text{(18)}& {\displaystyle }X=\left(\begin{array}{c}\text{IWP},\\ D_{\text{me}}\end{array}\right),Y=\left(\begin{array}{c}{\text{BTD}}_{\mathrm{2}-\mathrm{3}},\\ {\text{BTD}}_{\mathrm{1}-\mathrm{3}}\end{array}\right),\end{array}$$

where X is the vector matrix composed of the variables of the IWP and D_me to be solved. Y is the vector composed of the two BTDs and the uncertainty vector. The vector ϵ represents the uncertainties that are attached to the measurements (i.e. instrumental accuracy) and to the radiative transfer forward model (i.e.approximation errors in the radiative transfer model). Following Marks and Rodgers (1993), a good convergence can be obtained when the value of the cost function is lower than the size of the measurement vector. Since there is no robust a priori value for IWP and D_me, we selected an average value as the initial value for an a priori value of the IWP and D_me. Assuming that Y is the value of two BTDs measured by the aircraft, the inverse X is the minimum value for solving Eq. (19):

$$\begin{array}{}\text{(19)}& {\displaystyle }J\left(X\right)=min\mathit{\{}\Vert y-F(X)\Vert +\mathit{\gamma }\Vert X-X_{\mathrm{a}}\Vert \mathit{\}},\end{array}$$

where X_a is a vector consisting of the prior estimates of the IWP and R_e, γ is a Lagrangian operator, min denotes the minimal value of solving this function, and the value of X at the minimal value is considered to be the retrieved result. For the solution method of the non-linear problem in the inversion process, the Newton non-linear iterative method is usually used. It has a faster inversion speed and higher solution accuracy, and its iterative form follows Eq. (20):

$$\begin{array}{}\text{(20)}& {\displaystyle }{X}_{i+\mathrm{1}}=X_i-J^{\prime \prime }\left(X_i{)}^{-\mathrm{1}}{J}^{\prime }\right(X_i),\end{array}$$

where the i represents the number of iterations, and the superscripts denote the first-order and second-order derivations. Then the iterations start from the a priori initial estimate until the convergence criterion is satisfied or when the number of iterations satisfies the required number of iterations. Then the iteration is stopped, and the solution results are obtained. The final inversion results are validated using the results from the CoSSIR-MCBI algorithm.

4 Results

4.1 Single-scattering property database for the Voronoi, sphere, and column ICS models

Figure 2 compares the extinction efficiency, single-scattering albedo, and asymmetry factor, varying with the SZP for the Voronoi, column, and sphere ICS models at 325 and 874 GHz. For small ice particles with SZPs less than 0.1, the single-scattering properties are small and barely influenced by the shape of the ice particles. This is because the single-scattering properties of ice particles are close to Rayleigh scattering when the SZP is small. As particle size increases, particle scattering is predominantly Mie scattering, and the sensitivity of the single-scattering properties to the ice crystal habits becomes pronounced so that the ice crystal shape contributes to the large differences for large particle sizes. The single-scattering properties of the Voronoi and column models vary more smoothly than those of the sphere model. As shown in Fig. 2, the smooth curve of the Voronoi model reflects the notion that the FDTD and DDA methods provide good continuity for the transition between different size parameters of the Voronoi model. For large particles (SZP>3) at 874 GHz, the Voronoi model has a higher extinction efficiency than the other two models. At 325 GHz, the single-scattering albedo of the Voronoi model is close to the sphere model and is higher than the column model. At 874 GHz, the single-scattering albedo of the Voronoi model is close to the column model and is higher than the sphere model. The Voronoi model has a higher asymmetry factor than the other two models at 325 and 874 GHz. On the one hand, the higher extinction efficiency and single-scattering albedo of the Voronoi ICS model for large particles are possibly due to the multifaceted shapes of the Voronoi ICS model, which can result in significant side and backward scattering and increase the scattered energy. On the other hand, the higher asymmetry factor of the Voronoi ICS model for large particles is possible because the scattered energy is dominated by diffraction. The diffracted energy is concentrated in the forward direction, leading to a large asymmetry factor of the Voronoi ICS model.

Figure 2 The extinction efficiency, single-scattering albedo, and asymmetry factor as functions of the SZP for the Voronoi (solid blue line), sphere (dashed red line), and column (dashed green line) ICS models with a refractive index of 1.78+0.005i in the (a, c, e) 325 GHz and 1.78+0.015i in the (b, d, f) 874 GHz frequencies.

Figure 3 shows the contours of the single-scattering properties for the Voronoi, sphere, and column ICS models over 20 terahertz wavelengths (0.03–0.3 cm) and 31 particle sizes (1–270 µm). The sharp changes in the single-scattering albedo for the Voronoi, sphere, and column ICS models can be seen from small to large ice particles. For small particles in the low-frequency channels (long wavelengths), their scattering is mainly Rayleigh scattering with significant absorption effects. The absorption energy is proportional to the volume of ice particles and is barely affected by the ice particle shape. In the high-frequency channels (short wavelengths), the Mie scattering plays a leading role, and the scattering function plays a dominant role. Especially for large ice particles, the single-scattering albedo is close to one, and the influence of the ice particle shape becomes obvious. Figure 4 shows the differences in extinction efficiency, single-scattering albedo, and asymmetry factor between the Voronoi and the other two ICS models as functions of the effective radius and terahertz wavelength. Overall, the extinction efficiency of the Voronoi model is higher than the other two ICS models. The asymmetry factor of the Voronoi model is higher than the other two ICS models in the low-frequency channel for large ice particles.

Figure 3 The comparison of the extinction efficiency, single-scattering albedo, and asymmetry factor as functions of 31 particle sizes (1–270 µm) and 20 terahertz wavelengths (0.03–0.3 cm) for the (a, d, g) Voronoi, (b, e, h) sphere, and (c, f, i) column ICS models.

Figure 4 The (a, c, e) Voronoi minus sphere and the (b, d, f) Voronoi minus column ICS model differences in (a, b) extinction efficiency, (c, d) single-scattering albedo, and (e, f) asymmetry factor as functions of 31 particle sizes (1–270 µm) and 20 terahertz wavelengths (0.03–0.3 cm).

Figure 5 displays the scattering phase functions of the Voronoi, sphere, and column ICS models. The variation in the scattering phase function for the Voronoi model tends to be more dramatic compared to the sphere and column ICS models. The scattering phase function is axisymmetric about the 90^∘ scattering angle for small ice particles, which can be approximated as Rayleigh scattering. For small ice particles, the scattering phase function shows extreme values in the forward (0^∘) and backward (180^∘) directions and shows minimal values in both side directions (90 and 270^∘). As the ice particle size increases, the forward scattering is significantly larger than the side and backward scattering and gradually tends to be more Mie scattering. For the sphere model in the high-frequency channels, the forward scattering increases with the increase in particle size. In the low-frequency channels, the forward scattering remains almost constant with the increase in particle size. Compared to the Voronoi and sphere ICS models, the column ICS model shows smaller phase functions. For the Voronoi model, the low values of the scattering phase function move towards large scattering angles with the increase in particle size.

Figure 5 The scattering phase functions for ice particles with four sizes (R_e=30, 71, 107, and 153 µm) for the (a, d) Voronoi, (b, e) sphere, and (c, f) column ICS models with a refractive index of 1.78+0.005i in the 325 GHz and 1.78+0.015i in the 874 GHz, respectively.

4.2 Bulk-scattering properties of ice clouds in the terahertz region

Based on the single-scattering properties of the Voronoi, sphere, and column ICS models, the Voronoi, sphere, and column ICS schemes are developed in the terahertz region with the integration over both PSDs and terahertz frequencies. The calculated bulk-scattering properties of ice clouds include the mass extinction coefficients, single-scattering albedo, asymmetry factor, and mass absorption coefficients. Figure 6 compares the calculated bulk-scattering properties of the Voronoi, sphere, and column ICS schemes over 14 408 PSDs at four terahertz frequencies (325, 448, 664, and 874 GHz). The bulk-scattering properties of ice clouds depend on the effective diameter and the terahertz frequency. The single-scattering albedo increases linearly with the effective diameter for the sphere ICS scheme. For the Voronoi ICS scheme, the single-scattering albedo increases for the effective diameter smaller than 100 µm and approaches 0.95 with effective diameters exceeding 100 µm. The column ICS scheme has the lowest single-scattering albedo compared to the other two schemes for small effective diameters and gradually exceeds the other two schemes with the increase in the frequency and effective diameter. The mass extinction coefficients obtained from the three schemes show a uniformly positive correlation with the effective diameter and the terahertz frequency. At four terahertz frequencies, the mass extinction coefficients of the Voronoi ICS scheme are close to the column ICS scheme for all effective diameters. The sphere ICS scheme has the highest mass extinction coefficients compared to the other two schemes for all effective diameters and four terahertz frequencies. The Voronoi ICS scheme has a similar asymmetry factor to the sphere ICS scheme at low frequency and becomes the highest compared to the other two schemes at high frequency. The column ICS scheme has the lowest asymmetry factor compared to the other two schemes for all effective diameters. For all effective diameters, the mass absorption coefficients of the sphere ICS scheme have a stronger absorption effect than the other two schemes and decrease with increasing effective diameter. The Voronoi ICS scheme has the lowest mass absorption coefficients compared to the other two schemes for effective diameter smaller than 100 µm and becomes higher than the column ICS scheme with increasing effective diameters. Overall, the bulk-scattering properties for all schemes increase with increasing frequencies. This is mainly because the scattering properties of ice clouds are more significant as the wavelength decreases.

Figure 6 The comparison of the parameterized single-scattering albedo, mass extinction coefficients, asymmetry factor, and mass absorption coefficients as functions of effective diameters for 325, 448, 664, and 874 GHz for the sphere (blue line), Voronoi (orange line), and column (green line) ICS models.

4.3 Sensitivity results

We discuss the sensitivity of the TOC BTD_2–3 and BTD_1–3 to the IWP and D_me based on the Voronoi, sphere, and column ICS schemes, respectively, in the RSTAR model. Figure 7 shows that the BTD_2–3 and BTD_1–3 are positively correlated with the IWP and D_me. For the Voronoi and column ICS models, the BTD_2–3 and BTD_1–3 show monotonically increasing relationships with the increase in IWP. For ice clouds with the D_me larger than 200 µm, the 50 to 1000 µm particle sizes can lead to approximately 0–20 K BTD_2–3 and 0–30 K BTD_2–3. The BTD_2–3 and BTD_1–3 increase with the increase in the D_me and are close to a constant value for D_me larger than 600 µm. As shown in Fig. 8, the Voronoi ICS scheme has higher BTD_2–3 and BTD_1–3 compared to the sphere and column ICS scheme. In summary, the TOC BTD_2–3 and BTD_1–3 have a strong sensitivity to the IWP and D_me, especially for moderate to large ice particles and large IWP. The different ICS schemes vary considerably in their modelling of terahertz brightness temperatures.

Figure 7 The BTD_2–3 and BTD_1–3 for the (a, d) Voronoi, (b, e) sphere, and (c, f) column ICS models as functions of the IWP and D_me, respectively.

Figure 8 The BTD_2–3 and BTD_1–3 differences for the (a, c) Voronoi-minus-sphere and (b, d) Voronoi-minus-column ICS models as functions of the IWP and D_me, respectively.

Figure 9 exhibits the variation in the TOC BTDs at different IWPs and D_me for 380, 640, and 874 GHz frequencies based on the RSTAR model. The BTD_1–3 has a strong sensitivity to the IWP and increases with the increase in the IWP, and the slope shows that the increase is pronounced when the IWP is less than 600 g m⁻². The BTD_2–3 has a strong sensitivity to the D_me and increases with the increasing D_me. By comparing the LUTs of the Voronoi, sphere, and column ICS models, it is found that at the same BTD_2–3, the sphere LUT has the highest IWP compared to the column and Voronoi LUT, while the Voronoi LUT has the lowest. The D_me of the sphere LUT is smaller than that of the column and Voronoi LUT for the same BTD_1–3. This is mainly due to the higher single-scattering albedo and asymmetry factor of the Voronoi ICS model at low frequencies, resulting in stronger forward-scattering energy and a larger BTD than the sphere and column ICS model for the same IWP and D_me.

Figure 9 The LUT of BTD_2–3 and BTD_1–3 for the (a) Voronoi, (b) column, and (c) sphere ICS models varying with the logarithm of IWP and D_me. Grey dots in circles represent the randomly generated 2000 test data from the RSTAR model.

4.4 Inversion and validation results

To develop the inversion method from the CoSSIR measurements of the terahertz brightness temperature, we need to perform a quantitative simulation to test the accuracy of the algorithm. Figure 10 shows the validation results of the inversion algorithm using 2000 groups of BTD_2–3 and BTD_1–3 simulated by RSTAR using 2000 randomly generated IWP and D_me. The results show that most of the validation results of the red density lie on the 1:1 line. The correlation coefficients of the IWP and D_me are 0.94 and 0.99, and the mean absolute errors (MAEs) are 35.46 g m⁻² and 8.56 µm, respectively. The test results meet the formal scientific mission requirement of ice cloud remote sensing retrieval, according to Buehler et al. (2007). The high accuracy of the validation results proves the effectiveness and accuracy of the inversion algorithm.

Figure 10 The scatterplots of the randomly generated 2000 test data and the retrieved results. Panels (a) and (b) are for the IWP and D_me, respectively.

We adopt the CoSSIR-MCBI results as the benchmark for evaluating the retrieved IWP and D_me based on three schemes using the new inversion algorithm. To quantify the differences in the retrieval results between the three schemes and the CoSSIR-MCBI, we use the CoSSIR-MCBI results minus the retrieved results of the Voronoi, column, and sphere ICS schemes, respectively, to analyse their differences. Figure 11a and b show the CoSSIR-MCBI results on 19 July 2007 overlaid with the retrieval results of the IWP and D_me based on the Voronoi scheme. Red and blue dots, respectively, represent the Voronoi ICS scheme and CoSSIR-MCBI results. Figure 11a and b show that the overall performance shows good agreement for the Voronoi ICS scheme compared to the CoSSIR-MCBI results. The matching rates for the IWP and D_me are 80.3 % and 83.2 %, respectively. Figure 11c and d show the joint histogram of differences in the retrieved IWP and D_me between the Voronoi, sphere, and column schemes and the CoSSIR-MCBI results, respectively. Overall, the retrieved IWP and D_me parameters from the Voronoi ICS scheme agree well with the CoSSIR-MCBI results for all cases. Figure 12 shows the inversion results of IWP and D_me based on the optimal estimation algorithm from the CoSSIR aircraft-measured brightness temperature. For the Voronoi, sphere, and column ICS models, the correlation coefficients of the retrieved IWP are 0.87, 0.67, and 0.74, with MAEs of 22.38, 23.55, and 22.91 g m⁻² and RMSEs of 30.45, 35.22, and 32.64 g m⁻², respectively. The correlation coefficients of the D_me are 0.83, 0.64, and 0.76, with MAEs of 18.46, 21.19, and 19.91 µm and RMSEs of 24.57, 26.51, and 25.04 µm, respectively. Overall, the accuracy for the IWP and D_me based on the Voronoi ICS model is better than the other two ICS models compared to the CoSSIR-MCBI results. The sphere ICS model overestimates IWP and D_me compared with the validation data. According to the sensitivity results of Figs. 7 and 8, the Voronoi ICS scheme has higher BTD_2–3 and BTD_1–3 compared to the sphere and column ICS schemes, especially for large particles and IWP. This characteristic can be explicitly explained by the higher asymmetry factor of the Voronoi ICS model compared to the sphere and column ICS models. Thus, stronger forward-scattering energy can be detected for the Voronoi ICS model than for the other two models. The look-up table of the Voronoi ICS model can cover more IWP and D_me. The brightness temperature variations in the Voronoi-shaped ice clouds are more prominent and sensitive to the IWP and D_me. Therefore, the Voronoi ICS model results are better than the other two models. According to the look-up table as shown in Fig. 9, there are overlapping lines when D_me is small (D_me<40 µm) and large (D_me>140 µm). When BTD_2–3 and BTD_1–3 data fall under such overlapping lines of the look-up table, this overlapping region can lead to obtaining the same IWP and D_me when searching the look-up table.

Figure 11 The scattered red dots are the retrieved (a) IWP and (b) D_me (bottom row) validated by the results from the CoSSIR-MCBI algorithm (blue scattered dots). The joint histogram of differences in the retrieved (c) IWP and (d) D_me between the Voronoi, sphere, and column schemes and the CoSSIR-MCBI results, respectively.

Figure 12 The scatterplots of the retrieved IWP (a–c) and D_me (d–f) against the CoSSIR-MCBI results for the sphere (c, f), column (b, e), and Voronoi (a, d) ICS models.

5 Conclusions

In this study, we applied the irregularly shaped Voronoi ice crystal scattering (ICS) model to the ice cloud remote sensing retrieval of the ice water path (IWP) and particle size based on the Compact Scanning Submillimeter-wave Imaging Radiometer (CoSSIR) terahertz radiation measurements. The bulk-scattering property parameterization (Voronoi ICS scheme) in the terahertz region was developed based on the single-scattering properties of the Voronoi ICS database and 14 408 groups of particle size distributions from in situ observations. The Voronoi ICS scheme was applied to the atmospheric radiative transfer model RSTAR and compared with the sphere and column ICS schemes to carry out the sensitivity analysis. We analysed the sensitivity of brightness temperature differences between cloudy and clear skies (BTDs) at 380, 640, and 874 GHz to the IWP and particle size. Based on the sensitivity analysis results, we built three terahertz multi-channel BTD look-up tables (LUTs) for the Voronoi, sphere, and column ICS schemes using the RSTAR atmospheric radiative transfer model. Based on the three LUTs, we utilized the Gaussian–Newtonian non-linear optimization estimation method to retrieve the IWP and particle size from the CoSSIR terahertz radiation measurements. Finally, the retrieval results were evaluated by the IWP and median mass diameter (D_me) derived from the CoSSIR-MCBI algorithm. The main conclusions were obtained as follows.

The bulk-scattering properties of ice clouds in the terahertz region, including the mass extinction coefficients, single-scattering albedo, asymmetry factor, and mass absorption coefficients of the Voronoi ICS scheme, were applied to the RSTAR model and were compared with the sphere and column ICS schemes. The results showed that the Voronoi ICS scheme has a distinct feature of lower-absorption properties and higher asymmetry factors for larger particle sizes in the terahertz region. This feature could be related to the complex, multifaceted shape of the Voronoi ICS model and suggests that the Voronoi ICS scheme can produce relatively stronger forward scattering and fewer absorption effects compared with the sphere and column ICS schemes.

The sensitivity analysis showed that the BTDs between 640 and 874 GHz are sensitive to the IWP variation, and the BTDs of 380, 640, and 874 GHz are sensitive to the particle size. The atmospheric absorption peak near 380 GHz and the atmospheric window near 640 and 874 GHz can be effectively used for the IWP in the range of 50–200 g m⁻² and particle size of 50–300 µm.

A comparison of the results from the Voronoi, sphere, and column ICS models shows that the results from the Voronoi ICS model are better than the sphere and column ICS models. The IWP and D_me retrieved from the Voronoi ICS scheme showed higher consistency with CoSSIR-MCBI results than the other two schemes, with correlation coefficients of 0.87 and 0.83 for IWP and D_me, respectively. For the Voronoi ICS model, the mean absolute errors (MAEs) of the IWP and D_me are improved by 5.0 % and 12.8 %, and RMSE is improved by 13.5 % and 7.3 %, respectively. The sphere ICS scheme overestimates the IWP and D_me by up to the MAEs of 23.55 g m⁻² and 21.19 µm, respectively. This is mainly due to the differences in the absorption efficiency and asymmetry factor in the single-scattering properties of ice particles, which have a significant impact on the description of the scattering and radiative properties of ice clouds in the terahertz region.

In conclusion, the analysis of terahertz BTDs of 380, 640, and 874 GHz exhibits obvious sensitivity to the IWP and particle size of ice clouds, which could complement visible–infrared and microwave spectra. The present work provides the potential utility of inferring the IWP and particle size of ice clouds using BTDs of 380, 640, and 874 GHz. With the LUT for BTDs of 380, 640, and 874 GHz, the retrieval of terahertz ice cloud properties from airborne measurements based on the irregularly shaped Voronoi ICS models is newly developed. We find that the retrieval results based on the Voronoi ICS scheme present a better agreement with the CoSSIR-MCBI algorithm than the sphere and column ICS schemes. This study confirmed that the Voronoi ICS model has ice cloud inversion capabilities in the terahertz region, which may provide a reference for future use in aircraft-based and satellite-based terahertz ice cloud remote sensing applications.

Appendix A: Coefficients for Voronoi ICS scheme used in this study are tabulated in the Appendix

Table A1 Coefficients in the fitting of terahertz mass extinction coefficients (m² g⁻¹).

Table A2 Coefficients in the fitting of terahertz single-scattering albedo.

Table A3 Coefficients in the fitting of terahertz asymmetry factor.

Table A4 Coefficients in the fitting of terahertz mass-averaged absorption coefficients (m² g⁻¹).