Investigation of ice cloud modeling capabilities for the irregularly shaped Voronoi ice scattering models in climate simulations

Both weather–climate models and ice cloud remote sensing applications need to obtain effective ice crystal scattering (ICS) properties and the parameterization scheme. An irregularly shaped Voronoi ICS model has been suggested to be effective in remote sensing applications for several satellite programs, e.g., Himawari-8, GCOM-C (Global Change Observation Mission–Climate) and EarthCARE (Earth Cloud Aerosol and Radiation Explorer). As continuation work of Letu et al. (2016), an ice cloud optical property parameterization scheme (Voronoi scheme) of the Voronoi ICS model is employed in the Community Integrated Earth System Model (CIESM) to simulate the optical and radiative properties of ice clouds. We utilized the single-scattering properties (extinction efficiency, single-scattering albedo and asymmetry factor) of the Voronoi model from the ultraviolet to the infrared, combined with 14 408 particle size distributions obtained from aircraft measurements to complete the Voronoi scheme. The Voronoi scheme and existing schemes (Fu, Mitchell, Yi and Baum-yang05) are applied to the CIESM to simulate 10-year global cloud radiative effects during 2001–2010. Simulated globally averaged cloud radiative forcings at the top of the atmosphere (TOA) for Voronoi and the other four existing schemes are compared to the Clouds and the Earth's Radiant Energy System Energy Balanced and Filled (EBAF) product. The results show that the differences in shortwave and longwave globally averaged cloud radiative forcing at the TOA between the Voronoi scheme simulations and EBAF products are 1.1 % and 1.4 %, which are lower than those of the other four schemes. Particularly for regions (from 30^∘ S to 30^∘ N) where ice clouds occur frequently, the Voronoi scheme provides the closest match with EBAF products compared with the other four existing schemes. The results in this study fully demonstrated the effectiveness of the Voronoi ICS model in the simulation of the radiative properties of ice clouds in the climate model.

1 Introduction

Ice clouds cover about 20 %–30 % of the global area (Rossow and Schiffer, 1991; Wang et al., 1996; Stubenrauch et al., 2013), and they strongly affect the earth's energy budget and climate system mainly due to their optical and radiative properties (Liou, 1986, 1992; Baran, 2012; Ramaswamy and Ramanathan, 1989). The radiative properties of ice clouds mainly depend on their optical properties (e.g., scattering albedo and optical thickness), which are significantly influenced by the microphysical properties (e.g., ice particle sizes and habits) of ice clouds (Baran, 2009; Yang et al., 2015, 2018). Based on the accurate knowledge of ice particle sizes and habits, the single-scattering properties (e.g., extinction efficiency, single-scattering albedo and asymmetry factor) of ice particles can be calculated by using the light scattering computational method to develop the ice crystal scattering (ICS) model/database. Combined with the single-scattering properties of the ICS model/database and size distributions, the optical properties of ice clouds can be simulated through the parameterization scheme, and the radiative properties of ice clouds can be further simulated based on the radiative transfer theory.

At present, satellite remote sensing and weather–climate models are two effective ways to understand the ice cloud optical and radiative properties through the ice cloud optical property parameterization scheme in radiative transfer models. However, numerous field observations (Rossow and Schiffer, 1999; Lawson et al., 2006; Heymsfield et al., 2017; Lawson et al., 2019), e.g., the First International Satellite Cloud Climatology Project (ISCCP) Regional Experiment (FIRE-I) in 1986 and 1991 (Rossow and Schiffer, 1999) and the European cirrus experiment in 1989 (Liou, 1992), have shown that ice clouds contain a large variety of non-spherical ice particle sizes and habits, which can lead to inaccurate simulations of the optical and radiative properties of ice clouds in nature. Our understanding of how ice particle habits affect the optical and radiative properties of ice clouds is still limited (Heymsfield and Miloshevich, 2003; van Diedenhoven et al., 2014; van Diedenhoven, 2018). The current insufficient knowledge of the ice cloud microphysical properties and the ice particle single-scattering properties contributes to inadequate representation of the optical properties in the parameterization scheme, which can directly lead to uncertainties in the simulated radiative properties of ice clouds (Zhang et al., 2015; Yang et al., 2015, 2018; Yi et al., 2017; van Diedenhoven and Cairns, 2020). Thus, the accurate representation of the microphysical properties of ice clouds and ice particle single-scattering properties is essential for the parameterization of ice cloud optical properties and the study of the radiative properties of ice clouds in both satellite remote sensing and weather–climate models. In terms of current light scattering computational methods, it is still difficult for one specific method to accurately calculate the single-scattering properties for non-spherical particles with a different size parameter (SZP), which is defined as the ratio of the equivalent-volume sphere's circumference dimension (or π times particle maximum diameters) to the incident wavelength (Nakajima et al., 2009; Baran, 2012; Yang et al., 2015). The existing light scattering computational methods for non-spherical particles can be generally divided into the approximation method (AM) based on the ray-tracing techniques (Wendling et al., 1979), as well as the numerical simulation (NM) method based on the approximate solutions of Maxwell equations. The AM method is suitable for non-spherical particles with very large SZPs. The geometrical optics approximation (GOA) method is a typical AM method. This method can capture the halo phenomenon of large hexagonal ice particles in the visible wavelength. However, the AM method has difficulties in accurately simulating the single-scattering properties for particles with small and moderate SZPs. The NM method is suitable for particles with small SZPs and can be divided into the volume- and surface-based methods depending on how Maxwell equations are solved. The volume-based NM method includes the finite-difference time domain (FDTD) (Yee, 1966; Yang and Liou, 1996b) and discrete dipole approximation (DDA) methods (Draine and Flatau, 1994; Yurkin and Hoekstra, 2007). A typical method of the surface-based method is the T-matrix method (Havemann and Baran, 2001; Mishchenko and Travis, 1998). However, the NM method requires discretization of the whole volume/surface of the scatterer and has rather high computational demands (Nakajima et al., 2009), so it is difficult to efficiently calculate the single-scattering properties for particles with moderate and large SZPs. Later, the invariant imbedding T-matrix (II-TM) (Bi et al., 2013a; Bi and Yang, 2014) and physical-geometric optics hybrid (PGOH) method (Bi et al., 2011) were developed for particles with small to moderate SZPs. Combined with the advantages of the AM and NM methods, several improved GOA methods including the geometric optics integral equation (GOIE) (Yang and Liou, 1996a; Ishimoto et al., 2012) and improved geometric-optics method (IGOM) (Yang and Liou, 1995, 1996a; Bi et al., 2010) have been developed. The GOIE and IGOM methods are useful for particles with moderate SZPs, and therefore they can bridge the gap between the AM and NM methods. Thus, the light scattering computation of particles with different SZPs can be completed by a combination of the AM and NM methods.

