Bulk hydrometeor optical properties for microwave and sub-mm radiative transfer in RTTOV-SCATT v13.0

Abstract. Satellite observations of radiation in the microwave and sub-mm spectral regions (broadly from 1 to 1000 GHz) can have strong sensitivity to cloud and precipitation particles in the atmosphere. These particles (known as hydrometeors) scatter, absorb and emit radiation according to their mass, composition, shape, internal structure, and orientation. Hence, microwave and sub-mm observations have applications including weather forecasting, geophysical retrievals and model validation. To simulate these observations requires a scattering-capable radiative transfer model and an estimate of the bulk optical properties of the hydrometeors. This article describes the module used to integrate single-particle optical properties over a particle size distribution (PSD) to provide bulk optical properties for the Radiative Transfer for TOVS microwave and sub-mm scattering code, RTTOV-SCATT, a widely-used fast model. Bulk optical properties can be derived from a range of particle models including Mie spheres (liquid and frozen) and non-spherical ice habits from the Liu and Atmospheric Radiative Transfer Simulator (ARTS) databases, which include pristine crystals, aggregates and hail. The effects of different PSD and particle options on simulated brightness temperatures are explored, based on an analytical two-stream solution for a homogeneous cloud slab. The hydrometeor scattering "spectrum" below 1000 GHz is described, along with its sensitivities to particle composition (liquid or ice), size and shape. The optical behaviour of frozen particles changes in the frequencies above 200 GHz, moving towards an optically thick and emission-dominated regime more familiar from the infrared. This region is previously little explored but will soon be covered by the Ice Cloud Imager (ICI).


σ e (D g )n g (D g ) dD g . (1) The integration is done over a size range D min to D max that will be discussed in Sect. 3.2. Note that in this work the prime on the PSD notation n g (D g ) indicates that it has been rescaled to account for numerical integration issues and the limited size range, a process referred to as renormalisation (Sec. 3.2.2, Eq. 17).
To represent scattering in active and passive microwave/sub-mm radiative transfer in a fast model like RTTOV-SCATT 95 requires also the bulk scattering and backscatter coefficients β s and β b , both in [m −1 ], and the dimensionless bulk asymmetry parameter g which summarises the mean direction of scattering (strictly, g is the phase-function weighted mean of the cosine of the scattering angle; see Petty, 2006, for full definitions). These bulk properties are computed from the single-particle scattering and backscatter cross-section σ s (D g ) and σ b (D g ) [m 2 ] and single-particle asymmetry g single (D g ) [dimensionless], again by integrating over the PSD: 100 β s = Dmax Dmin σ s (D g )n g (D g ) dD g ; (2) g = 1 β s Dmax Dmin g single (D g )σ s (D g )n g (D g ) dD g . the bulk extinction coefficient β e . In an exception to the SI policy used elsewhere, the units of the extinction coefficient are [km −1 ]; the dimensionless bulk single scattering albedo (SSA) ω 0 = β s /β e ; the dimensionless bulk asymmetry parameter g; if the targeted sensor is a radar, also the bulk radar reflectivity factor Z = (10 18 /z 0 )β b , in [mm 6 m −3 ]. See Appendix A 115 for full definition.
The data files contain bulk optical properties for a set of possible hydrometeor types. The default configuration of the table generator is given in Table 1, which provides 5 hydrometeor types representing rain, "snow" (referring to precipitating particles in stratiform cloud), "graupel" (referring to all ice particles in convective cores), cloud water and cloud ice (referring to suspended frozen particles). This maps onto typical hydrometeor representations in global forecast models. However the total 120 number of hydrometeor types in RTTOV-SCATT is unlimited and this could be used to build up more complex representations (for example, there could be different hydrometeor types for tropical and extratropical ice cloud). Each hydrometeor type is defined by a set of physical options, with the main options illustrated in Table 1; these will be described in more detail in the rest of this article. The default settings for frozen hydrometeors were obtained from a multi-dimensional parameter search in order to produce the best fits between ECMWF modelled brightness temperatures and Special Sensor Microwave Imager 125 Sounder (SSMIS) observations (Geer, 2021b). The settings for rain and cloud water have been inherited from Bauer (2001).
Each hydrometeor type needs to be associated with one of the "placeholder" types listed in Table 2. This gives the hydrometeor a descriptive name and indicates whether it is frozen or liquid. It also associates a density and size range D min to D max , which are mainly relevant when the Mie sphere approximation is used to compute the optical properties. If optical properties are taken from a scattering database, then the density and integration size range are instead implied by the settings for that 130 database. Any placeholder type can be associated with any hydrometeor from the databases. The number of hydrometeor types in the lookup tables is unlimited, so the placeholder types can be re-used multiple times. Their names should not be taken too literally, such as in the use of the "snow" placeholder type with large-scale snow, and the "graupel" placeholder for frozen particles in convective clouds. The "totalice" placeholder was designed for the Met Office global model, which represents all ice particles as a single species (e.g. Doherty et al., 2007). 135 To save time, the table generator first calculates optical properties for the full set of possible channels (currently 136) and then it composes these channels into sets corresponding to satellite instruments (where currently 34 are represented). The list of channels and instruments is defined in a configuration file, along with the chosen hydrometeor types and physical assumptions as summarised in Table 1. One output file is then produced per instrument, containing the optical properties for each channel, tabulated as a function of water content, temperature and hydrometeor type. The relevant grids are fixed at:
-For temperature, 70 points, one every Kelvin, from 204 K to 273 K for frozen hydrometeors and from 234 K to 303 K for liquid hydrometeors (frozen or liquid is defined according to Table 2).
The channels for any target instrument can be specified in one of two ways. The normal way is the "condensed" approach, designed for unpolarised optical properties. Instruments such as conical-scanning microwave imagers may separately measure 145 vertical and horizontal polarisations at one frequency, but in the condensed representation this is represented as one channel in the tables. This eliminates redundancy in the data file and leaves it up to the radiative transfer model to remap the condensed set of channels onto the actual channel set of the instrument. An alternative "full" option is for use with polarised optical properties (Sec. 3.3) and in this case, the full channel set, including each distinct polarisation, is represented in the data file.
Given that over the width of a channel, bulk scattering properties vary slowly with frequency, there is no attempt to represent eteors represented in the hydrotables. The hydrometeor bulk optical properties are retrieved from the hydrotables as a function of temperature, in-cloud water content, and channel. These are then summed over the set of hydrometeors, together with the gas optical properties driven mainly by oxygen and water vapour absorption, to provide the total bulk optical properties of each layer in the model (see e.g. Bauer, 2001). These profiles are then input to the solver for scattering radiative transfer, which in the case of RTTOV-SCATT uses a delta-Eddington approach (for further details, see Bauer et al., 2006). Sub-grid heterogeneity 160 of cloud fields is represented through an effective cloud fraction (Geer et al., 2009).
2.1 Single-particle optical properties Single-particle optical properties are derived using either Mie theory, which assumes spherical particles, or from scattering databases that summarise the properties of non-spherical particle habits, which have been computed using more sophisticated methods such as DDA. Currently, the Liu (2008) and ARTS  databases are available within the RTTOV-165 SCATT hydrotable generator.

