Parametrizations of Liquid and Ice Clouds’ Optical Properties in Operational Numerical Weather Prediction Models

: Parametrization of radiation transfer through clouds is an important factor in the ability of Numerical Weather Prediction models to correctly describe the weather evolution. Here we present a practical parameterization of both liquid droplets and ice optical properties in the longwave and shortwave radiation. An advanced spectral averaging method is used to calculate the extinction coefﬁcient, single scattering albedo, forward scattered fraction and asymmetry factor ( β ext , (cid:36) , f, g), taking into account the nonlinear effects of light attenuation in the spectral averaging. An ensemble of particle size distributions was used for the ice optical properties calculations, which enables the effective size range to be extended up to 570 µ m and thus be applicable for larger hydrometeor categories such as snow, graupel, and rain. The new parameterization was applied both in the COSMO limited-area model and in ICON global model and was evaluated by using the COSMO model to simulate stratiform ice and water clouds. Numerical weather prediction models usually determine the asymmetry factor as a function of effective size. For the ﬁrst time in an operational numerical weather prediction (NWP) model, the asymmetry factor is parametrized as a function of aspect ratio. The method is generalized and is available on-line to be readily applied to any optical properties dataset and spectral intervals of a wide range of radiation transfer models and applications.


Introduction
Clouds have the main influence on the radiative fluxes in the atmosphere and on the entire atmospheric energy budget. The increase of aerosol number concentration can increase cloud droplets number concentration and retardation in their growth [1]. This effect may suppress drizzle and reduce the cloud droplet size, which would cause an increase in cloud albedo (cooling effect) [2,3]. Cloud radiative feedback (CRF), particularly of water clouds, is the dominant uncertainty for climate sensitivity [4,5].
The main goal of this work is to re-calculate the optical properties of both water droplets and ice particles enabling a simple implementation in both numerical weather prediction (NWP) and climate models (CLM). The benefit is measured by the non-hydrostatic compressible COSMO model (http://www.cosmo-model.org). operational models. The common difficulties are the different wavelength intervals used and the high variability in particle-size distributions (PSDs) and shape [6]. Cirrus clouds containing pure ice crystals cover roughly 20-40% of the earth's troposphere and are of great impact on the earth's energy budget [6][7][8]. These clouds have a warming effect at the top of the atmosphere of approximately 5 W/m 2 in the global annual average that was confirmed by numerous modeling and observational studies [9,10]. The optical properties of ice hydrometeors were for many years one of the unsolved problems of climate models. One of the main reasons for this difficulty is the various shapes of ice crystals such as plates, columns, bullet rosettes, and aggregates of these. The size range of ice particles is also much larger compared to water droplets [11]. As to date, there is no general analytical solution for the scattering and absorption of non-spherical ice particles. Nevertheless, ray scattering calculations and other models based on geometric optics can be used to calculate the single scattering optical properties of ice particles with a large range of size parameter in clouds [12,13]. Fu and Liou (1993) [11] followed the work of Takano and Liou (1989) [14] and used a light scattering program including a geometric ray-tracing program to calculate the single-particle optical properties. The mass extinction coefficient, single scattering albedo, and expansion coefficients of the phase function were calculated for 11 observed PSDs (references can be found in Fu and Liou, 1993) [11]. These authors parametrized the optical properties for 18 broadband intervals as a function of ice water content (IWC) and mean effective size, defined as: Where D and L are the randomly oriented hexagonal particles' width and length respectively and n(L) is the PSD. The single-particle aspect ratio D/L used here was based on observational data [15]: Since most atmospheric ice particles are much larger than most of the solar radiation wavelengths, geometric optics approximations are sufficiently accurate to model ice particles and radiation interactions [16].
The generalized effective size (D ge ) is more suitable than D e f f for ice optical properties parametrizations. This quantity is the ratio between the volume and the statistically averaged projected area (A) of the randomly oriented needle particle (see Equations 2.3, 3.10, and 3.11 in Fu, 1996) [11]: The generalized effective size (D ge ) can be related to the effective radius using the formula: The size parameter (α) is defined by the particle size to wavelength ratio: where λ is the wavelength, r va = 3V/3A is the equivalent radius for a non-spherical particle by conserving both its volume (V) and projected area (A) [17]. Ice crystals have a large variety of shapes, sizes, and complexity, which depends mainly on temperature and humidity. The volume extinction coefficient and single scattering albedo for randomly oriented hexagonal ice columns can be parameterized by the generalized effective size and the IWC [15,17]. These parametrizations can be well applied to other ice particle shapes since the IWC, along with D ge , conserves both the projected area and volume of non-spherical ice particles which largely determine the ice particle extinction and absorption [15]. As for the asymmetry factor, ray scattering simulations for many types of ice crystals showed that it is less dependent on size rather on the aspect ratio D/L [16,18] or more generally, the effective aspect ratio [18]: Fu (2007) [18] also showed that the asymmetry factor bulk parameterization based on the aspect ratio (AR) is less sensitive to the detailed ice shape. Fu (1996) [15] calculated the optical properties of cirrus clouds for 28 PSDs. His work was based on both conventional geometric ray-tracing for large particles with a size parameter greater than 200 and improved geometric ray tracing by Yang and Liou (1995) [19] for smaller particles. These PSDs were taken from several in situ aircraft observation campaigns in mid-latitudes and tropical regions. The optical properties were parameterized as a function of D ge and IWC, which for hexagonal particles is: where ρ ice is the bulk ice density (0.917 g·cm −3 ). The parametrizations for volume extinction coefficient and single-scattering albedo were adequate for various randomly-oriented nonspherical particles although the randomly-oriented hexagonal particles were assumed. Note that horizontally oriented ice particles appear infrequently:~12% in mid-latitude and near zero in the tropics (see Table 3 in Thorsen and Fu 2015 [20]). The parameterization of the asymmetry factor as a function of D ge is due to the dependence of the aspect ratio on the ice particle size used (i.e., Equation (2)), which might not be generally valid [18]. In infrared (IR) wavelengths, an ice particle's absorption is the dominant effect but one cannot neglect the scattering effect [21]. Fu et al. (1998) [17] used a linear combination of Mie theory, anomalous diffraction theory (ADT), and the geometric optics method (GOM) to calculate the IR radiative properties of ice crystals and to address the wide range of size parameter. The finite-difference time-domain (FDTD) method [22] was used as a reference calculation. These authors used the same 28 PSDs as used by Fu (1996) [15] dividing the distributions into 38 discrete size bins in the range of L min = 1 µm to L max = 3 mm. The bulk optical properties were parameterized as a function of IWC and D ge for wavelengths in the range between 3.969 µm and 100 µm. Compared to the reference results, the averaged relative error was only~2% for the emissivity of cirrus clouds. Yang et al. (2000) [7] published their work on the single scattering optical properties of six ice crystals habits that are usually found in in situ measurements (plate, columns, hollow column, planar and spatial bullet rosettes, and aggregates) and spectral bands in the range 0.2 µm-5 µm. They used GOM and FDTD to calculate the extinction and absorption coefficients, and the asymmetry factor and forward peak of the phase function (the latter is only significant for short wave intervals), which were parameterized as a function of effective size. In their study, they fitted the diameter of equivalent area D A and the diameter of equivalent volume D v in terms of the maximum dimension D max for each crystal shape where the D max /D eff ratio depends on the ice particle habit. Baum et al. (2005) [23] have investigated the bulk properties of ice particles using PSDs data taken from in situ field campaigns in five different locations over the tropics and  mid-latitudes namely FIRE-I, FIRE-II, ARM, TRMM, and CRYSTAL FACE field campaigns  (see their Table 1 therein). There were 1117 PSDs parameterized by fitting to a gamma distribution of the form: The bulk characteristics of IWC and the median mass diameter D m were represented analytically as a function of µ Λ and N 0 . To overcome the differences between calculated and in situ measured IWC for all field experiments, the size range was divided into parts and for each, a mixture of habits (droxtals, bullet rosettes, solid columns, and aggregates) was assumed. Baum et al. (2011) [24] calculated the bulk optical properties of ice cloud in the shortwave regime based on 12 thousand PSDs from both northern and southern hemisphere field campaigns containing nine habits with different roughness levels. From the comparison between their models and cloud-aerosol lidar with orthogonal polarization (CALIOP) data, they concluded that rough surface particles are best for modeling a mixture of habits that usually exist in frozen clouds. The wavelengths used were in the range of 0.4-2.24 µm in a later work [25], they expanded their results for 445 discrete wavelengths in the range 0.2-100 µm using the Yang et al. (2013) [12] data set. The latter used the Amsterdam discrete dipole approximation (ADDA), the T-matrix, and the improved geometric optics (IGOM) methods in the range of 0.4-100 µm. They included 11 ice particle habits in the size range of 2 µm to 10 mm divided into 189 size bins with three roughness levels. This publication provided the most up to date data library for the widest range of shapes, sizes, and wavelengths.
As the main goal of this work was to implement an efficient, simple, and fast scheme for NWP and CLM models, the single-particle data by Fu (1996) [15] and Fu et al. (1998) [17] was chosen as the base of our parametrizations for ice optical properties.

Water Droplets
The effective radius of water droplets is defined as the ratio of the third to the second moments of a droplet size distribution [26]: There are large uncertainties in simulated cloud cover, liquid water content (LWC), droplets size distributions, and r e from the NWP and CLM [27]. On top of these are the uncertainties in the optical properties of these hydrometeors namely the volume extinction coefficient (β ext ), asymmetry factor (g), and single scattering albedo ( ) at each cloud layer. Fortunately, LWC and r e are usually enough for a reasonable approximation of the cloud's optical properties even without knowing the exact particle size distributions (PSD) in the cloud [28]. Based on the Delta-Edington approximation for vertically inhomogeneous atmosphere [29], Slingo and Schrecker (1982, SS82 hereafter) [30] calculated a 24-band parametrization for g, , and β ext as a function of LWC and r e in the shortwave range of 0.25-4 µm. Slingo (1989) [31] used SS82 parametrizations and further developed the equations for the transmissivity, reflectivity, and absorptivity. These parametrizations were implemented into Helios, the UK Met Office global circulation model (GCM) radiation subroutine that was later replaced by the Edwards and Slingo (1996) [32] scheme. This author also made comparisons with Stephens (1984) [33] parametrization vs. aircraft observations and showed positive results. Hu and Stamnes (1993) (hereafter HS93) [28] published their work on liquid cloud optical properties based on Mie calculations. They expanded the spectral range to the longwave (LW) compared to previous works which mainly addressed the short wave (SW) of the solar radiation. Their work included computation of g, and β ext for 74 wavelengths between 0.29-150 µm and parametrized them for three droplet size ranges: 2.5-12 µm, 12.5-30 µm, and 30-60 µm. The parametrization was defined by non-linear square fitting to the exponential functions: where a i , c i are the fitting coefficients. HS93 also showed that for a fixed effective radius, the optical properties change only slightly (1%-2%) even under large variations in the PSD width, shape, skewness, and modality (number of modes in the distribution). Even changing the PSD from lognormal to a generalized gamma distribution with different widths showed only up to 8% differences in the optical properties of the clouds. The parameterization used was as accurate as exact Mie calculations but three orders of magnitudes faster. Later, Lindner and Li (2000) [34] calculated cloud droplets' optical properties for 36 wavelengths in the infrared spectrum between 4 µm and 1000 µm. Based on the exact Mie calculation, they used the least square fitting technique for the bulk optical properties for particles in the size range of 2-40 µm. While HS93 fitted the optical properties to exponential functions and used three size bins, Lindner and Li (2000) [34] defined a single parametrization for the whole size range and used rational functions of the form: Since HS93 includes both SW and LW spectral ranges, we adopted the HS93 fitting functions but used a single parametrization for the whole size range.

Ice Particles
Raw data for the single-particle extinction and scattering cross-sections (σ ext , σ sca ), asymmetry factor (g p ), and the forward scattering δ function peak (f δ ) were taken from Fu (1996) [15] and Fu et al. (1998) [17]. In this work, all particles were assumed randomly oriented hexagons particles of length L and diameter D, which are related in Equation (2). Although a more recent state-of-the-art work on a wider range of habits was available, most notably by Yang et al. (2013) [12], there were two main reasons for our choice. First, the parametrizations of the optical properties in terms of IWC, generalized effective sizes, and mean aspect ratio could be applied to ice clouds containing ice particles with complicated shapes [15,17,18] by conserving ice particle ice volume, projected area, and mean aspect ratio. For example, as shown in Fu (2007) [18], for the SW, this parameterization well reproduced the asymmetry factors of complicated ice particles such as bullet rosettes, aggregates with rough surfaces, and fractal crystals and agreed well with observations. Second, the diagnosis of the habit-mixture composition of ice particles at each grid-box containing ice was not straightforward while the computational cost could be high. Although scientifically reasoned, it was impractical for present-day operational numerical weather prediction models.
Further discussion on the bulk optical properties of cirrus clouds, shape consideration, and their effects on climate models can be found in Baran et al. (2014Baran et al. ( , 2016 [35,36]. For SW radiation, L and D were discretized in 30 non-equidistant sizes from 8.75 to 3100 µm [15]. The data covered 200 wavelengths from 0.21 to 5 µm, which we used for SW. For the thermal infrared range (LW) the particle axis ratios were the same, but the data set was extended by 10 smaller particle sizes down to L = 2 µm. There were 36 wavelengths from 3.969 µm to 100 µm [17].
The equations for the transition from single-particle data to cloud optical properties were: where β ext , β sca are the volume extinction and scattering coefficients, ω is the single scattering albedo, g is the asymmetry factor, and f d is the delta fraction (g p and f δ are the single-particle asymmetry factor and delta fraction respectively). β ext is the mass extinction coefficient. An ensemble of particles was defined with a PSD, which was assumed to follow a four-parameter modified gamma distribution MGD [37]: where N 0 is the scale parameter and µ, υ, Λ control the shape of the distribution, M i is the size distribution ith moment, and L is the mean particle length. A simpler gamma distribution was obtained by setting υ = 1 hence: We composed our PSD ensemble by a systematic variation of the PSD parameters within a wide range. We changed L in the range of 5-3000 µm in 500 logarithmically equidistant steps, while µ was changed from 0 to 14 in steps of 1. Λ was computed from L and µ using (2.10) resulting in 7500 PSDs. Applying the described ensemble of PSDs into Equation (3) with D(L) from Equation (2), resulted in a range of 5-600 µm in the generalized effective size (D ge ). The aspect ratio AR was computed using Formula 5 and results were in the range of 0.22-1.0. In the NWP model (i.e., COSMO), D ge and AR could be evaluated from the mass-size relation used by the cloud microphysics scheme, for example, m = αL b where m is the particle mass. See Smith et al. (2016) [38] for further discussion on the issue.
The technique of spectral averaging over wide bands, representative for typical ice water paths (IWP) encountered in NWP-models, depends on the parameter one desires to compute. For β ext , the physically appropriate definition equation for spectral averaging would be as follows: where λ is the wavelength, S(λ) is the Planck's function for an appropriate radiation source temperature, P(∆z) is a certain probability density function (PDF) of the vertical model layer thickness ∆z, f is the IWC PDF, and IWC ∆z is the ice water path (IWP). Nevertheless, Equation (26) was not practical because of the unclear definition of the "mean" ice water path IWC∆z. Instead, we propose to first separately average over wavelength and then average the result β ext λ over IWC and ∆z. The first step is: For numerical quadrature (integration), we employed the trapezoidal rule on the non-equidistant original wavelength discretization of the single-particle scattering data. The second step based on this result is: For quadrature, we used an equidistant discretization of IWC and ∆z with constant increments ∆IWC and ∆(∆z) and also employed the trapezoidal rule: The IWC PDF f was defined as a triangular function in the interval [0, IWC max ]: where: The level thickness function was defined as a uniform distribution in the range [∆Z min . ∆Z max ]: where: For our purposes IWC max = 1 gm −3 , ∆IWC = 0.01 gm −3 , ∆Z min = 100 m, ∆Z max = 600 m and ∆(∆z) = 10 m were chosen. For the single scattering albedo calculation we used the same numerical method as in (27) but the following formula: The same is done for the asymmetry factor using the formula: The formula for the forward peak f d is similar to (34). We parametrized β ext , ω, g, f d as a function of D ge and AR by fitting rational functions division using the Padé approximation form: where x can be D ge or AR and a i , b i are fitting coefficients. This approximation is a common practice for these parametrizations [39] and requires a non-linear fitting algorithm like Levenberg-Marquardt. It is very flexible and can reproduce nearly any asymptotic behavior if N and M are chosen adequately. If N = M, an asymptotic convergence to a constant for x → ∞ results. If N < M, the function converges to 0. We tried many other fitting functions, also the ones available in previous publications, but the mentioned method led to the best results for all spectral intervals. For β ext D ge , we chose N = 3; M = 4, for ω D ge and f d (AR), we defined N = 3; M = 3 while for g(AR) we put N = 3; M = 2. The reasoning for choosing the AR and not D ge as the fitting parameter for both the asymmetry factor and delta function is thoroughly discussed by Fu (2007) [18]. For simplicity reasons we used g(AR) also for the LW intervals since scattering is less important in the thermal spectrum and we assumed the radiation model would not be sensitive to the difference between g(AR) and g D ge . Another reason for choosing g(AR) was the lack of confidence in the extrapolation of the fitting functions of g R e f f for particles larger than the largest calculated point (R e f f = 370 µm or D ge = 570 µm) in the LW regime. Unlike the asymmetry factor, both the extinction coefficient and the single scattering albedo fits can be safely extrapolated to particles larger than that in the entire spectral range. These parametrization fits were applied to the eight spectral bands of the COSMO model radiation scheme (Ritter and Geleyn, 1992-RG92 thereafter) [40]. An example of the parametrization fit for three spectral intervals: 0.25-0.7 µm (SW), 1.53-4.64 µm (SW) and 20-104 µm (LW) can be seen in Figure 1. We could see β ext had an exponential decrease as a function of size in the SW, while in the LW intervals, a peak was seen around an effective size of 30 µm. A single scattering albedo of 1 for the entire range of particle sizes in the 0.25-0.7 µm (SW) interval meant that 100% of extinction was due to scattering. For longer wavelengths, this amount was decreasing adding more absorption contribution to the extinction. A similar peak in single scattering albedo was seen in the LW as was seen in the β ext curve. As for the asymmetry factor, we could see an almost linear decrease as a function AR which was a surprisingly simple description for an ensemble of 7500 PDFs and a very large range of particle sizes. Contours of the optical properties parametrizations can be seen in Figure 2. The extinction efficiency in this plot is calculated using the relation: The asymmetry factor presented in Figures 1 and 2 were calculated for particles with smooth surfaces, where the forward scattered peak (delta fraction f d ) was large because of the plane-parallel surfaces of the ice particles. However, real particles might have rough surfaces, where there is more diffuse scattering. Therefore, an asymmetry factor for rough surfaces was defined as [18]: The results for f d and g r for the 0.25-0.7 µm and 1.53-4.64 µm spectral intervals are shown in Figure 3. We note here that g r was calculated following the same methodology applied for g but using Equation (37). Although the differences between g and g r were less than + 3% for the spectral intervals selected in this work, a small but important difference will be presented in the results section. We assumed here that the scattering in the LW was less important (although not negligible) and that uncertainties in its asymmetry factor would not significantly affect the radiative energy budget.  In the original RG92-formulation, the scattering phase function is approximated by the so-called δ-Eddington approximation, which is a sum of a general two-term Legendrepolynomial expansion of the phase function plus a delta-forward peak, weighted by the forward scattered fraction f [41]. This approximation was "tuned" to the original phase function by matching g, which is the first moment of the phase function, and by requiring that f = g 2 , which follows from the matching of the second moment of the phase function with that of the Greenstein function.
For ice particles in the visible spectral range, we changed the RG92 formulation in a way that f was an independent free parameter [18]: For rough surfaces particles and: For smooth surfaces, for which a directly transmitted fraction f d is added (planeparallel surfaces). We did, however, limit f so that it could not be larger than g 2 (AR) .

Liquid Droplets
For the quasi-spherical water droplets parametrization, we defined the effective radius as: This definition is motivated by the large-size limit where β ext is directly proportional to R −1 e : lim Where A geo the geometric cross-section of a spherical droplet, LWC is the liquid water content, n w , ρ w are the number concentration and bulk density of water, respectively. The constant coefficients c 1 , c 2 depend on the constant mass-size-relation parameters and the size distribution parameters.
Based on Mie calculations performed by HS93 with a few typing corrections in the original tables, and using the spectral averaging technique described in the previous section, we applied the calculation to the 8 RG92 wavelength bands. While for ice particles we had the raw single-particle data for certain wavelengths, here we employed the original fit functions as a function of R e for discrete equidistant R e values. The same choices for ρ w,max = 1 gm −3 , ∆z min = 100 m and ∆z max = 600 m as was chosen for ice, were used here. We used the same fitting ansatz as for ice using formula (34) but chose N = 3; M = 4, for β ext (R e ) while for ω(R e ) and g(R e ) we defined N = 3; M = 3. The three size sub-ranges of fitting used by HS93 were replaced with one formula that fit all ranges and was valid in the range 2.5 µm → ∞. An example of the parametrization fit can be seen in Figure 4. A full map of the optical properties in all sizes and wavelengths can be seen in Figure 5.

Model Description
The optical properties parametrizations obtained in Section 2 were implemented in two NWP models COSMO [42][43][44][45][46] and ICON [47,48]. In the next section, the sensitivity tests are performed using the COSMO model. Hence, a more detailed description of this model is provided below. COSMO is based on the primitive thermo-hydrodynamic equations that describe non-hydrostatic compressible flow in a moist atmosphere. The model uses a rotated geographical Arakawa C-grid and generalized terrain-following height coordinates with a Lorenz vertical staggering method. Its vertical extension reaches 23.5 km (about 30 hPa) with 60 model levels in the atmosphere. It uses a two-time-level integration scheme based on the Runge-Kutta method of third order for the prediction of the three Cartesian wind components u, v, and w, the pressure perturbation p from a hydrostatic base state and the temperature perturbation T . A fifth-order upwind scheme is used for their horizontal advection. The default version of the COSMO model uses 1-moment bulk microphysical parameterization. In addition to the mass fraction of water vapor, the microphysical prognostic variables are liquid water contents of cloud water and rain, and ice water contents of cloud ice, snow, and graupel. As shown in Section 2, the parametrizations of the optical properties require the knowledge of the effective radius of cloud water and rain r e , and of the generalized effective size of cloud ice, snow, and graupel D ge . r e in (2.28) and D ge [17,49] can be approximated using the ratio of water content to the number concentration using the form R e = c 1 LWC n w c 2 where c 1,2 depends on the constant mass-size-relation parameters and the size distribution parameters. Some parametrizations of D ge considers both temperature and IWC for D ge . However, the number concentration of the hydrometeors is not predicted by the model and needs to be diagnosed. In our new implementation, the number concentration of liquid hydrometeors was a tuning parameter, while the number concentration of solid hydrometeors was parametrized as a function of temperature (following Cooper, 1986 [50] for cloud ice and Field et al., 2005 [51] for snow). The default COSMO model RG92 optical properties parametrizations are based on Mie calculations by Stephens (1989) [52] and the parametrizations proposed by Slingo and Schrecker (1982) [30] for water clouds. For ice hydrometeors, only a parametrization based on IWC only is implemented (COSMO internal documentation). In the following section, we analyze the effect of the new optical properties parametrizations on the radiation transfer and demonstrate the important role of particle sizes r e and D ge .

Sensitivity Analysis
In this section, we examine the impact of the new optical properties parametrizations on the radiation transfer and compare it with the default COSMO RG92 scheme under different cloudiness situations. For this purpose, we have utilized the idealized version of the COSMO model [53] to simulate radiation transfer through two types of clouds: warm stratus and cirrus. The characteristics of the simulated clouds are summarized in Table 1. In order to analyze the effect of clouds' width and microphysical properties on the radiation transfer, seven simulations of warm stratus and three simulations of cirrus were performed. These simulations were intended to test the radiation transfer parametrizations only, and not to mimic real weather events or observation campaigns. However, the variety of simulated conditions (Table 1), together with various pollution conditions (described below), properly represented the expected effect of the new parametrization in real weather events. The implementation of new parametrizations generally influenced the shortwave downwelling solar radiation (or global radiation hereafter) both directly and indirectly. The change in radiative heating lead to a different evolution of the simulated cloud, which in turn influenced the radiation transfer. For simplification, we ignored this feedback mechanism and removed the radiative heating term in the equation for temperature tendency. All the idealized COSMO simulations were performed for 6 hours while the solar zenith angle was fixed to 30 • . The domain was chosen to be 100 × 100 km 2 with a flat and uniform lower boundary surface and periodic boundary conditions and a horizontal resolution of 0.025 • (~2.5 km). Both deep and shallow convection Tiedtke-type schemes were switched off. Each simulation was initiated with vertical profiles of temperature and humidity suitable for stratiform cloud formation at the chosen heights (Table 1). During the run period, the simulated clouds slightly changed due to evaporation, particle sedimentation, and other processes. However, they continuously covered the entire simulation domain homogeneously. For each simulation, the output fields (global radiation, total water path, etc.) at 5 min steps were averaged over the domain. Ignoring the first 50 min stabilization period, we obtained 51 values for each field for each simulation. Figure 6a compares the global radiation of RG92 and new parametrization for the seven experiments of warm stratus (see Table 1). The global radiation is shown versus the liquid water path (which is combined from seven experiments, 51 points each). Obviously, the global radiation was smaller for thicker clouds. One could see that in the new version, the clouds were optically thinner than in RG92. The global radiation difference was larger for thicker clouds, reaching 20%.  Table 1) marked by the dashed circle. The red dots from left to right correspond to cloud droplets number concentrations of 1000, 500, 200, 100, 50, 10 cm −3 . The range of variation in global radiation due to change in effective radius is marked by dashed lines.
As explained in Section 3, the new scheme parametrized the effect of cloud water number concentration on the effective radius, which influenced the optical properties and the global radiation. Thus, for the warm stratus of experiment 5 (see Table 1), six more runs were performed with various number concentrations. Figure 6b shows the effect of cloud water number concentration on global radiation in the new scheme. In contrast to RG92 where the effective radius was not dependent on number concentration, the new scheme showed a significant effect of number concentration and effective radius on the global radiation. Having the same liquid water path, a pristine cloud (number concentration of 10 cm −3 ) may double the global radiation compared to a polluted one (number concentration of 1000 cm −3 ). The immense effect on the cloud optical depth due to cloud number concentration variations was thoroughly discussed by Lu and Seinfeld (2005) [54].
Similar to the warm stratus presented in Figure 6a, Figure 7a shows the global radiation versus ice water path for the three experiments of cirrus (see Table 1). In addition to the default RG92 (blue dots), four more versions of the new parametrization are presented (see Section 2.1 for details): the first assumes rough ice particles surfaces (red dots), the second is similar to the first but adding snow particles into consideration (lime green), the third is similar to the second but assuming smooth surfaces (magenta), and the last is similar to the previous but using the forward scattered fraction approximation of f = g 2 instead of the Fu (2007) [18] formulas 37 and 38 (cyan). Obviously, the global radiation was smaller for thicker clouds. In contrast to the warm stratus, the cirrus cloud in the new versions was optically thicker than in RG92. The largest difference from RG92 occurred when assuming rough ice particles and considering snow effects on radiation, where the effect on the global radiation reached 15% (for thick cirrus-experiment 3). Note that considering snow significantly reduced the global radiation (from red to lime green). Note also that the assumption of smooth surfaces made the clouds optically thinner (from lime green to magenta and cyan), and that there was no significant effect of the different assumptions for the forward scattered fraction f (cyan and magenta). In blue dots the default COSMO radiation scheme RG92, red dots represent simulation which assumes rough ice particles surfaces and ignoring snow particles, lime green dots is the result when snow particles are considered, the same but assuming smooth surfaces are in magenta and the last, in cyan, is similar to the previous but using the forward scattered fraction approximation of f = g 2 instead of the Fu (2007) [18] formulas. (b) Global radiation in the third version of the new scheme (red dots) versus cloud-averaged ice particles effective diameter for experiment 2 of cirrus (see Table 1) marked by the dashed circle. The range of variation in global radiation due to change in ice particles effective diameter is marked by dashed lines.
As explained in Section 3, the new scheme parametrized the influence of ice effective diameter D ge on the optical properties and the global radiation. As mentioned, D ge itself was calculated in the model from the ice water content and the number concentration that was a function of temperature. In order to highlight the sensitivity of global radiation to D ge , we implemented an option to artificially determine the constant D ge value at all cloudy grid points. Thus, for the cirrus of experiment 2 (see Table 1), six more runs were performed with various D ge values between 5 and 200 microns. The new version with the assumption of rough needles ignoring the snow effect (red dots in Figure 7a) was chosen for these runs. Figure 7b shows the influence of D ge on the global radiation in the new scheme. In contrast to RG92 where D ge was not introduced, the new scheme showed a significant effect of D ge on global radiation. Having the same ice water path, a cirrus cloud with large particles (an effective diameter of 200 microns) could yield 10% higher global radiation than an identical cloud with small particles (effective diameter of 5 microns).
The important quantity for climate applications is the net heating/cooling effect of the clouds. Usually, during the daytime, the SW incoming radiation heats the atmosphere, which is then cooled down via LW outgoing radiation during night time. We cannot discuss this quantity here since our 6-hour long simulations assumed a fixed solar zenith angle of 30 • . However, we can check the effect of the new scheme on the net radiation at the top-of-atmosphere (incoming minus outgoing). Figure 8 presents the net top-of-atmosphere radiation difference between the new scheme and RG92 versus liquid water path for seven experiments of warm stratus (Figure 8a) versus liquid water path for three experiments of cirrus (Figure 8b). Positive (negative) values indicated the warming (cooling) effect by the new scheme with respect to RG92. In Figure 8b, the differences (with respect to RG92) were calculated for the same four versions of the new scheme as described in Figure 7. One could see that the new scheme had a warming effect with respect to RG92 for warm stratus, and a cooling effect for cirrus. Generally, these effects were enhanced (in absolute values) for optically thicker clouds. Note that the assumptions regarding the roughness of the ice particle surfaces were of similar importance to the snow for the net radiation (compare cyan, magenta, and green dots in Figure 8b). It was nicely seen that the consideration of snow enhanced the net cooling, which could be explained by an increase in the cloud albedo.

Real Cases Evaluation
In this section, we compare the default and the new schemes using numerical weather prediction (NWP) simulations with fully interactive radiation. We utilized the convectionpermitting configuration of the COSMO model over the eastern Mediterranean (domain 25-39 E/26-36 N) with 2.5 km horizontal resolution and 60 vertical levels reaching 23.5 km. The simulations were initialized and driven by the Integrated Forecast System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF). For each of the COSMO versions, we performed 5 months (July-November 2019) of 78 hours simulations twice a day. Hourly forecasted global radiation was verified against 17 automatic radiation stations over Israel in cloudy conditions. Only cloudy events with observed and simulated cloud cover above 95% were selected for global radiation verification. For observations, the cloudiness products from CMSAF (https://www.cmsaf.eu) were used. The verification was performed for low and high cloudiness situations, separately, as defined by the CMSAF cloud type algorithm. We have assumed that during the verification period high clouds are composed mainly of solid particles and low clouds from liquid droplets. As proof, we have used the simulations to retrieve the temperatures at the cloudy layers which correspond to the observed clouds. Indeed, the temperatures in high cloudiness situations were below −19 • C, while the temperatures in low cloudiness situations were above 6 • C. Figure 9 shows the results of the separate verification for low cloudiness (Figure 9a,b) and high cloudiness (Figure 9c,d). Figure 9a,c show the distribution of the error (bias) having a bin width of 50 W/m 2 . A positive error meant model overestimation. Figure 9b,d show the mean diurnal cycle of the global radiation in both observations (black), the default COSMO RG92 version (blue), and the new version (red). In the idealized simulations presented in Section 4, the global radiation by the two COSMO versions was compared in identical cloudiness conditions. In real case simulations with cloudy conditions, having the same 3D cloud cover, water content and particle effective radius in various model versions and reality are hardly achievable. Therefore, the verification conclusions regarding radiation transfer through clouds may be drawn only statistically based on many events, yet assuming that there are no biases in the forecasts of cloudiness depth. In order to reduce the noise originating in comparing the global radiation in cases of significantly different cloud depths, we ignored the events where the global radiation absolute error was above 200 W/m 2 due to cloud thickness differences (see Figure 9). Panels a and b show that the new version, having the updated clouds optical properties, reduced the mean bias from −7.1 to + 3 W/m 2 . In agreement with Figure 6, the new version resulted in slightly higher global radiation. However, the improvement was not statistically significant. The improvement was much more evident for high cloudiness, where the new version reduced the global radiation overestimation from 15.5 to 1.2 W/m 2 . The reduction in the global radiation by the new version was as expected from idealized COSMO runs (Figure 7) not only due to the change in ice optical properties parametrization but also due to the inclusion of snow in the radiation transfer scheme.

Conclusions
In this study, we re-calculate liquid water droplets and ice particles' optical properties in clouds using former works available in the literature regarding these hydrometeors. For water droplets, a unified and general mathematical formulation was used for both the short wave and long wave regimes. For ice particles, thousands of generalized gamma distributions were assumed to simulate a broad range of ice clouds. An advanced spectral integration method was introduced offering more practical use in NWP and climate models.
Sensitivity analysis using COSMO ideal framework to simulate stratus water clouds showed a small impact of optical properties parametrization on clouds optical thickness and as expected, the 5 months of real COSMO simulations have indicated a relatively small improvement of the new scheme compared to the default parameterization. Nevertheless, a much greater effect on the results was obtained for variation in the effective radius. This emphasizes the advantage of implementing a realistic cloud number concentration calculation. A possible solution could be to use prognostic aerosols in cloud formation schemes. The same argument holds for the ice particles' effective size and the need for precise simulations of nucleation processes in cirrus clouds. As for ice clouds, we showed a much greater impact of the new ice particles' optical properties parametrizations on radiation. We also showed that both the inclusion of snow particles and the smooth/rough surface assumptions can have a significant effect on the net surface radiation of a few tens of W/m 2 . The calculation of the energy balance of the Earth is strongly dependent on ice and water optical properties and the hydrometeor's effective size. The large (>50 W/m 2 ) uncertainties can be critical in any climate model but maybe masked out due to the partial cloud cover of the Earth or by cancelation effects such as between liquid and ice particles as demonstrated in our work. A web application for ice and liquid hydrometeors optical properties parametrizations based on this work was created and available for the public on the COSMO webpage (www.cosmo-model.org/content/tasks/operational/ims/ cloudOptics/default.htm).