With the development of the ICS model/database, numerous parameterization schemes of the ice cloud optical properties have been developed for use in ice cloud remote sensing and weather–climate model applications (Yang et al., 2015, 2018). In terms of weather–climate model applications, Fu (1996) developed a parameterization scheme (referred to as the Fu scheme hereafter) using the GOA-based ICS database for the randomly oriented hexagonal particle (Takano and Liou, 1989). The Fu scheme was subsequently applied to the Fu–Liou radiative transfer model for use in the climate models (Fu, 1996, 2007). Mitchell et al. (1996b, 2006) used the modified anomalous diffraction approximation (MADA) method (Mitchell and Arnott, 1994) to generate an ICS database for a habit mixture and completed a parameterization scheme (referred to as the Mitchell scheme hereafter) combined with the bimodal size distributions (Mitchell et al., 1996a). The Mitchell scheme was then employed in the National Center for Atmospheric Research Community Atmosphere Model (CAM). Yang et al. (2000) used the IGOM and FDTD methods to develop an ICS database for six ice particle habits. However, this database contains several inconsistencies in the spectral regions caused by differences in particle habits and computational methods. Later, Yang et al. (2013) utilized the Amsterdam DDA (Yurkin et al., 2007; Yurkin and Hoekstra, 2011), T-matrix (Mishchenko et al., 1996) and improved IGOM (Bi et al., 2009) methods to generate a spectrally consistent ICS database for 11 ice particle habits. Yi et al. (2013) employed the ICS database of Yang et al. (2013) and developed a parameterization scheme (referred to as the Yi scheme hereafter) for use in the CAM. For ice cloud remote sensing applications, Baum et al. (2005a, b) used the ICS database of Yang et al. (2000) to develop a parameterization scheme (referred to as the Baum-yang05 scheme hereafter) for the MODIS Collection 5 ice cloud product. C.-Labonnote et al. (2000, 2001) and Doutriaux-Boucher et al. (2000) developed an ICS database for the inhomogeneous hexagonal monocrystal (IHM) model containing embedded inclusions (air bubbles and aerosols) and developed a parameterization scheme for use in the ice cloud retrievals from the French satellite Polarization and Directionality of the Earth's Reflectance (POLDER) measurements (Deschamps et al., 1994). Ishimoto et al. (2012) and Letu et al. (2016) developed an ICS database by using a combination method of the FDTD, GOIE and geometric-optics method (GOM) for an irregularly shaped Voronoi model based on in situ microphysical measurements. Letu et al. (2016, 2020) demonstrated that the Voronoi model can effectively retrieve the ice cloud microphysical properties from satellite measurements. Furthermore, the Voronoi 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 programs (Illingworth et al., 2015), which will be launched in 2023. These studies demonstrated the superiority of the Voronoi model in the ice cloud remote sensing applications. However, the performance of the Voronoi model in the climate model simulations has not been investigated quantitatively.

Motivated by the abovementioned situations, this study aims to quantify the effects of the Voronoi model on the optical and radiative properties of ice clouds in a climate model in comparison with the other ice cloud optical property schemes (Fu, Mitchell, Yi and Baum-yang05). To achieve this goal, we develop an ice cloud optical property parameterization scheme (referred to as the Voronoi scheme hereafter) for the Voronoi model. The Voronoi scheme and the other schemes are employed in the Community Integrated Earth System Model (CIESM) (Lin et al., 2020) to simulate shortwave and longwave fluxes at the top of the atmosphere (TOA). The CIESM-based simulations from the five schemes are compared with the Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) products. This study addresses the following questions through the comparison. How different is the Voronoi scheme from the other schemes? How and to what extent does the Voronoi scheme affect the radiative effects of ice clouds? What are the possible reasons for the impacts of the Voronoi scheme?