Mie spheres
The Mie theory is solved using an iterative method that computes a set number of terms from the infinite Mie series, using recursion relations to evaluate the required polynomials (see e.g. Ulaby et al., 1981). These calculations depend only on the size parameter x g = πD g /λ (where D g is the diameter of the sphere and λ is the wavelength) and the complex refractive index 170 of the material composing the sphere, n = √ , where is the complex permittivity. For spheres composed of liquid water the permittivity models of Liebe (1989), Kneifel et al. (2014) and Rosenkranz (2015) are available. These models were evaluated in a weather forecasting context by Lonitz and Geer (2019 scattering brightness temperature depressions at lower microwave frequencies, but it generates insufficient scattering at higher frequencies, failing to reproduce the observed behaviour (Geer and Baordo, 2014). This is primarily due to excessive forward scattering from the idealised spherical particle (Sect. 4). However, the Mie capability is still used, where appropriate, to fill gaps in the scattering databases. Specifically, this supports an optional extension of the size ranges below the smallest available particles from the databases. It is also used to fill the gap below 3 GHz where the Liu (2008) database does not provide data. resulting estimates of melting particle optical properties are placed in the 273 K temperature bin of the lookup tables. Melting particles can increase microwave brightness temperatures by 2 to 8 K over radiatively cold surfaces, mainly at frequencies of 37 GHz and below (Bauer, 2001). The equivalent bright band effect would be important for simulating radar reflectivity.
However, DDA calculations from partially melted ice aggregates show that sphere-based models perform poorly . The representation of non-spherical melting particles is a matter of ongoing research. Hence, melting particles 200 and bright band effects are not represented by default in the hydrotable generator; this awaits the availability of realistic nonspherical melting particles in scattering databases.

Liu (2008) non-spherical frozen particles
Use of the Liu (2008) scattering database for non-spherical ice particles revolutionised the quality of microwave scattering simulations made by RTTOV-SCATT (Geer and Baordo, 2014). Table 3 lists the options, of which the sector snowflake was the 205 previous default choice for snow. The other options are a dendrite snowflake and a variety of hex-plates, columns and rosettes.
These habits are geometric models of ice particles which are rescaled and then discretised into a 3D grid of polarizable points for input to DDA scattering calculations. To create the database, single-particle scattering properties have been averaged over a large number of random orientations, and these averages have been tabulated at a range of particle sizes (Liu, 2008, their Table 3, with D max being the largest size available from the database and D min the smallest, but bounded at 100 µm on the assumption that the Field et al. (2007) 215 PSD will be used (see Sect. 3.1).
The geometric particle model also provides the link between the particle size (the geometric diameter or maximum dimension D g ) and its mass m, via the mass-size relation: Coefficients a and b describe the mass-size relation m = aD b g and are in SI units; see the code for full numerical precision. ARTS IDs are unique to RTTOV and do not correspond to Eriksson et al. (2018); Liu IDs do correspond to Liu (2008). *ARTS standard habits with IDs from 10 to 15 are a mixture of two habits, with the small size range covered by thick plate, long column, block column, short column, gem cloud ice and 8-column aggregate respectively. There is some ambiguity in the fitting of these coefficients to a particle model (Geer and Baordo, 2014). In this work, a and b coefficients appropriate to the Liu particle models have been taken from Kulie et al. (2010, their Table 1). This choice means that some particles are affected by a slightly unrealistic choice of ice density (Geer and Baordo, 2014), but these values have been retained with the aim of consistency with earlier results.
The optical properties of the Liu particles are illustrated later in Fig. 9 (or see also Kulie et al., 2010;Geer and Baordo, 2014). A limitation of the Liu database is the relatively low diversity among the bulk optical properties achievable using the 225 different habits. For example, the 4-, 5-and 6-bullet rosettes give similar results to the sector snowflake. Then, there is a big gap to the next-most scattering particle, the 3-bullet rosette, and another big gap to the intensely scattering hex. plates and columns which provide similar results as each other. However, these five hex. particles with b = 3 provide a uniquely strong bulk extinction that cannot be obtained from the ARTS database (further discussion in Sect. 4).

230
The ARTS scattering database  was created to support sub-mm frequencies and to provide a broader range of non-spherical ice particles, including a variety of aggregates and densely rimed particles (e.g. hail and graupel). Table 3 summarises the options available. The current default frozen particles in RTTOV-SCATT are based on ARTS particles (see Table 1).
The ARTS database provides optical properties at 34 frequencies from 1 to 886 GHz, at three temperatures (190, 230, and 235 270 K), and at least 34 sizes per habit. The ARTS standard habits are used here; these simplify the application of the database by ensuring a full coverage of size, temperature and frequency. The size issue is that in the underlying database, the smallest size of aggregate habits can exceed 200 µm, so where necessary the standard habits consist of a habit mix. In these cases, the small size range is covered by a single crystal habit having similar shape as the constituents of the aggregate (see Table. 3). For example, the "Large Plate Aggregate" habit is complemented with the "Thick plate" habit. To avoid discontinuities, there is a 240 linear transition between the two habits over a certain size range. These "mixed" standard habits are throughout named after the habit covering the main size range and contain at least 42 particle sizes. Remaining standard habits are essentially a copy of the original ones. To improve temperature coverage in the standard habits, points at 210 and 250 K are added by a second order interpolation, in order to decrease the error by subsequent linear temperature interpolation. Finally, due to limitations in DDA, there are some gaps in the database for combinations of large size and high frequencies. In the standard habits, these 245 gaps are filled by copying data from lower frequencies that should be a better approximation than setting the values to zero.
When producing the RTTOV-SCATT hydrotables, the standard habit data are interpolated to the required temperature, frequency and particle size using 3D linear interpolation, with linear extrapolation also permitted up to a limit. This is used to provide values outside the available temperature range, but extrapolation in frequency or size is not used. This is because the frequency range already covers all current instruments, and the size range of the available particles from Table 3 is also used as 250 the integration range D min to D max . The bulk optical properties of the ARTS particles will be explored in the rest of this work: see in particular Figs. 2, 9, 10 and 12. In most places in this work, the shortest unambiguous name (such as "ARTS plate") is used. Computations have been done for a water content l = 1×10 −4 kg m −3 and temperatures T = 253 K for frozen particles and T = 283 K for liquid particles. Frequency steps are logarithmically distributed with 20 points per decade. Figure 1 shows the spectral variation of bulk optical properties across the microwave and sub-mm frequencies. These properties 255 have been generated from the default 5-hydrometeor configuration (Table 1) for a water content l = 1×10 −4 kg m −3 . A Mie sphere snow particle is also included to support discussions in Sec. 4. Rain and cloud water have broadly similar extinction per mass of particles (panel b), but because the rain particles are large enough to be in the Mie regime across most of the frequency range, they have much increased single scattering albedo (up to 0.5, panel a), asymmetry (up to 0.8, panel b) and radar reflectivity (up to 27 dBZ, panel d). Cloud water starts to depart from the Rayleigh regime above around 500 GHz, with 260 non-zero values of the single scattering albedo and asymmetry. Moving to snow, graupel and cloud ice, these have generally much lower extinction than the liquid particles below 100 GHz, but this reverses above around 300 GHz. Snow and graupel provide substantial scattering above around 20 GHz (SSA > 0.3), becoming almost purely scattering particles above 150 GHz (SSA 0.95). Cloud ice has an order of magnitude less extinction than snow and graupel at 100 GHz, and much less scattering (lower SSA) across the whole frequency range. Compared to snow and graupel, this arises mainly from the choice of PSD, 265 which gives much smaller numbers of large particles. Figure 1 hence shows the "spectral signatures" of hydrometeors and illustrates the utility of making measurements across the whole of the microwave and sub-mm in order to characterise the physical details of cloud and precipitation particles.

Bulk optical properties
The only difference between snow and graupel in the default configuration (Table 1, Fig. 1) is the use of, respectively, an ARTS large plate aggregate and an ARTS column. The primary resulting difference is the asymmetry, with graupel giving 270 less forward scattering between 50 GHz and 500 GHz. This allows the graupel to generate deeper brightness temperature depressions (see Fig. 8 later). This greater "scattering" ability led to its selection as a reasonable representation of convective snow (Geer, 2021b).
In Fig. 1 the frozen particles have small oscillations with frequency, particularly obvious in the radar reflectivity at lower frequencies. This is a result of interpolating away from the original temperature, size and frequency steps in the ARTS database.

275
This should not be an issue near the frequencies of typical satellite channels, since the ARTS frequencies have been chosen with this in mind . Figure 2 illustrates the full range of frozen particle representations available from the ARTS database, along with the Mie sphere. The ARTS particles fall into two classes. The first is less dense particles with branched shapes including rosettes, snowflakes and most of the aggregates. These are shown with solid lines and typically generate smaller SSA, extinction, 280 asymmetry and radar reflectivity. The second class is denser and more compact particles including pristine crystals, densely rimed particles (graupel and hail) and the 8-column aggregate. These are shown with dashed lines and typically generate higher values of all the optical properties. Further discussion, and comparison to the Liu (2008) particles, is made in terms of brightness temperature in Sec. 4.

Importance of mass-size relation 285
The mass-size relation (Eq. 5, specified by the a and b coefficients from Table 3) plays an important role in controlling the bulk optical properties derived from non-spherical frozen particles. In the microwave and sub-mm, the primary control over the single-particle optical properties is the particle's mass (e.g. Eriksson et al., 2015). Hence the mass-size relation already describes a lot about how particle size (as specified by the assumed PSD) maps onto optical properties. Further, as will be shown in Sec. 3.1, the mass-size relation also affects the shape of the PSD itself.

290
To summarise the available mass-size options, Fig. 3 shows those of the ARTS particles, illustrated using the effective particle density: This is the mass of the particle m(D g ) divided by the volume of a sphere of diameter D g . For spherical particles, the effective density and true density are equal. The true particle mass is reported in the particle databases and is the mass of ice used in 295 the DDA calculation; these can vary slightly from the fitted mass-size relation (see e.g. Eriksson et al., 2018, their Fig. 11).   Within the hydrotable generator, it is the mass-size relation that is used to estimate the particle mass where required (primarily, in the derivation of the PSD, Sec. 3.1) but this is an approximation. In Fig. 3 the particles with b close to three (the hail, graupel and 8-column aggregate) have almost constant effective density as a function of particle size. Most of the other particles have   densities may be generated, since that is how they are used in Eq. 9 (later). The legend is ordered by the effective density at 1×10 −3 m.  Eriksson et al. (2018) is similar but is based on the reported particle masses, rather than the mass-size relation.
b closer to 2 and hence their effective density decreases strongly with size. Some mass-size relations generate non-physical 300 super-dense particles when taken out of their validity range (as mentioned in Sec. 2.1.1, the assumed density of pure ice in RTTOV-SCATT is 917 kg m −2 ). This is relevant because when the PSD is fitted analytically to the water content (Eq. 9; Sec. 3.1) any superdense region will be included. But as explored in Sects. 3.2.4 and 3.2.2 this is of little practical relevance, and even in the sub-mm the bulk extinction is insensitive to ice particles smaller than 100 µm.
In an ideal world, users would impose their own constraints on the mass-size relation. For example, a certain mass-size 305 relation may be assumed within the physics of the forecast model which supplies the cloud and precipitation profiles, and it may be intended to achieve microphysical consistency throughout the modelling chain. Further, there are observational constraints, and for example the Brown and Francis (1995) mass-size relation gives a good description of midlatitude stratiform ice cloud (Hogan et al., 2012); this is shown on Fig. 3. However, there is currently no way of decoupling the mass-size relation from the DDA particle choice in the hydrotable generator; this could only be achieved by choosing an appropriate (and probably 310 different) database particle for each size bin -in other words, a particle ensemble approach, which is not yet supported. The Mie sphere does allow a free choice of mass-size relation (Sec. 2.1.1) but obviously brings many other drawbacks, so this is not advised. However, Fig. 3 shows that the available DDA particles span a wide range of mass-size possibilities; the dendrite particle in the Liu database (Table 3) would almost exactly match the Brown and Francis (1995) mass-size relation, for example. However, microwave radiances have their strongest sensitivity to convective snow particles, both in the tropics 315 and midlatitudes, so the appropriate mass-size relation remains poorly known. Hence the dominant approach is to explore all potential DDA particle choices and to use the one that provides the best fit between model and observations (e.g. Geer, 2021b, and references therein). Interestingly, the best choices in that work, reflected in the default RTTOV-SCATT configuration (Table 1) seem to be particles with b around 2.0 -2.4; particles with b closer to 3 seem to work poorly as a description of convective snow.

320
The sensitivity of optical properties to the mass-size relation is further illustrated by the sector snowflake, which is present in both ARTS and Liu (2008) databases and has almost identical optical properties as a function of particle size. However, as shown in Table 3, the a and b coefficients used with the Liu and ARTS databases are different, due to different, but equally valid methodological choices in fitting those coefficients to the particle masses within the databases (see Geer and Baordo, 2014, appendix B for further explanation). These small differences still have a significant effect on the bulk optical properties. The

325
ARTS sector snowflake provides less scattering than the Liu equivalent, and simulations of very thick clouds using the ARTS sector snowflake and the Field et al. (2007) PSD can be up to 20 K warmer around 300 GHz (shown in Fig. 9 later). Using an identical a and b it is possible to eliminate this difference. However, it is preferred to retain the values previously used with the Liu database to ensure full back-reproducibility, but for the ARTS database to use the coefficients that are supplied with that database.

330
Further discussion on the importance of the mass-size relation is found in Sects. 3.1 and 3.2.4. In this work, it is important to realise that when the bulk optical properties of a particle habit are discussed, this is the net result of both the physical characteristics of the individual particles and the effect of the corresponding mass-size relation on the PSD.
This follows the universal framework of Petty and Huang (2011). The version of the MGD used here is based in geometric diameter or maximum dimension, D g , consistent with the majority of available PSD formulations. n g (D g ) is the number The moments of the MGD, M k , can be derived analytically (Petty and Huang, 2011):

345
where the Gamma function Γ(z) = ∞ 0 x z−1 exp(x) dx arises naturally from the integration of the MGD. This is computed in the table generator by means of a built-in Fortran function.
The PSD is fitted to the hydrometeor water content l. Where a power law mass-size relation is known, and a and b are its coefficients (Eq: 5), l is proportional to the bth moment of the PSD: 350 Typically all but one parameter of the MGD is prescribed and the remaining "free parameter" is adjusted to fit the hydrometeor water content. The table generator allows either N 0 or Λ to be the free parameter since the other two are less mathematically convenient. These are hence computed from Eqs. 8 and 9 as follows: where p = (µ + b + 1)/γ. There are a couple of issues with the analytical approach: first, any numerical integration of the PSD, such as to obtain the bulk optical properties, is necessarily done over a limited size range D min to D max (see Eq. 1); second, some particles with b < 3 in the mass-size relation can generate non-physical super-dense small particles (Sect. 2.2.1 and Fig. 3). In the hydrotable generator these are partially dealt with via "renormalisation", an empirical rescaling of the PSD 360 described in Sec. 3.2.2.
The PSDs available in the table generator are summarised in Table 4. The Marshall and Palmer (1948, MP48) PSD is used for rain in the default configuration and is also a possibility for snow. This PSD has µ = 0 and γ = 1, producing what is classed as an exponential distribution; in the table generator, fixed values of N 0 are specified for liquid or frozen hydrometeors and Λ is the free parameter (see Table 4). It is optionally possible to add a temperature dependent N 0 (Panegrossi et al., 1998, appendix) 365 which represents the collection of smaller droplets by larger drops during sedimentation.
The "gamma" distribution, where γ = 1, is often used for cloud water or cloud ice. Typically µ and Λ are prescribed and N 0 becomes the free parameter. The default configuration for cloud water follows this approach with fixed parameters that ensure cloud water particles are in the Rayleigh regime at microwave frequencies (see Table 4). The pre-v13 equivalent is retained for back-comparison purposes; this used an alternative power law fit to Eq. 10, but the resulting difference in bulk optical 370 properties is minimal. Prior to v13.0, the gamma distribution was also used for cloud ice, also with an alternative formulation for Eq. 10. An MGD implementation is examined later in this section, labelled "MGD A". However, the default PSD for cloud ice at v13.0 ("MGD B") was identified by parameter estimation (Geer, 2021b) and is similar to the Heymsfield et al. (2013) PSD where its distribution becomes close to exponential (µ = 0).   Table 4) made the particles very small and the simulated ice cloud was almost invisible at frequencies of 183 GHz and below. Geer (2021b) explored other options, hoping to make cloud ice more visible, as seen in observations (e.g. Doherty et al., 2007;Hong et al., 2005). A number of aircraft-based PSDs were tested (H13 S, F07 M, MH97) but all produced too much scattering, even when ), an ice water content of 1×10 4 kg m −3 , and temperatures of 223 K (dashed) and 263 K (solid). PSDs have been split across two panels for clarity but for reference F07 T is present in both. At 223 K, the H13 S and H13 C distributions are nearly identical below n g (Dg) = 10 7 m −4 . See Table 4 for PSD name abbreviations. PSDs have been renormalised (Sec. 3.2.2) and all have been integrated using the new approach (Sec. 3.2). the particle type was chosen to generate as little scattering as possible (e.g. the ARTS large column aggregate). Figure 4 shows 390 that these PSDs can generate significant numbers of larger particles, particularly at warmer temperatures, which must be the cause of the excess scattering. To fill the gap between the previous gamma configuration (e.g. MGD A, too few large particles) and the aircraft-based PSDs (too many large particles), new PSDs were created, inspired by the low-temperature part of H13.
The configuration labelled here as MGD B (see Table 4) was ultimately chosen as the cloud ice default. Once sub-mm data from ICI is available, it will be seen whether this is indeed a physically reasonable choice; however it was shown to do a 395 reasonable job in representing observations at 183 GHz.
An issue with many ice PSDs, and particularly evident with F07 and MH97 on Fig. 4, is the presence of a "small mode" of ice particles. The aircraft measurements on which these PSDs were based were subject to probe shattering (Korolev et al., 2011) and optical effects (O'Shea et al., 2020) that, it is now thought, create a spuriously large number of small particles.
Hence the small-size mode of these distributions might be non-physical. The H13 PSD is intended to be free from at least the 400 probe shattering effect (Heymsfield et al., 2013).

Field-type PSDs
The Field et al. (2005Field et al. ( , 2007 PSDs are based on a "universal" rescaled PSD Φ 23 (x 23 ), which is a function of a non-dimensional Here, M 2 and M 3 are the second and third moments of the PSD (Eq. 8) but any pair of moments could have been used. The universal PSD is parametrised as the sum of exponential and gamma PSDs in x 23 which 405 gives the resulting PSD a characteristic population bulge in the smaller sizes (Fig. 4).
The size-based PSD is recovered by To evaluate the PSD hence requires knowledge of M 2 and M 3 , or equivalently any other pair of moments; these are obtained by an empirical relation that converts one moment to any other (e.g. Eq. 3 in Field et al., 2007). The water content l provides 410 M b through Eq. 9; this is first converted to M 2 and then M 2 is used to obtain M 3 . The universal PSD is not itself temperature dependent, but Field et al. (2007) provides two parameterisations, one for tropical and one for midlatitude conditions. The temperature dependence arises through the empirical relation between moments, so that the F05 and F07 PSDs generate smaller particles at lower temperatures (Fig. 4).
There are some issues to consider with the Field PSDs, as well as just the small-particle mode noted before. First, the aircraft 415 observations on which they were based did not measure particles with D g smaller than 1×10 −4 m (100 µm). The universal PSD can be used to extrapolate to smaller sizes; the hydrotable generator allows this for the F05 PSD. Field et al. (2007) recommended more strongly not to extrapolate, so the table generator terminates the F07 PSD at D min =1×10 −4 m (Fig. 4).
When the above procedure is followed to define a size-based PSD from the ice water content, it is assumed that it is valid with an integration over sizes from 100 µm to infinity. The numerical integration of the resulting n g (D g ) and particle mass m(D g ) 420 (following Eq. 9) should recover the original ice water content, but instead the results can be very different; this is covered in Sec. 3.2.2. The F07 T PSD has proved useful in fitting real observations and it is vital to the default configuration of the table generator; however users need to be aware of these complex issues.

McFarquhar and Heymsfield (1997) PSD
The McFarquhar and Heymsfield (1997) PSD is based on mass-equivalent diameters D e where D e = 6m πρ ice Here m is the mass of the particle and ρ ice is the density of solid ice (note the table generator uses ρ ice = 917 kg m −3 compared to ρ ice = 910 kg m −3 in the original work). Similar to the Field PSDs, it has two modes (Fig. 4b): the first represents particles smaller than D g =1×10 −4 m (100 µm) using a gamma distribution in D e . Larger particles are represented by a lognormal distribution, also in D e ; this cannot be represented in the universal MGD framework of Petty and Huang (2011). There is no 430 hard cutoff between the distributions, rather they are summed for all D e from 0 to ∞, observing that the two distributions do not have a big overlap. To adapt the PSDs to the specified ice water content l, the water content is first split into two parts representing the small and large particles. This is done based on an empirical relation (McFarquhar and Heymsfield, 1997, Eq. 5). The small-and large-particle PSDs are then dependent on the partial masses (their Eqs. 3 and 4). The parameters of both PSDs are also temperature dependent (their Eqs. 7-12), producing behaviour broadly similar to Field et al. (2007, see Fig. 4b).

435
A PSD based on geometric diameter is recovered by the conversion where dDe dDg has been evaluated numerically. The MH97 PSD is less sensitive to the choice of mass-size relation and hence less sensitive to variations in the particle habit (see Fig. 12). This is not, it is thought, because it is based in mass-equivalent diameter D e , as hypothesised by Eriksson et al. 440 (2015), but because it puts so much of the mass in the small particle mode (Ekelund et al., 2020b). As with the Field PSDs, this small-particle mode may be physically incorrect and may have been generated by probe-shattering or optical effects (Korolev et al., 2011;O'Shea et al., 2020).

Integration methods
The core of the hydrotable generator is the numerical integration over the PSD to produce the bulk optical properties, as 445 described earlier (Eqs. 1 to 4). This section first describes the more technical aspects of the integration: numerical integration methods, renormalisation, and diagnostics (Sects. 3.2.1, 3.2.2, and 3.2.3 respectively). Then Sect. 3.2.4 explores the scientific importance of the integration, illustrating the size range of particles that contribute to the bulk optical properties, and helping to explain the impact of different microphysical choices.

Numerical integration 450
In previous versions, numerical integration was done at fixed steps in particle size D g , using a rectangle rule integration, centred on the integration point. The current version also offers an improved integration using the trapezium rule, and with log-spaced integration points in D g to better resolve the small size ranges. The number of integration points is fixed at 100 and is the same in both methods. The integrations in Eqs. 1 to 4 use a PSD that has been renormalised to conserve integrated mass, n g (D g ); this is described in Sect. 3.2.2. For reasons to be explained, a mix of the old and new integration techniques is used in the default 455 configuration ( Table 1).
The integration of optical properties is done over the truncated range D min to D max . For Mie spheres, the integration range is given in Table 2. For particles from the Liu or ARTS databases, the integration size range is taken from Table 3, with two exceptions. First is that the size range can optionally be extended down to the relevant D min from Table 2. If selected, the relevant optical properties are computed using Mie theory, assuming this is valid for particles smaller than the minimum 460 particle size (see Sect. 2.1). If that option is not selected, then a minimum size D min = 1×10 −4 m is applied when the F07 PSD is used with the ARTS shapes, to avoid extrapolating the PSD. Note that as described in Sect. 2.1, for the implementation of the Liu database in the hydrotable generator, the D min = 1×10 −4 m constraint was imposed unilaterally in Table 3, meaning that it affects the Liu particles no matter which PSD is chosen (this behaviour is undesirable but is preserved for back-compatibility).
The numerical integration methods are illustrated in Fig. 5 using the computation of the implied water content (Eq. 15, next 465 section). The old method is represented by the stepped red line; its grid was too coarse to resolve sharp PSD features in the small size ranges. The trapezium rule used in the new method is represented by straight lines between log-spaced integration points (indicated by the black crosses in the example with the D min = 1×10 −4 m (100 µm) cutoff). The new method is a more exact representation of the integration range D min to D max , whereas in the old method the centred bins extended all the way down to D g = 0, indeed fractionally beyond in some cases. This means that with the old method, even if the nominal D min 470 was significantly above zero, the integration was still roughly representing the full size range of particles from zero to infinity.
Hence if the intention were to exclude the smallest particles from the PSD, such as when the 100 µm cutoff is used with the F07 PSD, then the old integration scheme does not fully achieve this aim. The importance of this is examined in the next section.

Renormalisation
An important test of the numerical integration is whether the specified water content l, used to specify the PSD, can be recovered 475 in the implied water content when the particle mass and PSD are numerically integrated across the chosen integration range D min to D max : The reconstructed water content may be different due to deficiencies in the numerical integration, if the chosen particle size range omits part of the PSD, or if there are inaccuracies in the conversion from water content to PSD parameters.

480
A renormalisation factor r can be computed: In order not to lose or gain mass, the PSD is renormalised as follows: n g (D g ) = r n g (D g ).
All the PSDs used in this work are renormalised, with the exception of those shown in Fig. 5. In most cases the renormalisation 485 is minor, with |r − 1| less than 0.03, often much smaller. However there are some exceptions. As shown in the figure, the F07 T PSD requires relatively large amounts of renormalisation, with |r − 1| of 0.13 and 0.04 in this example (using the new integration with the 100 µm cutoff, for the cold and warm temperatures respectively). The other main execptions are the MGD A PSD, which has |r − 1| = 0.4, and the F05 PSD, which has |r − 1| = 9.1 at the low temperature setting where the problem is worst (for the same example, but not shown on Fig. 5). Since these are the PSDs with the largest number concentrations in the 490 smallest sizes, this illustrates how the failure of Eq. 15 is typically due to a large amount of mass, and/or sharply peaked PSD features in the very smallest size ranges. This can make numerical integration difficult.
In the case of the F05 and F07 PSD, the extrapolation of the PSD below D min = 1×10 −4 m exposes the question of whether to use this portion of the Field-type PSDs. Comparing panels Figure 5a and b shows that any issues are most relevant for the smallest water contents. A question is whether to truncate the F07 PSD at D min = 1×10 −4 m as recommended (Field However the problem is worse for smaller water contents. Overall, with the F07 PSD, the least renormalisation is generated with the old integration basis (see Sec. 3.2) and with the 100 µm cutoff; hence even at v13.0 these are the approaches used in the default configuration of hydrometeors based on the F07 PSD. This has the advantage of retaining comparability of the results with earlier work (e.g. Geer and Baordo, 2014). But it is important to realise that the results coming from the F07 PSD are dependent on these choices.

505
Renormalisation is always active in the table generator, but to alert the user to any significant issues, it will throw an error if the order of magnitude of renormalisation |log 10 (r)| exceeds a pre-set threshold. For hydrometeors using the MP48 and MGD PSDs, the thresholds are 0.05 or less, showing they are not much affected. For the hydrometeors using the F07 PSDs, the threshold has to be 0.5. However, the largest renormalisations are for the smallest water contents, meaning the issue does not in most cases have a significant influence on the final simulated brightness temperatures.

Diagnostic mode
As illustrated in this subsection, there are many complexities to the apparently simple task of numerical integration of bulk optical properties or mass, particularly since many PSDs put significant mass in the smallest size ranges, where the particles are unimportant in the microwave and sub-mm radiative transfer. Since it has not been possible to demonstrate or test every combination of options provided by the tool, users may wish to check the quality of the integrations for themselves. If the 515 amount of renormalisation required is large, this is an early warning of problems, but even better is to make use of the new diagnostic mode, which writes out an additional diagnostic text file during the generation of the lookup tables. For a chosen particle ID (from Table 2), temperature, frequency, and water content, the diagnostic mode outputs the values of key parameters at each integration point: D g , D m (D g ), m(D g ), N g (D g ), β e (D g ), β s (D g ), β b (D g ), g single (D g ) and the "extagrand", the numerical representation of β e N g (D g )dD g , which is summed to create the final bulk integrated extinction. The resulting bulk 520 values are also provided, along with the renormalisation factor r to be able to recreate N g (D g ). The new diagnostic mode was heavily used in the development of v13.0 and in the preparation of figures for this paper. Figure 6 illustrates the integration of extinction (Eq. 1) for frozen particles at 190.31 GHz using the F07 T PSD. The integration combines the single-particle extinction (β e (D g ), panel a) and the PSD (N g (D g ), panel b), so that each integration element gives 525 a contribution of β e N g (D g )dD g to the bulk extinction (panels c and d), referred to here as the "extagrand". The elements are logarithmic in this example (the "new" integration has been used) and the size axis of the plot is logarithmic; hence the bulk extinction is proportional to the area under the curves in panels c and d. All panels have been normalised: the extinction and the size distribution have been respectively divided and multiplied by the single-particle mass (based on Eq. 5), presenting Eq. 1  Figure 6. Integration of single-particle extinction over the PSD: (a) single-particle extinction per unit particle mass βe(Dg)/m(Dg); (b) mass distribution m(Dg)N g (Dg); (c) contribution to total extinction ("extagrand") normalised by water content βe(Dg)N g (Dg)dDg/l for ARTS ICON hail; (d) as (c) but for ARTS plate aggregate. The F07 T PSD has been used with the "new" integration, at a temperature of 253 K and a frequency of 190.31 GHz. The extended integration range for cloud ice (down to 5×10 −6 m) has been used for illustrative purposes. Dotted lines on (a) correspond to, first, the largest sizes for which the plate aggregate is super-dense; second, the largest size for which Rayleigh scattering would be an appropriate model for the ICON hail particle. as follows:

Converting mass to extinction
This has two aims; first to normalise quantities that would otherwise vary over more than 10 orders of magnitude; second and most importantly, to focus on the key process in the computation of bulk optical properties, which is to convert hydrometeor mass to bulk extinction. Further, to more easily compare the results with different water contents in panels c and d, the extagrand has been normalised by the respective hydrometeor water content. Panel b shows that the effect on the F07 PSD of increasing 535 the hydrometeor water content is not just to increase the overall mass of particles, but also to significantly increase the maximum particle sizes included.
For the ICON hail particle and a water content of l = 1×10 −2 kg m −2 , almost all the extinction is generated by particles with D g between 3×10 −4 m and 5×10 −3 m, in other words particles with a maximum dimension of around 1 mm. This corresponds both to the peak in the mass-weighted PSD and the peak in the per-mass extinction. This peak in per-mass extinction could be 540 called the "resonance" zone: particles with sizes that are a little larger than the wavelength give particularly large extinction even without normalisation by mass (see e.g. Petty, 2006, their Fig. 12.4). The ARTS plate aggregate is a less dense particle and the resonance zone is found at larger particle sizes. Hence the size range contributing to the bulk extinction is between 3×10 −4 m and 1×10 −2 m. The extension of the PSD to slightly larger particles (because the plate aggregate model implies different parameters in the mass-size relation) also contributes to this. But the range of particle sizes which contribute to the 545 bulk extinction is much smaller than the range of the mass-weighted PSD. In other words, there is a significant amount of particle mass with sizes smaller than 1 mm that is mostly or completely "invisible".
Within the Rayleigh scattering regime it is broadly the mass of ice, and not the particle shape, that controls the single-particle scattering properties 2 . For example, Fig. 12 of Eriksson et al. (2018) shows optical properties of non-spherical ice particles converging for x e < 0.5, where the size parameter x e = πD e /λ is based on the effective (mass-equivalent) particle size, not 550 the maximum dimension. This is a more relaxed definition of the Rayleigh regime than often suggested (x e < 0.1 is typical) but using this, the ARTS ICON hail particle departs the Rayleigh regime above D g = 2.6×10 −4 m at 190 GHz. However, even in the small particle limit, the extinction per mass shown on Fig. 6a is different between the ICON hail and the plate aggregate. This is a potentially confusing aspect of using the maximum dimension D g as the x-coordinate. Because the ICON hail particle is significantly denser, it thus has higher mass for the same D g and hence even after normalisation by particle 555 mass, it still has a disproportionate effect on the radiation field for the same particle size D g . Hence even in the Rayleigh regime, particle morphology still needs to be taken into account when mapping from particle size D g to particle mass; this is not a completely obvious point given that Rayleigh and Mie sphere optical properties are typically described in terms of sphere diameter, rather than mass. This also further illustrates the importance of the mass-size relation of the particle model (Eq. 5; Table 3) in determining the bulk optical properties. Interestingly, in the mass-weighted viewpoint of Fig. 6b, changing from the Water contents are coloured as indicated in the legend. The dashed black line indicates the limit of the Rayleigh regime for this particle; above this, non-spherical optical properties must be computed using DDA or equivalent method. mass-size relation of hail (exponent b = 2.99) to that of the plate aggregate (b = 2.26) has only a secondary effect on the PSD shape.
A minor issue is that some particles with an exponent in the mass-size relation b < 3 (Eq. 5; Table 3; Fig 3) can be "superdense" at small sizes, in other words that the mass-size relation implies a particle density that is higher than that of solid ice.
This affects the plate aggregate (but not the ICON hail) below around D g = 2.7 × 10 −5 m. However, in the computation of the 565 single-particle optical properties, whether Mie theory or DDA, particle densities are in practice not allowed to exceed those of solid ice. This could in theory result in an incorrect calculation of bulk extinction, but Figure 6 shows that even for the smallest water contents, any issue with representing superdense particles is irrelevant from the point of view of the optical properties, since particles as small as these are invisible. However, the treatment of small particles does affect the distribution of mass within the PSD, and hence can affect the bulk optical properties through renormalisation as explored in Sec. 3.2.2. 570 Figure 7 explores the frequency dependence of the range of particle sizes that are optically relevant; this is based on the ARTS ICON hail in order to show a particle which departs the Rayleigh regime at relatively low frequencies, here 10 GHz.
The transition to non-Rayleigh scattering is associated with nearly an order of magnitude increase in the minimum particle size that contributes to the bulk optical properties (indicated by the 2.5 percentile of the integration here). However (and as also suggested by Fig. 6) the maximum particle size is more controlled by the PSD, and hence also the water content, and is 575 less variable with frequency. This means that the size range contributing to the bulk optical properties is particularly squashed in the "resonant" region, which occurs just above the Rayleigh regime. One broad conclusion is that sophisticated models for non-spherical particle scattering are always required to correctly simulate ice hydrometeor optical properties at microwave and sub-mm frequencies, even for very small water contents, and even for PSDs that do not generate such large particles as the F07 T PSD (compare Fig. 4). The regions of the frequency and particle size spectrum where an approximate solution (such 580 as Rayleigh or Mie) would be valid are those where the particles would be mostly invisible anyway. The ICON hail is the most dense available ARTS particle (Fig. 3) and would generate the most scattering from particles that are small in D g . Hence this confirms that sub-100 µm (D g < 1 × 10 −4 m) ice particles are irrelevant to the radiative transfer up to 886 GHz. For an alternative view, excluding the constraint from the PSD, but finding similar results, see Pfreundschuh et al. (2020, their Fig. 5).

585
The Liu and ARTS particles available in the hydrotable generator (and obviously the Mie sphere) represent only randomly oriented particles. However, ice hydrometeors are often preferentially oriented, as revealed by polarisation signatures in the high-frequency channels of microwave imagers (e.g. Defer et al., 2014;Gong and Wu, 2017). The ARTS database has recently been extended with a more advanced representation of particle orientation, giving particles a preferred canting angle, but retaining random orientation in azimuth . However, this generates optical properties that are fully polarised, 590 i.e. scattering can transfer energy from one polarisation to another. To model such fully-polarised optical properties would require the whole Stokes vector to be represented in the radiative transfer, but this is not available in a scalar fast model like RTTOV-SCATT. However, it is still possible to represent much of the effect of preferential orientation on microwave imager brightness temperatures, using an approximate method. This is done by scaling the bulk extinction β e , as generated from fully randomly oriented particles, according to a polarisation ratio ρ (Barlakas et al., 2020): Here, β e is increased by the proportion α to provide the scattering coefficient for use in horizontally (H-) polarised channels β e,H . Similarly, it is reduced by α to provide the scattering coefficient for vertically (V-) polarised channels β e,V . This description is more than just a tuning factor: αβ e describes the bulk scattering cross-section for linear polarisation in fully-polarised radiative transfer, which represents the differences in the extinction between V-and H-channels (Barlakas et al., 2020). and will require further scientific development. Further, an approach for radar backscattering needs to be developed. The "full" channel representation provides a framework for future support of single-particle optical property databases based on oriented particles (e.g. Brath et al., 2020). Although the bulk optical properties (e.g. Fig. 2) are already informative, their effect on radiation fields is both situationdependent and a complex function of the optical properties. For example the effect of scattering on cloud-top brightness temperatures depends not just on the scattering coefficient but also the phase function, as summarised here by the asymmetry parameter. Further, even a relatively small amount of thermal emission within a cloud can substantially increase its brightness temperature compared to a purely scattering case. Hence there is a need for a standardised and simplified way to compare the 620 cloud-top brightness temperatures arising from different choices in computing the bulk optical properties. To do this we use the two-stream solution for the radiance at the top of a uniform cloud layer taking into account both scattering and thermal emission/absorption within the cloud (Appendix B): Here, I ↑ (τ = 0) is the upwelling radiance at the top of the cloud layer, which is assumed isotropic within each hemisphere in the two-stream approximation. The vertical coordinate is optical depth τ which is 0 at the top of the cloud and τ * at the 630 bottom, so τ * is the optical thickness of the cloud. ∆z is the geometric thickness of the cloud. The downwelling radiation at the top of the cloud is 0, and there is a steady source of upwelling radiation at the bottom of the cloud, I ↑ (τ * ) = I 0 . The upwelling radiance at the top of the cloud is given by Eq. 20 and is made up of below-cloud radiation that has been scattered or directly transmitted (the I 0 term) plus thermal emission from within the cloud, either scattered or directly transmitted to the top of the cloud (the B 0 term, where B 0 is the Planck function at the temperature of the cloud, which is uniform throughout).

635
The additional terms Φ, Υ and r ∞ are dependent only on the cloud's geometric thickness and the basic optical properties: the extinction β e , the SSA ω 0 and the asymmetry parameter g from the lookup tables. The terms Φ and Υ do not have a particular geophysical interpretation, but r ∞ is the cloud-top albedo, of most relevance to solar radiation. The emitted radiation at the top of the cloud is hence a function of the three optical properties, plus the below-cloud upwelling radiation I 0 , thermal emission inside the cloud B 0 (and thus the cloud temperature) and the geometric thickness of the cloud ∆z.  Table 1) present in a 2 km thick layer with a water content of l = 1×10 −3 kg m −3 . The cloud temperature is 253 K if frozen (snow, graupel or cloud ice) or 283 K if melted (rain, cloud water). Upwelling brightness temperature below the cloud is 280 K (solid lines) or 100 K (dashed lines).
The cloud just described is an approximate but compact description of typical situations in microwave and sub-mm radiative transfer. Gas absorption and emission have been neglected and the bottom boundary of the cloud is assumed black, so that any radiation leaving the cloud downwards can be forgotten -for example radiation reflected from the surface is ignored. This would still be a good representation of a cloud in the upper troposphere in any part of the spectrum where water vapour or oxygen absorption blocks visibility of the surface, yet is not a significant source of emission at the level of the cloud itself. It 645 is also a good representation of radiative transfer over land surfaces, where the surface is mostly black. It would be trivial to add gas absorption within the cloud and the surface-reflected term could be included but with additional complexity. But these would be a distraction from the simple standardised comparison of hydrometeor optical properties that is intended.  Table 1. In each case the cloud is 2 km thick and the water content l = 1e −3 kg m −3 . This is quite a heavy cloud of 2 kg m −2 , but this is helpful in more clearly differentiating the types of hydrometeor. The cloud temperature is 253 K if frozen (snow, graupel or cloud ice) or 283 K if melted (rain, cloud water). A below-cloud upwelling brightness temperature of 280 K could represent a window channel over land or a lower-peaking water vapour channel (solid lines). In this case, scattering from rain generates brightness temperature depressions peaking at around 40 K, and starting above 10 GHz. Cloud water is strongly 655 absorbing, but the situation has minimal thermal contrast, so it provides only a tiny boost to brightness temperatures. Scattering from the frozen hydrometeors is much more effective, generating depressions up to 200 K above 50 GHz for snow and graupel, and above 100 GHz for cloud ice.

Overview of hydrometeor choices
In Fig. 8, a below-cloud upwelling brightness temperature of 100 K (dashed lines) could represent a window channel over ocean (at lower frequencies, in horizontally polarised channel, ignoring the surface reflection) or a cirrus cloud above a very 660 strongly scattering cloud placed lower in the atmosphere. Here, thermal emission from the rain and cloud water is the main effect above around 5 GHz. Above around 50 GHz, the frozen hydrometeors become visible, with snow and graupel generating up to 30 K brightness temperature depression even below the 100 K of upwelling radiation. The standard "graupel" configuration is a little more scattering than the "snow" particle, as intended (Geer, 2021b). The effect of cloud ice in this scenario is to warm the brightness temperatures, since cloud ice has relatively low SSA. This warming effect of, for example, cirrus 665 over a strongly scattering lower cloud has surprised a number of investigators (e.g. Xie et al., 2020;Barlakas et al., 2020).
Where the dashed and solid lines join is where the clouds become optically thick and where below-cloud radiation becomes irrelevant. An interesting feature is that around 300 GHz to 500 GHz, snow and graupel clouds produce their lowest scattering TB depressions. Above this frequency of maximum scattering, these clouds start emitting more radiation again and brightness temperatures are higher. 670 Figure 9 shows brightness temperatures for the standardised ice or snow cloud using the F07 tropical PSD with all options from the ARTS and Liu databases. The amount of brightness temperature depression from the ARTS particles roughly follows the amount of extinction in Fig. 2 and hence the progression broadly from low density aggregates and snowflakes to dense rimed particles. The Evans snow aggregate is again the least scattering, with TB dropping only to 150 K, and again the ICON hail is the most scattering particle, dropping to 90 K. The frequency of maximum scattering varies from around 500 GHz for 675 the Evans snow aggregate to around 200 GHz for the ICON hail. Between the extremes of light aggregate and dense hail, ARTS database provides a good spectrum of potential scattering properties. As mentioned before, the Liu database has bigger gaps and some of the particle shapes are almost redundant. The least-scattering Liu particle is the dendrite snowflake, but this still generates a relatively deep TB depression (as low as 130 K at 340 GHz). By contrast the ARTS database provides the column aggregate and Evans aggregate, which produce even less scattering. This provides an important new capability to reproduce 680 cloud ice signatures. Hence the ARTS column aggregate was required in the study of Geer (2021b) and is used in the new default RTTOV-SCATT configuration (Table 1).
At the strong scattering end in Figure 9b, the Liu columns and thick hex plate are more strongly scattering than the ICON hail from the ARTS database, giving TBs as low as 70 K around 180 GHz. Hence these Liu shapes continue to provide a capability that is not available from the ARTS database. However it is interesting that the biggest variations in brightness temperature 685 depression are in frequencies below 200 GHz. The studies of Geer and Baordo (2014) and Geer (2021b), based on SSMIS observations, have already benefitted from a spectral region of great sensitivity to particle habit and PSD. In these studies, the most strongly scattering particle models, those which generate deep TB depressions even at 50 GHz, have been decisively rejected, at least as a means of representing the effect of snow on TBs representing a global model grid box on the 10 km -100 km scale.  Now is a good point to tidy up some loose ends from Geer and Baordo (2014), who rejected the Mie soft sphere as a viable model for snow particles at microwave frequencies. This was partly based on the excessive scattering below 100 GHz, in common with the dense non-spherical particles just described. However the soft sphere also provided insufficient TB depressions at 183 GHz; it is very clear in Fig. 9a, that it is an outlier compared to more-realistic non-spherical particle models. Geer and Baordo suggested that the primary problem of the soft sphere was not its overall level of extinction or scattering, but its unusu-695 ally high asymmetry parameter, and hence stronger forward scattering (see Fig. 2: at 183 GHz and l = 1×10 −3 kg m −3 , the asymmetry of the soft sphere is around 0.9, compared to around 0.3 for the ARTS sector snowflake). The hypothesis was not confirmed in the study, but the standardised cloud model helps to clarify. The black dashed line in Fig. 9a shows TBs generated using the Mie bulk optical properties but with the asymmetry parameter from the ARTS sector snowflake. This hybrid particle gives a TB signature that is similar in some parts to the ARTS 8-column aggregate, making it indistinguishable from the 700 DDA-based non-spherical particles. Hence it is the strong forward asymmetry of the Mie soft sphere that makes the difference.
This result shows the utility of the standardised cloud model for assessing bulk optical properties, and the importance of the asymmetry parameter (and more broadly the phase function) in determining the brightness temperature. There is reasonable sensitivity to particle shape across all the higher frequencies, with a 50 K difference between the most and least scattering shapes. For the thick clouds, TB depressions move to lower frequencies and get deeper still, with the ICON hail 710 giving at TB of 50 K at 100 GHz. At these mid-frequencies, the sensitivity to the particle shape of thick clouds is very strong, with around 150 K difference between the most and least-scattering particle. This frequency of maximum TB depression is not fixed, however, and increases to 350 GHz for the same particle in the thin cloud case. Conserving the integrated water content, but changing the water content by an order of magnitude in either direction (e.g. l = 1×10 −2 kg m −3 or 1×10 −4 kg m −3 , not shown) the broad layout of these figures remains similar. Three interesting things appear above the frequency of maximum 715 TB depression: First, the brightness temperatures start to become less sensitive to particle shape (and hence also particle size)

Optical behaviour of frozen hydrometeors in the microwave and sub-mm
which is particularly evident in the clustering on Fig. 10b; Second, there is an inversion in the ordering of scattering, with the dense particles (pristine crystals, hail and graupel) generally giving slightly warmer brightness temperatures than the less-dense aggregates and rosettes; Third, sensitivity to the integrated water content starts to disappear, with the thin cloud giving similar brightness temperatures at 900 GHz to the thick cloud.

720
To further explain the transition in cloud optical properties as the frequency increases, we can return to Eq. 20. The cloud brightness temperature results from two terms: radiation transmitted by the cloud, and radiation emitted by the cloud. The terms in square brackets multiplying I 0 and B 0 are hence the transmittance and the emissivity of that cloud. These terms are shown in Fig. 11 along with the total brightness temperature. This is based on the strongly-scattering ARTS ICON hail particle with the F07 T PSD, in order to show as much as possible of the transition within the available frequency range. The transmittance 725 of the 2 km thick cloud goes to zero towards 1000 GHz, in other words radiation from below can no longer pass through the cloud. However, the total brightness temperature bottoms out at around 100 K, rather than dropping to zero, due to increasing thermal emission from within the cloud. As illustrated in Fig. 1, albeit with a different particle type, the single scattering albedo never reaches 1 (which would imply complete scattering or zero thermal emission); instead it reaches a maximum somewhere around 500 GHz and starts to drop off at higher frequencies. This, along with the steady increase in overall extinction, must  contribute to the continuous increase in thermal emission as the frequency increases. The clouds do not need to be optically thick to show this effect; Fig. 11 also shows the 0.2 km cloud, which even at 884 GHz has a transmittance of 0.5, but still has sufficient thermal emission to contribute significantly to the brightness temperature. Hence the spectral region around 500 GHz seems to be the turnaround point beyond which cloud optical behaviours tend towards those more familiar from the infrared, with increasing optical thickness, and cloud emission becoming the dominant optical process.  Fig. 9a, the main features of the ARTS particles are preserved. MH97 generates a more compact spread of brightness temperatures across the ARTS particles, but it does not much change the overall order of scattering. The Evans 740 snow aggregate, for example, is still the least-scattering of the particles. This compact spread is attributed to the MH97 PSD putting a relatively high proportion of particles in the smaller size ranges (Ekelund et al., 2020b) and hence keeping more of the particles in the Rayleigh regime. However the bunching of the different particle models still occurs above 500 GHz, suggesting this is not strongly affected by the particle size. Hence this supports a fairly universal transition towards an emission-dominated, optically thick regime above 1000 GHz, where sensitivity to particle size and shape appears to become smaller.

Conclusion
This work has summarised the process of generating bulk hydrometeor optical properties based on physical assumptions about the sizes, masses, habits and orientations of cloud and precipitation particles. It documents the hydrotable generator for microwave and sub-mm scattering radiative transfer in version 13.0 of RTTOV (Radiative transfer for TOVS, Saunders et al., 2018Saunders et al., , 2020, a widely used satellite simulator, hopefully with relevance to the users of many similar tools. The work 750 has overviewed the bulk optical properties and brightness temperatures generated by Mie spheres and two databases of nonspherical ice particles (Liu, 2008;Eriksson et al., 2018). andHuang, 2011;Heymsfield et al., 2013) and the core process of numerical integration across them. This process maps from a chosen particle size distribution and particle model, through to individual particle masses and other physical properties, and 755 finally through to bulk optical properties. This process has a number of issues. For example, very small (< 100 µm) particles are invisible at these frequencies but can affect the bulk optical properties through renormalisation adjustments to the PSD, relating back to the underlying question whether current PSDs do a good job of representing small particles (Korolev et al., 2011;O'Shea et al., 2020). Furthermore, the effect of small differences in the particle mass-size relation has been highlighted in the results using the ARTS and Liu sector snowflakes, which have exactly the same single-particle optical properties, but 760 generate significantly different brightness temperatures due to slightly different mass-size relations.
To illustrate the available options, a standardised homogeneous layer cloud was proposed. This is based on a simple twostream analytical solution and converts the bulk optical properties into brightness temperatures, which are easier to interpret.
In particular this resolves the trade-offs between the absolute level of scattering (the scattering coefficient) and the asymmetry parameter (summarising the shape of the scattering phase function) in the amount of brightness temperature depression that 765 is generated. A further aspect of the standardised homogeneous layer cloud is that it illustrates the balance between radiation coming into the cloud from below, that may be scattered or absorbed, versus thermal emission from the hydrometeors themselves. This is particularly important for frozen particles above 200 GHz where, as we have shown, the cloud optical properties start to move towards an emission-dominated, optically thick regime that is much less sensitive to particle size and shape, a regime more familiar from the behaviour of clouds in the infrared.

770
The standardised cloud helps further investigate the soft Mie sphere as a representation of ice particles. This has unusually strong forward scattering (high asymmetry) compared to most non-spherical particle models representing the same ice mass (e.g Eriksson et al., 2015, their Fig. 3). This unusually strong asymmetry is also found in the bulk scattering properties. The standardised cloud model shows that the spherical model cannot generate sufficient brightness temperature depressions at higher frequencies. It generates a "scattering spectrum" that is very different to any more physically reasonable model. Just 775 by replacing the asymmetry of the soft Mie sphere with that from a non-spherical model, it becomes a much more plausible representation falling within the range of non-spherical models. This is further proof, if any is needed, of the problems with the Mie sphere representation of snow; Kuo et al. (2016) have made a very similar point. Only within the Rayleigh regime do spherical and non-spherical particles generate similar optical properties as a function of the particle mass. However, even at 10 GHz a significant fraction of ice particles can be outside the Rayleigh regime (Fig. 7); hence there really is no alternative to 780 taking account of the microphysical characteristics of realistic non-spherical ice particles when simulating observations in the microwave and sub-mm regions.
An underlying assumption in this work (and many others that follow the same approach) is that clouds can be represented using a single particle shape model to cover all instances of a highly heterogeneous class of particles like snow -a "one shape fits all" approach. An important aspect of this approach is that the particle shape model defines the mass-size relation; 785 this affects the bulk scattering properties both by its influence on the shape of the PSD and simply by defining the mapping between maximum dimension D g and the particle mass, which means that particles of nominally the same size can generate very different single-particle optical properties. This applies even within the Rayleigh regime. The use of mass-equivalent diameter D e could remove the latter issue (e.g. Eriksson et al., 2015) but this is of less practical relevance as long as the PSDs which are used to map from water content to particle size (and hence mass) are defined in terms of maximum dimension D g . A 790 more physically-based approach would be to consider an ensemble of particle habits (e.g. Kulie et al., 2010;Baran et al., 2011), and might impose a mass-size relation thought to be an appropriate description of certain hydrometeors, such as midlatitude stratiform ice cloud (e.g. Brown and Francis, 1995;Hogan et al., 2012). However, RTTOV has the job of modelling satellite observations with reasonable accuracy across the entire globe; microwave observations are strongly sensitive to many different cloud regimes and particularly to those such as deep convection, where even basic details of the microphysics remain poorly 795 known. Forecast models are currently incapable of representing the full range of microphysical parameters needed to constrain the hydrometeor optical properties; therefore simplification and generalisation is unavoidable, with parameter estimation being the ideal method for identifying the best microphysical assumptions (Geer, 2021b, with an appropriate Bayesian framework this could incorporate prior expert microphysical knowledge, such as the realism of particle habits, PSDs and mass-size relations).
In the event that a much richer microphysical representation is required, RTTOV-SCATT allows an unlimited number of 800 hydrometeor categories that could be used to represent almost any level of microphysical complexity. Hence the tools are already available, but the underlying issue is whether is possible to appropriately specify this complex microphysics.
There are many ways to further improve the representation of bulk scattering from hydrometeors, most of them applicable generally: -RTTOV (and many other codes) use a different mechanism and different physical assumptions for the optical properties 805 of clouds in the IR and solar regions. A major future development should be to produce optical properties from the microwave to the UV using the same approaches as much as possible; this requires non-spherical scattering databases to be expanded to cover the whole range, a process that is only starting (e.g. Yang et al., 2013;Ding et al., 2017;Baran et al., 2018).
-More work is needed to unpick the tight coupling between the particle shape model, mass-size relation and the particle 810 size distribution. Better ways of mapping bulk hydrometeor mass to particle ensembles would be very welcome.
-More work needs to be done to standardise and/or interface the physical assumptions, such as shape, PSD and mass-size relation, with assumptions made in atmospheric models. This would particularly support efforts to learn better physical models of cloud and precipitation, directly from observations (e.g. Schneider et al., 2017;Geer, 2021a).
-A global effort is ongoing to standardize and package many of the different scattering databases that are becoming 815 available (see Kneifel et al., 2018); once complete, interfaces will be added to allow the user a much broader choice of particle models. This would allow the use of non-spherical oriented raindrops, for example (e.g. Ekelund et al., 2020a).
-Current PSDs have significant limitations, such as the small-particle bulge that is now thought to be unphysical (Korolev et al., 2011;O'Shea et al., 2020). Even modern PSDs based on large amounts of aircraft data, like Heymsfield et al. tation of the "cloud ice" category in global models (Geer, 2021b). As with the simple exponential PSD now used as the default for cloud ice, it may in some cases be better to infer PSDs through parameter estimation, using the parameters of the standardised Modified Gamma Distribution (MGD, Petty and Huang, 2011).
-A representation of non-spherical melting hydrometeors (e.g. Johnson et al., 2016;Leinonen and von Lerber, 2018) needs to be included. This will be particularly important for simulating radar backscatter in the melting layer (the bright band).

825
It would also improve brightness temperature simulations, which may otherwise be too low in channels and situations sensitive to the melting layer, possibly by up to 8 K (Bauer, 2001). This survey has also revealed a few issues more specific to the RTTOV hydrotable generator: -The lowest temperature bin, at 203 K, is at least 20 degrees higher than the lowest tropospheric temperatures. Hence, in a future version, it is necessary to extend the lower temperature range of frozen particles permitted by the table generator.

830
Note that, although the lowest temperature point in the ARTS database is 190 K, Eriksson et al. (2018) showed that extrapolation would be reasonably accurate down to 150 K.
-The default representation of rain has not been updated since Bauer (2001); apart from the aforementioned points on shape and orientation, it is worth examining how well the rain PSD is represented, and potentially updating this.
-The code is operated offline, with a lookup table approach, but it is fast enough to be operated online within the radiative 835 transfer code, especially as weather forecasting systems get more capable (e.g. English et al., 2020). This will be an aim of future development work.
To summarise, the future evolution of the code may be towards online operation, so that parameters of the PSD and even particle shape can be estimated directly from the observations, rather than prescribed in a "one shape fits all" approach.
Code and data availability. This work is based on RTTOV v13.0, which is available for free from https://nwp-saf.eumetsat.int/ to users who 840 have registered and agreed to the licence conditions. The licence is available at https://nwp-saf.eumetsat.int/site/software/licence-agreement/ with the main condition being that onward redistribution of the code is not permitted. For the purposes of anonymous review of the current manuscript, the code of the hydrotable generator is temporarily available at http://nwp-saf.eumetsat.int/downloads/james/rttov13_ hydrotable_gen_for_review_only.tar.gz. This code should not be redistributed or used for any other purpose.
Author contributions. All authors contributed to the code, science and data embedded in the RTTOV-SCATT hydrotable generator. AG wrote 845 the majority of the text, with reviews and contributions from all other authors.
Competing interests. No competing interests are present permittivity model (see Sec. 2.1). If these configuration choices are not sufficient, users would need to directly change the code according to the way the radar reflectivity factor is defined for any particular instrument; this might need to be improved in future.
Appendix B: Two-stream slab cloud with scattering and thermal emission and absorption The two-stream approximation for radiative transfer is proposed here as a way to compare different sets of bulk optical proper-865 ties. To simplify the unpolarised full scattering radiative transfer equation, the two-stream approximation assumes the radiance field is constant in each hemisphere (as a function of azimuth and zenith angle) and is hence described by just two variables, the upward and downward radiances I ↑ and I ↓ (Thomas and Stamnes, 1999;Petty, 2006). A second assumption is that backscattering from one hemisphere into another is proportional to the asymmetry parameter g. With these assumptions, the following radiative transfer equations can be defined, which are averaged over all azimuth and zenith angles in each hemisphere: This follows the derivation in Petty (2006) but does not drop the Planck function B, so it represents thermal emission as well as scattering. The single scattering albedo ω 0 , asymmetry g and temperature (and hence Planck function B) are assumed constant Here, A and D are constants deriving from c 4 to c 4 that still need to be determined. r ∞ is the cloud albedo, a function of ω 0 and asymmetry g defined earlier. The constants are found by imposing the boundary conditions I ↓ (τ * ) = 0 and I ↑ (τ * ) = I 0 , in other words zero downward radiation at the top of the cloud, and a fixed amount of radiation upwelling from below the cloud.
The boundaries are hence assumed black bodies, unaffected by the radiation emitted from the cloud. This gives: These can be easily solved for A and D, and the radiation emitted at the top of the cloud I ↑ (0) can be determined from Eq. B7 to provide Eq. 20. Setting B = 0 recovers the two-stream pure-scattering example from Petty (2006, his equations 13.39 and 13.40) . Setting I 0 = 0 and g = 0 (no external radiation sources, isotropic scattering) recovers the "imbedded source" solution from Thomas and Stamnes (1999, their equations 7.69 and 7.70).