This paper is organized as follows. Sections 2 and 3 introduce the data and methodology used in this study, and Sect. 4 demonstrates the influence of the Voronoi model on the cloud radiative properties through the radiative transfer model and climate model. Section 5 presents the summary and conclusions of this study.

Figure 1 Single-scattering properties (extinction efficiency, single-scattering albedo and asymmetry factor) of the Voronoi model from the composite method based on the FDTD, GOIE and GOM methods at the wavelengths of (a) 0.64 µm and (b) 2.21 µm.

2 Data

2.1 Single-scattering property database for the Voronoi model

In this study, the single-scattering property database of the Voronoi model developed by Ishimoto et al. (2012) and Letu et al. (2016) is used in the parameterization process. The single-scattering properties including the extinction efficiency, single-scattering albedo and asymmetry factor of the Voronoi model from the ultraviolet to the infrared are utilized to calculate the shortwave and longwave optical properties of ice clouds. The single-scattering properties of the Voronoi model are calculated by a combination of methods of FDTD, GOIE and GOM. This combination method implements a treatment of particle edge effects well (Ishimoto et al., 2012). With this treatment, the gaps between the results calculated by the FDTD and those by the GOIE are relatively small, which can lead to consistency in the single-scattering properties of ice particles. Figure 1a and b show the extinction efficiency, single-scattering albedo and asymmetry factor for Voronoi models that vary with SZPs at wavelengths of 0.64 and 2.21 µm, respectively. Note that the FDTD and GOIE methods are used for small (SZP < 40) and moderate (SZP < 300) ice particles, respectively, and the GOM method is used for large (SZP > 300) particles. The extinction efficiency at wavelengths of 0.64 and 2.21 µm has a peak value when the SZP approximately equals 10 and decreases to be a constant value of 2 with increasing SZPs larger than 100. The single-scattering albedo at both wavelengths is close to 1, which is related to the high values of the real part in the refractive index. The asymmetry factor decreases with the increasing SZPs at the wavelength of 0.64 µm. At the wavelength of 2.21 µm, the asymmetry factor increases for the SZPs smaller than 10 and larger than 300.

2.2 Aircraft-based measurements of the particle size distributions

To generate the parameterization scheme of ice cloud optical properties for application in the climate model simulations, particle size distributions (PSDs) of the Voronoi model need to be assumed in ice clouds. In this study, we utilized 14 408 PSDs derived from in situ aircraft-based measurements obtained in 11 field campaigns (available http://stc-se.com/data/bbaum/Ice_Models/microphysical_data.html, last access: 1 February 2021) (Heymsfield et al., 2013). These data confirm that the particle phase is unambiguously ice after filtering by cloud temperature ($T\le -\mathrm{40}$ ^∘C). For the fitting of PSDs for the Voronoi model, we adopt the gamma distribution form as follows:

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

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

Figure 2 Variations in ice particle size distributions for different temperatures.

As shown in Fig. 2, the particle concentration decreases with increasing L for all ranges. When the temperature is between −60 and −55^∘C, the particle concentration is the largest, and it decreases most sharply with increasing L. When the temperature is between −45 and −40^∘C, as well as between −65 and −70^∘C, the particle concentration is the smallest, and it decreases slowly with increasing L.

2.3 Satellite data used in the validation

To evaluate the cloud radiative effects for different schemes of the ICS model, we adopted the CERES EBAF Ed4.1 products (available https://ceres.larc.nasa.gov/data/, last access: 6 May 2021) (Draine and Flatau, 1994; Doelling et al., 2016) from 2001 to 2010 as validation data. The “toa_cre_sw_mon” and “toa_cre_lw_mon” EBAF products are used for comparisons with the simulated shortwave and longwave cloud radiative effects from all five schemes. The “toa_cre_sw_mon” and “toa_cre_lw_mon” products are monthly mean shortwave and longwave cloud radiative effects at the TOA, and they are calculated as all-sky fluxes minus total region clear-sky fluxes for shortwave and longwave spectra. The “toa_sw_all_mon” and “toa_lw_all_mon” EBAF products are monthly mean all-sky outgoing shortwave and longwave fluxes at the TOA, and they are used for comparisons with simulated upwelling shortwave flux at the TOA (FSUTOA) and upwelling longwave flux at the TOA (FLUTOA) from all schemes. The “sfc_sw_down_all_mon” and “sfc_lw_down_all_mon” are monthly mean all-sky downwelling shortwave and longwave fluxes at the surface, and they are used for comparisons with simulated downwelling shortwave flux at the surface (FSDS) and downwelling longwave flux at the surface (FLDS) from all schemes. The spatial and temporal resolution of EBAF data is ${\mathrm{1}}^{\circ }\times {\mathrm{1}}^{\circ }$ latitude by longitude and monthly means.

Figure 3 Flowchart of the investigation of ice cloud modeling capabilities for the irregularly shaped Voronoi models in climate simulations.

3 Methodology

The main flowchart of this study is described in Fig. 3. Firstly, we developed the parameterization scheme of the Voronoi ICS model by using the aforementioned single-scattering properties of the Voronoi ICS database and large amounts of PSDs. Then, the Voronoi scheme and the other four existing schemes (Fu, Mitchell, Yi and Baum-yang05) were evaluated through simulations of shortwave upward and downward flux profiles in the rapid radiative transfer model for general circulation models (RRTMG). The RRTMG was utilized to understand how the different optical properties of the five schemes influence the upward and downward fluxes under several idealized conditions. Furthermore, all schemes were applied to the CIESM to simulate global shortwave and longwave cloud radiative forcing at the TOA from 2000 to 2010, among which the first year was removed for reaching the equilibrium state and the last 10 years were evaluated by EBAF products. The CIESM was employed to evaluate the effectiveness of the Voronoi ICS model in the simulations of ice cloud radiative properties compared with the other four schemes in the climate system.

3.1 Parameterization of ice cloud optical properties

To better understand the ice cloud modeling capabilities of the Voronoi model in the climate model and explain how ice clouds play a role in the climate system, it is necessary to introduce the main scattering parameters to evaluate the Voronoi ICS model through the parameterization scheme. The main radiative transfer processes can be simply attributed to extinction, scattering and absorption coefficients (Liou, 1986, 1992), which are calculated by Eq. (2).

$$\begin{array}{}\text{(2)}& {\mathit{\beta }}_{\text{e,s,a}}=\underset{L_{\text{min}}}{\overset{L_{\text{max}}}{\int }}{\mathit{\sigma }}_{\text{e,s,a}}n\left(L\right)\mathrm{d}L,\end{array}$$

where β_e,s,a is the extinction, scattering and absorption coefficients and σ is the cross section (See Table A1 for a list of acronyms). The single-scattering albedo and co-albedo can be defined as the ratio of the scattering and absorption coefficients to the extinction coefficient in the form of Eq. (3).

$$\begin{array}{}\text{(3)}& \mathit{\varpi }={\displaystyle \frac{{\mathit{\beta }}_{\text{s}}}{{\mathit{\beta }}_{\text{e}}}},\text{or}\mathrm{1}-\mathit{\varpi }={\displaystyle \frac{{\mathit{\beta }}_{\text{a}}}{{\mathit{\beta }}_{\text{e}}}},\end{array}$$

where ϖ and 1−ϖ are single-scattering albedo and co-albedo, respectively. Based on the extinction coefficient, the optical depth can be defined by Eq. (4).

$$\begin{array}{}\text{(4)}& \mathit{\tau }=\underset{z}{\overset{\mathrm{\infty }}{\int }}{\mathit{\beta }}_{\text{e}}\mathrm{d}z,\end{array}$$

where τ is the optical depth, and z is the outer boundary of the atmosphere. In the assumption of plane-parallel atmospheres, changes in the diffuse intensity penetrating from below the layer considering multiple scattering processes can be given by Eq. (5).

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& J\left(\mathit{\tau };\mathit{\mu };\mathit{\varphi }\right)={\displaystyle \frac{\mathit{\varpi }}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}I\left(\mathit{\tau };{\mathit{\mu }}^{\prime };{\mathit{\varphi }}^{\prime }\right)P\left(\mathit{\mu },\mathit{\varphi };{\mathit{\mu }}^{\prime },{\mathit{\varphi }}^{\prime }\right)\mathrm{d}{\mathit{\mu }}^{\prime }\mathrm{d}{\mathit{\varphi }}^{\prime }\\ & +{\displaystyle \frac{\mathit{\varpi }}{\mathrm{4}\mathit{\pi }}}{F}_{\mathrm{\Theta }}P\left(\mathit{\mu },\mathit{\varphi };-{\mathit{\mu }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}}\right)e^{-\mathit{\tau }/{\mathit{\mu }}_{\mathrm{0}}}+\left(\mathrm{1}-\mathit{\varpi }\right)B\left[T\left(\mathit{\tau }\right)\right]\end{array},\end{array}$$

where P is the phase function corresponding to a volume of ice particles. $P(\mathit{\mu }\mathit{\varphi };{\mathit{\mu }}^{\prime }{\mathit{\varphi }}^{\prime })$ denotes the redirection of the incoming intensity defined by (${\mathit{\mu }}^{\prime }{\mathit{\varphi }}^{\prime })$ to the outgoing intensity defined by (μϕ). I indicates the total (direct plus diffuse) radiance, B indicates Planck's function associated with thermal emissions, and Θ is the scattering angle. Therefore, the extinction coefficients, single-scattering albedo and phase function are fundamental driving parameters within the transfer of diffuse intensity.

Based on these principles, in this study, we completed the Voronoi scheme by using the single-scattering properties of the Voronoi model and 14 408 groups of PSDs data. The parameterization of the ice cloud optical properties for the Voronoi scheme are developed following Eqs. (6)–(15). Firstly, the spectral ice cloud optical properties (mass-averaged extinction coefficients, single-scattering albedo and asymmetry factor) for the Voronoi scheme are calculated for all PSDs given by Eqs. (6)–(9).

$$\begin{array}{}\text{(6)}& {\displaystyle }{D}_{\text{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{max}}}^{L_{\text{max}}}A\left(L\right)n\left(L\right)\mathrm{d}L}},\text{(7)}& {\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{(8)}& {\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{(9)}& {\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}}\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 volume and projected area of the Voronoi models. K_ext(λ) is spectral mass-averaged extinction coefficients (m² g⁻¹), ϖ(λ) is the spectral single-scattering albedo, and g(λ) is the spectral asymmetry factor. Q_ext, g and Q_sca are extinction efficiency, asymmetry factor and scattering efficiency for Voronoi models.

Then, based on the spectral bulk optical properties including K_ext(λ), ϖ(λ) and g(λ) of ice clouds, the band-averaged optical properties are calculated to apply the parameterization scheme in RRTMG and CIESM following Eqs. (10)–(12).

$$\begin{array}{}\text{(10)}& {\displaystyle }{\stackrel{\mathrm{̃}}{K}}_{\text{ext}}={\displaystyle \frac{{\int }_{{\mathit{\lambda }}_{\text{min}}}^{{\mathit{\lambda }}_{\text{max}}}{\mathit{\beta }}_{\text{ext}}\left(\mathit{\lambda }\right)E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{{\mathit{\lambda }}_{\text{min}}}^{{\mathit{\lambda }}_{\text{max}}}E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\text{(11)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\varpi }}={\displaystyle \frac{{\int }_{{\mathit{\lambda }}_{\text{min}}}^{{\mathit{\lambda }}_{\text{max}}}\mathit{\varpi }\left(\mathit{\lambda }\right)E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{{\mathit{\lambda }}_{\text{min}}}^{{\mathit{\lambda }}_{\text{max}}}E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\text{(12)}& {\displaystyle }\stackrel{\mathrm{̃}}{g}={\displaystyle \frac{{\int }_{{\mathit{\lambda }}_{\text{min}}}^{{\mathit{\lambda }}_{\text{max}}}g\left(\mathit{\lambda }\right)E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{{\int }_{{\mathit{\lambda }}_{\text{min}}}^{{\mathit{\lambda }}_{\text{max}}}E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\end{array}$$

where ${\stackrel{\mathrm{̃}}{K}}_{\text{ext}}$, $\stackrel{\mathrm{̃}}{\mathit{\varpi }}$ and $\stackrel{\mathrm{̃}}{g}$ are band-averaged and mass-averaged extinction coefficients, single-scattering albedo, and asymmetry factor for the Voronoi scheme, respectively. E is assigned by the solar constant provided by Chance and Kurucz (2010) for the shortwave spectrum and is replaced with the Planck function B(T) for the longwave spectrum, and T is an assumed cloud temperature of 233 K according to Liou (1992). The coefficients of the polynomial expressions of the ice cloud band-averaged optical properties as functions of D_e are determined in each band interval to develop the Voronoi scheme for shortwave and longwave spectra as shown in Eqs. (13)–(15).

$$\begin{array}{}\text{(13)}& {\displaystyle }{\stackrel{\mathrm{̃}}{K}}_{\text{ext}}=a_{\mathrm{0}}+a_{\mathrm{1}}/D_{\text{e}}+a_{\mathrm{2}}/D_{\text{e}}^{\mathrm{2}},\text{(14)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\varpi }}=b_{\mathrm{0}}+b_{\mathrm{1}}{D}_{\text{e}}+b_{\mathrm{2}}{D}_{\text{e}}^{\mathrm{2}}+b_{\mathrm{3}}{D}_{\text{e}}^{\mathrm{3}},\text{(15)}& {\displaystyle }\stackrel{\mathrm{̃}}{g}=c_{\mathrm{0}}+c_{\mathrm{1}}{D}_{\text{e}}+c_{\mathrm{2}}{D}_{\text{e}}^{\mathrm{2}}+c_{\mathrm{3}}{D}_{\text{e}}^{\mathrm{3}},\end{array}$$

where a, b and c are coefficients as functions of band intervals.

In terms of the other four existing schemes, the band-averaged optical properties of the Mitchell, Yi and Baum-yang05 schemes are developed as functions of D_e following Eqs. (13–15). Coefficients of the Mitchell scheme can be obtained from the CIESM. Values of coefficients for the Yi and Baum-yang05 schemes are listed in Appendix A (Tables A1, A2, A3, and A4) provided from Zhao et al. (2018). Coefficients of the Fu scheme (default scheme in RRTMG) are obtained from the existing ice cloud band-averaged optical properties from RRTMG. Formulation of the Fu scheme is similar to the Mitchell and Baum-yang05 schemes except for using the generalized effective diameter (Fu, 1996). The generalized effective diameter of the Fu scheme is unified into D_e for comparability.

Figure 4 The comparison of (top row) mass extinction coefficients, (middle row) single-scattering albedo and (bottom row) asymmetry factor as functions of effective diameter and wavelength for the (a) Voronoi, (b) Mitchell, (c) Fu, (d) Baum-yang05 and (e) Yi schemes.

Table 1 Shortwave and longwave bands set in the RRTMG.

3.2 RRTMG and CIESM simulation experiments

The version of the RRTMG used in this study is the current version of the radiative transfer code applied in the CIESM (Mlawer et al., 1997; Iacono et al., 2008; Clough et al., 2005; (available from http://rtweb.aer.com/, last access: 5 January 2021). RRTMG utilizes the correlated-k approach to calculate shortwave fluxes and heating rates efficiently and accurately for application in climate models. The version of RRTMG utilizes a two-stream method for radiative transfer calculation. RRTMG has 14 bands for the shortwave spectrum and 16 longwave bands (see Table 1). Since the wavelength range is from 0.2 to 15 µm for the single-scattering property database of the Voronoi ICS model, ice cloud bulk optical properties of the default scheme (the Mitchell scheme) remain unchanged when bands are larger than 15 µm. To quantify the radiative flux differences caused by the five schemes under the same conditions, we design an assumed ice cloud case in a standard tropical atmospheric profile in the RRTMG. The RRTMG settings are as follows: ice effective radius is set to 45 µm, the ice water path is set to 60 g m⁻², the ice cloud top pressure/height is between 125.1 and 245.5 hPa, the cloud fraction is 50 %, and the solar zenith angles is set to 60^∘. The vertical resolution is 60 levels for the standard tropics. The RRTMG is run by the five ice cloud schemes under the same conditions; thus relative difference in fluxes can be explained by difference among the five schemes. The five schemes are implemented in the CIESM to calculate the FSDS, FLDS, FSUTOA, FLUTOA, shortwave cloud forcing (SWCF) and longwave cloud forcing (LWCF). SWCF is defined as Eq. (16), LWCF is defined the same as SWCF but for the longwave spectrum.

$$\begin{array}{}\text{(16)}& \text{SWCF}=F_{\text{cloudy}}-F_{\text{clear}},\end{array}$$

where F_cloudy and F_clear are the difference between downward and upward fluxes for cloudy and clear conditions, respectively. The CIESM is operated in two ways: (1) the CIESM is run with the default Mitchell scheme for ice clouds and the default water cloud scheme to obtain SWCF and LWCF for the Mitchell scheme, and (2) the CIESM is run by using the other four schemes (the Voronoi, Yi, Baum-yang05 and Fu schemes) in place of the Mitchell scheme, along with the default liquid water cloud scheme. Liquid water clouds adopt a spherical particle model, whose single-scattering properties are derived from the Lorenz–Mie theory (van de Hulst, 1957). Because the CIESM is unable to separate ice clouds from liquid clouds efficiently, the SWCF and LWCF for the five schemes are under the same liquid water cloud parameterization, and their differences are attributed to different ice habits and their scattering and absorption properties within the five schemes. The CIESM run is integrated for 11 years in 1-month increments; the initial first year is used for state initialization and stabilization, and the runs from the last 10 years were utilized for analysis. Horizontal and vertical resolution of CIESM run experiment is set to ${\mathrm{1.9}}^{\circ }\times {\mathrm{2.5}}^{\circ }$ and 31 levels, respectively. The run is driven by prescribed climatological sea surface temperature and sea ice fraction with an annual cycle in the year 2000.

Figure 5 The (a) Mitchell, (b) Fu, (c) Baum-yang05 and (d) Yi schemes minus the Voronoi scheme differences (%) in (top row) mass extinction coefficients, (middle row) single-scattering albedo and (bottom row) asymmetry factor as functions of effective diameter and wavelength.

4 Results and discussions

4.1 Band-averaged optical properties of the ice cloud

Based on the integration over both PSDs and band intervals, band-averaged bulk optical properties of the Voronoi scheme are compared with the Mitchell, Fu, Baum-yang05 and Yi schemes in Fig. 4. The differences in ice cloud optical properties between the Voronoi and the other four existing schemes are shown in Fig. 5. Since the Fu scheme uses generalized effective diameter, the remaining four schemes use D_e. Horizontal axes in both Figs. 4 and 5 are unified into D_e for comparability. Parameterized mass extinction coefficients, single-scattering albedo and asymmetry factor are plotted as functions of the D_e from 10 to 150 µm and 14 selected bands.

Figure 6 Shortwave (a) upward fluxes, (b) diffuse downward fluxes, (c) downward total fluxes and (d) net fluxes for the Voronoi, Mitchell, Fu, Baum-yang05 and Yi schemes and clear conditions (blue line) for standard tropical atmospheric profile. Graphics in dotted black box are magnified and displayed in (e)–(h).

In Fig. 4, mass extinction coefficients obtained from the five schemes show a uniformly negative correlation with D_e. Mass extinction coefficients exceed up to 0.2 m² g⁻¹ for D_e smaller than 20 µm and approach 0 for D_e larger than 100 µm. This could be partly related to the majority of small particles in PSDs for the five schemes. Note that there is a minimum of mass extinction coefficients between 3.08 and 3.85 µm. It could be because the real part of the refractive index reaches the minimum at 3 µm (Warren and Brandt, 2008; Yang et al., 2013). This could weaken the scattering and extinction efficiency of ice particles. The single-scattering albedos obtained from the five schemes approach 1.0 in visible bands (0.2–0.78 µm). In near-infrared bands (0.78–3.85 µm), the single-scattering albedos are inversely proportional to D_e. This result is related to the large real part in visible bands and small imaginary part of the complex refractive index of ice particles. The asymmetry factor obtained from the five schemes increases with increasing wavelength for all D_e. From the visible to near-infrared band, the asymmetry factor increases with the increasing D_e.

Figure 7 The 10-year globally averaged CIESM-based TOA SWCF simulations for the Voronoi, Mitchell, Fu, Baum-yang05 and Yi schemes and corresponding CERES EBAF products during 2001–2010.

Figure 8 Same as the Fig. 7 but for TOA LWCF.

In Fig. 5, differences in mass extinction coefficients between the Voronoi scheme and the other four schemes show that the Voronoi scheme has slightly larger values than the Fu, Yi and Baum-yang05 schemes but is smaller than the Mitchell scheme. For the single-scattering albedo, the Voronoi scheme has a slightly larger single-scattering albedo than the Fu and Mitchell schemes and lower single-scattering albedo than the Baum-yang05 and Yi schemes in infrared bands. This result may be because large ice particles are closer to geometric optics and have a greater proportion of absorption than small ice particles. The low asymmetry factor of the Voronoi scheme is because the multifaceted shapes of the Voronoi ice model can result in significant side and backward scattering and reduced forward scattered energy. Since the impacts of different size distribution assumptions on the bulk optical properties of ice cloud parameterization are negligible (Heymsfield et al., 2013, 2017), differences in band-averaged bulk optical properties between the five schemes are originally rooted in different habits of ice particles and their single-scattering properties.

Figure 9 Comparison of 10-year zonally averaged annual mean (a) SWCF and (b) LWCF between the Voronoi, Mitchell, Fu, Baum-yang05 and Yi schemes and CERES EBAF. EBAF minus the five schemes' differences in (c) SWCF and (d) LWCF.

4.2 RRTMG simulation results

After the parameterization, band-averaged optical properties of ice cloud from the five schemes (Fu, Mitchell, Yi, Baum-yang05 and Voronoi) are subsequently parameterized as functions of D_e and 14 bands as shown in Fig. 4. To illustrate and quantify the influence of optical properties of ice cloud on its radiative effects, an ideal experiment is designed to test the response of radiative flux to the five ice cloud schemes under the same idealized conditions. Band-averaged optical properties for the five schemes are subsequently implemented in RRTMG to simulate radiative fluxes under prescribed ice clouds in standard tropical profiles (Anderson et al., 1986) which have a high proportion of ice cloud coverage (Massie et al., 2002; Stubenrauch et al., 2013). According to observations of Hong and Liu (2015), top and bottom pressures of ice cloud layer are set to 125.1 and 245.5 hPa, respectively, the D_e is set to 45 µm, and ice water paths equal 60 g m⁻².

Figure 10 Box analysis of zonal distributions of 10-year annual mean SWCF (a) and LWCF (b) difference between the Mitchell, Voronoi, Fu, Baum-yang05 and Yi schemes and EBAF products, respectively.

Shortwave radiative flux profiles of cloudy sky for the five schemes and clear-sky conditions are shown in Fig. 6. Specific comparison of the five schemes inside the dotted black region are enlarged and shown in the bottom row. In Fig. 6, upward fluxes of the five schemes gradually increase from cloud bottom to cloud top, reaching the maximum at the cloud top. The Voronoi and Mitchell schemes have higher upward fluxes and lower downward diffuse fluxes than the other three schemes. Figure 6 shows 6–30 W m⁻² differences in TOA upward fluxes, 10–40 W m⁻² differences in surface downward diffuse fluxes, 10–30 W m⁻² differences in surface net fluxes and 8–42 W m⁻² differences in TOA net fluxes owing to the five different ice cloud schemes. The radiative properties for the Voronoi scheme in shortwave fluxes can be explained by its lower asymmetry factor than the other four schemes, leading to smaller proportion of forward scattering and larger backward scattering and thus less shortwave flux reaching the ground and more upward flux for the Voronoi scheme compared with the other four schemes.

Table 2 TOA upwelling and downwelling surface fluxes (W m⁻²) and TOA SWCF and LWCF (W m⁻²) for the five schemes compared with CERES EBAF products. Numbers in parentheses are differences between the five scheme simulations and CERES EBAF products.

4.3 Climate model simulation results

As shown in RRTMG simulations in Sect. 4.2, the influence of the five ice cloud schemes (Voronoi, Yi, Mitchell, Baum-yang05 and Fu schemes) on the radiative effects is evaluated under standard tropical atmospheric profiles and with assumptions of idealized ice cloud microphysical properties as input data. Simulation results based upon the radiative transfer model are capable of showing differences in ice cloud radiative effects for the five schemes under some specific conditions but are unable to demonstrate comprehensive the performance of the five schemes corresponding to real and complex atmospheric situations. To study the ice cloud modeling capabilities of the five schemes in the climate model, the Voronoi, Yi, Baum-yang05 and Fu schemes are applied in the CIESM. Upwelling and downwelling fluxes and TOA SWCF and LWCF from CIESM simulations for the five schemes compared with CERES EBAF products are plotted in Figs. 7 and 8, and their globally averaged values are listed in Table 2. The results show that the Voronoi scheme produced a lower difference of approximately −0.45 W m⁻² (1.1 %) for the TOA shortwave cloud radiative forcing and −0.30 W m⁻² (1.4 %) for the TOA longwave cloud radiative forcing compared with the four existing schemes. For FSDS, the Voronoi scheme has the smallest downwelling fluxes at the surface and is the closest to the EBAF products due to the Voronoi scheme scatter having the least energy in the forward direction. For FSUTOA, the Voronoi scheme possesses the largest upwelling fluxes compared to the other four schemes due to its strong backward scattering. For the longwave spectrum, the effects of the Voronoi scheme on the FLDS and FLUTOA are negligible.

To discuss the influence of the five schemes on the global distributions of SWCF and LWCF, the zonal average analysis is shown in Fig. 9. Results show that the SWCF of CIESM simulations is stronger than the EBAF results in tropical regions. This might be because the CIESM overestimates ice and liquid cloud fraction in low latitudes compared with the observations (Eidhammer et al., 2014; Kay et al., 2016). Note that the Voronoi scheme exhibits weaker cooling effects in tropical regions than the other four existing schemes. According to the definitions of the SWCF and LWCF, the highest TOA upward fluxes of the Voronoi scheme can produce the lowest TOA net fluxes, which means the Voronoi scheme can produce the smallest negative TOA SWCF among the five schemes. Figure 10 displayed the distribution of differences between the five schemes and the CERES EBAF product. The difference boxes of the Voronoi scheme are the most concentrated on the zero line, and its statistical deviation is the smallest, which means the spatial distribution of cloud radiative effects of the Voronoi scheme is closer to EBAF products compared with the other four existing schemes.

5 Conclusions

The optical property parameterization (Voronoi scheme) of the Voronoi ice crystal scattering (ICS) model is investigated for simulations of the ice cloud radiative properties in the Community Integrated Earth System Model (CIESM). The Voronoi scheme is completed based on the single-scattering properties of the Voronoi ICS database and particle size distributions from in situ observations. The band-averaged optical properties of ice clouds including the mass extinction coefficients, single-scattering albedo and asymmetry factor of the Voronoi scheme are applied in the CIESM climate model and compared with those of the four existing schemes (Baum-yang05, Fu, Yi and Mitchell). The results show that the Voronoi scheme has a distinct feature of having the lowest asymmetry factor in the shortwave bands. This feature could be related to the complex multifaceted shape of the Voronoi ICS model and suggests that the Voronoi scheme can produce relatively stronger backward scattering compared with the other schemes.

Radiative properties of ice clouds are firstly assessed in the rapid radiative transfer model for general circulation models (RRTMG) in the CIESM. The profiles of upward and downward fluxes from different ice cloud schemes are simulated for the prescribed atmospheric condition. The RRTMG results show that the Voronoi scheme has the highest upward flux at the top of the atmosphere (TOA) and lowest downward flux at the surface when the solar zenith angle equals 60^∘. Therefore, the net flux of the Voronoi scheme is largest at the TOA and smallest at the surface compared with the other schemes, which is mostly due to its lowest asymmetry factor.

The five schemes (Baum-yang05, Fu, Yi, Mitchell and Voronoi) are then applied to the CIESM to simulate shortwave and longwave global total cloud radiative forcing at the TOA during 2001–2010. The accuracy of simulated 10-year global total cloud radiative forcing from different schemes is evaluated with the Clouds and the Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) product. The results show that the Voronoi scheme produced a lower difference for the TOA shortwave and longwave cloud radiative forcing compared with the four existing schemes. Especially for the region (from 30^∘ S to 30^∘ N) where the ice clouds occur frequently, the Voronoi scheme provides the closest match with the CERES EBAF product.

In conclusion, simulations of globally averaged shortwave and longwave cloud radiative forcing at the TOA from the five schemes are investigated through the EBAF product. We find that the Voronoi scheme present a better agreement with EBAF products than the other schemes, especially in the tropical region, which confirmed that the Voronoi ICS model has the possibility of having ice cloud modeling capabilities in the climate model of the CIESM.

Appendix A: Nomenclature

L	Particle maximum diameter (µm)
λ	Wavelength (µm)
SZP	Size parameter (unitless)
n(L)	Particle concentration (cm−3)
N0	Intercept coefficient of n(L) (unitless)
λ∗	Slope coefficient of n(L) (unitless)
μ	Dispersion coefficient of n(L) (unitless)
PSD	Particle size distribution defined by N0, λ∗, μ and L
TOA	Top of atmosphere
βe,s,a	Extinction, scattering and absorption coefficients
σe,s,a	Extinction, scattering and absorption cross section
θ	Inclination to the upward normal direction scattering angle
$\mathit{\mu },{\mathit{\mu }}^{\prime }$	Cosines of θ and incoming and outgoing intensity direction, respectively
$\mathit{\varphi },{\mathit{\varphi }}^{\prime }$	Incoming and outgoing intensity azimuthal angle with reference to the x axis, respectively
P	Phase function regulated by $\mathit{\mu },\mathit{\varphi },{\mathit{\mu }}^{\prime }$ and ϕ′
z	Upper limit of the outer boundary
τ	Optical thickness
I	Total (direct plus diffuse) radiance
B[T]	Planck's function
$J(\mathit{\tau };\mathit{\mu };\mathit{\varphi })$	Source function
De	Effective particle diameter
Qext,  sca(λ,L)	Extinction efficiency and scattering efficiency
V(L)	Ice particle volume (µm3)
A(L)	Average geometrical cross section (µm2)
$K_{\text{ext}}\left(\mathit{\lambda }\right){\stackrel{\mathrm{̃}}{K}}_{\text{ext}}$	Spectral and band-averaged mass extinction coefficients
$\mathit{\varpi }\left(\mathit{\lambda }\right)\stackrel{\mathrm{̃}}{\mathit{\varpi }}$	Spectral and band-averaged single-scattering albedo
$g\left(\mathit{\lambda }\right)\stackrel{\mathrm{̃}}{g}$	Spectral and band-averaged asymmetry factor
N	Cloud fraction
Fcloudy	Net fluxes of cloudy conditions
Fclear	Net fluxes of clear conditions
FSDS	Downwelling shortwave flux at the surface
FLDS	Downwelling longwave flux at the surface
FSUTOA	Upwelling shortwave flux at the top of atmosphere
FLUTOA	Upwelling longwave flux at the top of atmosphere
SWCF	Shortwave cloud forcing
LWCF	Longwave cloud forcing