Improved scattering radiative transfer for frozen hydrometeors at microwave frequencies

To simulate passive microwave radiances in allsky conditions requires better knowledge of the scattering properties of frozen hydrometeors. Typically, snow particles are represented as spheres and their scattering properties are calculated using Mie theory, but this is unrealistic and, particularly in deep-convective areas, it produces too much scattering in mid-frequencies (e.g. 30–50 GHz) and too little scattering at high frequencies (e.g. 150–183 GHz). These problems make it hard to assimilate microwave observations in numerical weather prediction (NWP) models, particularly in situations where scattering effects are most important, such as over land surfaces or in moisture sounding channels. Using the discrete dipole approximation to compute scattering properties, more accurate results can be generated by modelling frozen particles as ice rosettes or simplified snowflakes, though hexagonal plates and columns often give worse results than Mie spheres. To objectively decide on the best particle shape (and size distribution) this study uses global forecast departures from an NWP system (e.g. observation minus forecast differences) to indicate the quality of agreement between model and observations. It is easy to improve results in one situation but worsen them in others, so a rigorous method is needed: four different statistics are checked; these statistics are required to stay the same or improve in all channels between 10 GHz and 183 GHz and in all weather situations globally. The optimal choice of snow particle shape and size distribution is better across all frequencies and all weather conditions, giving confidence in its physical realism. Compared to the Mie sphere, most of the systematic error is removed and departure statistics are improved by 10 to 60 %. However, this improvement is achieved with a simple “one-size-fits-all” shape for snow; there is little additional benefit in choosing the particle shape according to the precipitation type. These developments have improved the accuracy of scattering radiative transfer sufficiently that microwave all-sky assimilation is being extended to land surfaces, to higher frequencies and to sounding channels.


Introduction
Microwave observations are widely used to infer atmospheric temperature and water vapour, particularly in numerical weather prediction (NWP, e.g.English et al., 2000).Increasingly, NWP centres are making use of these observations in cloudy and precipitating situations as well as in clear skies (e.g.Bauer et al., 2011).This helps to infer water vapour information in cloudy and precipitating areas and it also gives the possibility to assimilate the cloud and precipitation itself.When all observations, whether clear, cloudy or precipitating, are assimilated using the same scattering-capable radiative transfer model, this is often referred to as an "all-sky" approach (e.g.Bauer et al., 2010).However, it has been difficult to use cloud-and precipitation-affected microwave observations in situations where atmospheric scattering is most important, such as over land surfaces and in temperature and water vapour sounding channels (e.g.Baordo et al., 2012;Geer et al., 2012).This study aims to improve the quality of radiative transfer for NWP by improving the modelling of frozen hydrometeor optical properties.Observation minus forecast statistics from an NWP system will be used to objectively guide the choices of frozen hydrometeor particle model.
Optical properties of single particles have typically been estimated using Mie theory, with an ice or snow particle represented as a sphere containing a mixture of ice and air.It has often been necessary to tune the choice of particle size distribution and the sphere's density (often also a function of particle size) to improve the quality of simulations (e.g.Wiedner et al., 2004;Doherty et al., 2007;Sreerekha et al., 2008).Though tuning schemes can be effective at one frequency, they can make results worse at others.To obtain good results from Mie theory Surussavadee and Staelin (2006) went as far as representing frozen particles with spheres that had a different size and density at each frequency.It would be preferable to use a particle model with a closer link to physical reality.The discrete dipole approximation (DDA, Purcell and Pennypacker, 1973;Draine and Flatau, 1994) is becoming more widely used in microwave applications.The DDA represents a particle as a three-dimensional array of polarisable points and provides a better model of the optical properties of non-spherical particles than can be obtained from the Mie sphere (Kulie et al., 2010;Petty and Huang, 2010) or indeed from homogeneous spheroids (Leinonen et al., 2012).Databases are available that contain the pre-computed DDA optical properties of idealised ice and snow particles (e.g.Liu, 2004Liu, , 2008;;Hong, 2007;Hong et al., 2009), making it practical to incorporate discrete dipole results into the fast radiative transfer models required for data assimilation.
Mie spheres produce their worst results in areas of deep convection.Figure 1 shows Hurricane Irene on 25 August 2011 at four microwave frequencies, with the radiances represented in terms of brightness temperature (TB).The first column shows observations from TMI (TRMM Microwave Imager) and SSMIS (Special Sensor Microwave Imager Sounder); the second and third columns show simulations from the European Centre for Medium-Range Weather Forecasts (ECMWF) first-guess (FG) using either Mie sphere or DDA models for snow.At 10 GHz, where the observations sense principally rain emission, the simulations and observations are in reasonable agreement, though the model captures only the central core of the hurricane and does not capture the full intensity of the rain band to the north.At higher frequencies the Mie simulations are badly wrong.At 37 GHz the observations reveal the hurricane as an area of warm brightness temperatures in the range 260 K to 280 K, which could be achieved by emission (predominantly from water cloud) at an altitude of about 5 km.Instead, the Mie simulations show the hurricane as a "black hole" with a radiant temperature of around 210 K.In the Mie simulations, frozen hydrometeors cause an excessive depression in brightness temperatures, i.e. excessive scattering.By contrast, at 150 GHz, where scattering from upper-level ice and snow is expected to depress TBs (e.g.Hong et al., 2005), the Mie simulations do not provide enough scattering.The DDA simulations are based on the Liu (2008) sector snowflake, a shape that will later be identified as optimal by this study.Overall, moving to DDA has removed a lot of the excessive scattering at middle frequencies (i.e.TBs have become higher at 37 GHz and 52.8 GHz) and increased scattering at high frequencies (i.e.TBs have become lower at 150 GHz).However, problems with the frequency dependence of scattering from the Mie sphere are not limited to tropical convective areas.Kim et al. (2007) looked at winter light precipitation in the midlatitudes and found that a Mie sphere with the same physical parameters was unable to provide good results simultaneously at 89 GHz and 150 GHz.
Even if we can abandon unrealistic particle models like the Mie sphere, the problem remains that the particle sizes and shapes are poorly known and subject to enormous natural variability.While there have been observational studies of particle shapes and size distributions, they only represent case studies and they cannot provide sufficient guidance for an objective and globally applicable description of hydrometeor shapes and sizes.Also there have been many simulation studies, often showing apparently good agreement with observations.However, the quality of agreement has not always been quantified and the results have again been limited typically to case studies -such as a single midlatitude front -and to a small range of microwave frequencies.A rare example of a study with global, broad-frequency applicability is that of Kulie et al. (2010), who looked at constraining the choice of particle shape for DDA scattering computations using a combination of radar and passive microwave observations.Ultimately they did not recommend any particular shape because of the great variability of hydrometeor habits in the atmosphere.
In the current study we hope to find simple models for scattering from frozen particles that can improve or maintain agreement between model and observations across all weather conditions and all microwave frequencies from 10 GHz to 183 GHz.A rigorous methodology will be applied to quantify the fit between observations and model and to make sure this is improved everywhere: to make sure that by improving one aspect, we are not degrading another.Global, broad-frequency applicability is a necessity for NWP and moreover it should give greater confidence in the physical basis of the chosen particle models.
Although in an ideal world the choice of shape and size distribution would be situation dependent, this would add complexity to the radiative transfer model and make it harder to objectively validate these choices.It would be better to have an objectively tuned simple scheme than a complex one that sounded realistic but was not properly tuned.Though the use of a single size distribution and shape is an oversimplification, the resulting errors need not cause serious problems in the data assimilation context.As long as the errors are random rather than systematic, the poorer accuracy of cloud and precipitation radiative transfer can be accounted for with an observation error model that assigns bigger errors in cloudy and precipitating situations than in clear skies (e.g.Geer and Bauer, 2011).Furthermore, forecast models find it difficult to put cloud and precipitation in exactly the right place with the right intensity (e.g.Roberts and Lean, 2008;Fabry and Sun, 2010).In practice, the observation errors assigned in all-sky data assimilation are very large, reaching 20-40 K in convective situations (Geer and Bauer, 2011).These observation errors represent both radiative transfer error and "mislocation" error.In the presence of such large random errors, a radiative transfer error of 2 K is irrelevant and the real concern is with large systematic errors greater than 20 K, such as those in Hurricane Irene in Fig. 1.At this stage in the development of all-sky data assimilation, it is the most obvious problems that need fixing.As will be shown, the finer details of many cloud and precipitation radiative transfer issues do not matter when comparing current forecast models with real observations.The data assimilation system, radiative transfer model and observations are introduced in Sect. 2 and the computation of bulk optical properties is described in Sect.3. Methods for using the fit between NWP model and observations to find the best size distribution and particle shape are considered in Sect. 4. The results are presented separately for ocean and land surfaces in Sects.5 and 6.  2007) size distribution.Some Liu (2008) shapes have similar bulk scattering properties to others and are ignored: the five-bullet rosette is similar to the six-bullet; the short column is similar to the thick plate.

Radiative transfer model
Radiative transfer simulations are provided by RTTOV-SCATT, which is a fast model designed for assimilating microwave radiances in all-sky conditions (Bauer et al., 2006).It is a component of the wider RTTOV package (Radiative Transfer model for Television Infrared Observation Satellite Operational Vertical sounder; Eyre, 1991;Saunders et al., 2012).The radiative transfer equation is solved using the delta-Eddington approximation (Joseph et al., 1976).Transmittances for oxygen and water vapour are computed from regression tables driven by atmospheric predictors, just as in the normal RTTOV.Bulk optical properties for cloud water, cloud ice, rain and snow are taken from look-up tables that will be described in Sect.3.1.Ocean surface emissivity is computed by version 5 of FASTEM (English and Hewison, 1998;Liu et al., 2011;Bormann et al., 2012).Landsurface emissivity comes from the TELSEM atlas (Aires et al., 2011).The all-sky brightness temperature is computed as the weighted average of the brightness temperature from two independent sub-columns, one clear and one cloudy.The weighting is done according to the effective cloud fraction of Geer et al. (2009a) which provides a fast but approximate way to account for the effects of sub-grid variability in cloud and precipitation, particularly the beam-filling effect (e.g.Kummerow, 1998).

ECMWF system
ECMWF produces global forecasts and analyses using a 4D-Var data assimilation system (Rabier et al., 2000).Microwave imager radiances are assimilated directly in 4D-Var alongside many other conventional and satellite observation types.Observations drive the data assimilation system through the first guess (FG) departure d, which is the difference between real and simulated observations y o and y b : A simulated observation is computed as The background x b (t 0 ) is a forecast initialised from the previous analysis, with t 0 being the time of the start of the assimilation window.The nonlinear forecast model M[ ] propagates this atmospheric state forward in time.In this paper "first guess" refers to the complete forecast trajectory defined by M[x b (t 0 )] through the 12 h assimilation window.H [ ] is the observation operator which, in the case of microwave imager observations, selects the nearest model profile to the observation (in time and space) and then runs RTTOV-SCATT.
A bias correction c is included in the computation of the FG departures.This is essential to remove the systematic differences between simulations and observations that result from a combination of instrument, observation operator and forecast model biases.For microwave instruments, biases are inferred as functions of predictors including the scan angle, the surface wind speed and the layer thickness, though the exact set of predictors is channel dependent.Bias coefficients are derived within the analysis system using variational bias correction (VarBC, Dee, 2004;Auligné et al., 2007).There are no cloud-related predictors and the bias correction is not intended to represent cloud-or precipitation-dependent biases.The convention in this paper is to consider the bias correction part of the simulation.
Further details of the all-sky microwave imager assimilation at ECMWF are given by Bauer et al. (2010), Geer et al. (2010b) and Geer andBauer (2010, 2011).For assimilation, a wide range of quality control measures need to be applied, but a smaller set of restrictions will be applied here: observations are restricted to latitudes equatorward of 60 • ; scenes containing sea ice or coasts are removed; the surface temperature must be higher than 274 K over ocean and 278 K over land to help avoid sea ice and snow cover.Unlike in Fig. 1, observations are averaged (or "superobbed") in boxes of approximately 80 km by 80 km, in order to make the horizontal scales of observed cloud and precipitation more similar to their effective resolution in the model; the sensitivity of the results to this choice is examined later.
The ECMWF model has four prognostic hydrometeor types representative of large-scale cloud processes: cloud water, cloud ice, rain and snow; also cloud and precipitation fractions are provided on each model level (Tiedtke, 1993;Forbes et al., 2011).In addition to this, the convective rain and snow on each model level is diagnosed from a convection scheme which assumes that only 5 % of each grid box contains convection.There is no explicit representation of convective cloud.For input to RTTOV-SCATT, precipitation fluxes are converted to mixing ratios by assuming a distribution of particle sizes and fall speeds consistent with the computations in RTTOV-SCATT (Bauer, 2001).Convective and large-scale precipitation are added together.Hence the hydrometeor inputs to RTTOV-SCATT are the vertical profiles of cloud water, cloud ice, total rain and total snow, plus the effective cloud fraction, which is a hydrometeor-weighted average of the cloud, precipitation and convective fractions across all vertical levels (Geer et al., 2009b).This study also looks at splitting the snow category into large-scale and convective hydrometeor types to better account for differences in their microphysical characteristics.
To provide a set of FG model fields M[x b (t 0 )] and VarBC bias corrections c, the full assimilation system has been run from 1 to 30 June 2012 using 91 levels in the vertical and a horizontal resolution of approximately 40 km (T511 in spectral terms).In the following experiments, only the radiative transfer model H [ ] is varied when computing the departures d (Eqs. 1 and 2).Cycle 38r2 of the ECMWF system has been used.

Observations
In order to investigate frequencies from 10 GHz to 183 GHz, this study combines observations from TMI and SSMIS.Table 1 summarises the channels used.
TMI (Kummerow et al., 1998) on the Tropical Rainfall Measuring Mission (TRMM) has a relatively high spatial resolution (between 60 km and 6 km, depending on the channel) and channels between 10 GHz and 85 GHz.TRMM's inclined orbit was designed to sample the entire diurnal cycle of tropical precipitation, which limits the coverage to a band between roughly 40 • S and 40 • N. TMI observations have been obtained from NASA and are at version 6. Solardependent anomalies are present in the data (e.g.Gopalan et al., 2009) but they are accounted for in the ECMWF bias correction (Geer et al., 2010a).
SSMIS (Kunkee et al., 2008) has a slightly lower spatial resolution than TMI and no 10 GHz channel, but instead it provides temperature sounding channels in the 50 GHz oxygen line and moisture sounding channels in the 183 GHz water vapour line.Though there are a number of satellites with an SSMIS onboard, only Defense Meteorological Satellite Program satellite F17 (DMSP-F17) has been used, in line with ECMWF operational usage.The data have been preprocessed to eliminate calibration anomalies following Bell et al. (2008).There are still some anomalies of order 0.2 K visible in the FG departures in the 50 GHz channels.These anomalies are one of the main reasons the SSMIS 50 GHz channels are not assimilated at ECMWF.However, as will be demonstrated later, these anomalies are not large enough to affect the results of this study, with its focus on scattering signals of order 20 K.

Computation
To solve the radiative transfer equation, RTTOV-SCATT needs to know the bulk optical properties of the atmosphere at each model level.Given the optical properties of a single particle as a function of its maximum dimension D, i.e. the diameter in the case of a sphere, bulk scattering properties are computed by integrating across the size distribution, N (D).From the extinction and scattering cross sections σ e (D) and σ s (D) and the asymmetry parameter g(D), it is possible to compute the extinction coefficient β e , scattering coefficient β s , and average asymmetry parameter g bulk : (3) In practice, RTTOV-SCATT represents the scattering coefficient through the single-scattering albedo (SSA), ω o = β s /β e .
To determine the size distribution N (D) from the hydrometeor water content it is necessary to know the particle mass as a function of its maximum diameter, m(D) = aD b . (6) In the case of a sphere with constant density, for example, b = 3 and a is determined by the particle density.The water content (hydrometeor mass per unit volume) is given by It is then necessary to define a functional form for the size distribution N (D).One parameter in the size distribution is left free so that it can be adjusted according to the water content.In other words, given l, a and b, the free parameter of the size distribution can be determined.Size distributions are usually designed in a way that makes for a convenient analytic solution (e.g.Appendix A; Marshall and Palmer, 1948;Field et al., 2007;Petty and Huang, 2011).A common way to look at Eq. ( 7) is to see that the water content defines the bth moment of the size distribution.
Since the generation of the bulk optical properties is computationally demanding, they are pre-tabulated for each hydrometeor type as a function of temperature, frequency and water content.Given the water content of each hydrometeor type present in the layer, the final bulk optical properties of the layer are obtained by summing over hydrometeor types in a manner analogous to the integrals in Eqs.(3)-( 5).More information can be found in Bauer (2001), but note that the melting layer effects described in that paper are not included here.
Up until now in RTTOV-SCATT, all hydrometeors have been modelled as spherical particles using Mie theory.Cloud water and cloud ice have been modelled using a gamma size distribution (e.g.Petty and Huang, 2011) and a constant density.Rain and snow have used a Marshall and Palmer (1948) size distribution and again a fixed density.See Appendix A for details.Frozen particles are assumed to be made up of ice inclusions in an air matrix, with the dielectric properties combined according to the approach of Bohren and Battan (1982) and Fabry and Szyrmer (1999).Replacing this with the standard method of Maxwell-Garnett (1904) makes very little difference to the simulated brightness temperatures, and it is unlikely that other schemes, such as those discussed by Petty and Huang (2010), would make that much difference either.
The authors have added into RTTOV-SCATT the facility to use optical properties for non-spherical hydrometeors from the Liu (2008) database.This is available with RTTOV version 11.In the Liu database, optical properties are tabulated as a function of frequency, temperature and particle size for a variety of hexagonal ice columns and plates, rosettes composed of between three and six orthogonal "bullets" and two simple snowflake models, the "sector" and "dendrite".See Liu (2008) for further information.The particles have been assumed to be randomly oriented, so the optical properties are the average over a large number of random orientations.In the computation of the bulk optical properties for RTTOV-SCATT, the particle shape determines the a and b coefficients of the mass-size relation (Eq.6).Coefficients a and b appropriate to the Liu shapes have been taken from Table 1 of Kulie et al. (2010).Some minor issues around that choice are described in Appendix B.
To simulate bulk optical properties from the Liu (2008) shapes, the Field et al. (2007) size distribution has been chosen.It (or its predecessors) have been a typical choice in recent studies (e.g.Doherty et al., 2007;Kulie et al., 2010;Di Michele et al., 2012).This size distribution exists in tropical and midlatitude versions; for simplicity the tropical version has been used globally.In this study, where we are searching for a more optimal model for scattering properties, we make the deliberate decision to keep the size distribution fixed while varying the particle shape.Things would have become far too complex if the size distribution were also allowed to vary.However, sensitivity to the choice of size distribution will be examined later.
A final point of detail can be made.Integrals like Eqs. ( 3), ( 4) and ( 5) are in practice computed numerically and the integration range is truncated rather than running from 0 to infinity.We do not believe this should affect the results substantially, but further information is given in Appendix C.

Comparison
Figure 2 shows the bulk optical properties of the Liu (2008) shapes with the Field et al. (2007) size distribution, at 52.8 GHz as a function of snow water content.Mie sphere results are also shown, both for the Marshall-Palmer and Field et al. (2007) size distribution.Plates and columns tend to have an SSA and extinction larger than or about the same as the Mie sphere.Bullet rosettes and snowflakes tend to have much lower SSA and extinction.Only by looking at the asymmetry can the Liu shapes be clearly distinguished from the Mie sphere.Mie theory gives very strong forward scattering for larger size parameters, i.e. asymmetry approaching 1, so high snow water contents have strong forward scattering.The Liu DDA computations produce more balanced forward and backward scattering, i.e. asymmetry is always in the range 0 to 0.2. Figure 3 looks at the variation with frequency and shows that these conclusions are broadly true for all frequencies between 50 GHz and 183 GHz.However, at lower frequencies the Mie sphere with Marshall-Palmer size distribution produces more scattering and extinction than any Liu shape, consistent with the excessive scattering exhibited by the Mie simulations in Fig. 1 at 37 GHz and 52.8 GHz.
At ECMWF, before starting to work with the DDA shapes of Liu (2008), a number of attempts were made to improve the accuracy of the radiative transfer model while continuing to represent snow hydrometeors using Mie spheres, taking inspiration from other studies (e.g.Wiedner et al., 2004;Surussavadee and Staelin, 2006;Doherty et al., 2007;Johnson et al., 2012).However, the ECMWF attempts were not successful and there seems little point in detailing them here; as mentioned in the introduction, improving results at one frequency gives greater problems at another.To give just one example, using the Field et al. (2007) size distribution with the Mie sphere would help moderate the amount of scattering at lower frequencies but, as can be inferred from the reductions in extinction and SSA in Fig. 3, it makes the problem of under-scattering at high frequencies even worse.The Mie sphere results that follow will be based on the Marshall-Palmer size distribution and the fixed 100 kg m −2 density, so that comparisons can be made directly to the old RTTOV-SCATT.
4 Choosing the best DDA shape for NWP

Overview
As explained previously, the aim is to find the best Liu (2008) particle shape or shapes using observation minus forecast statistics as an objective criterion.The search is split into a coarse search and a fine search.The coarse search concentrates on the main issue, which is the representation of scattering from snow hydrometeors.It seeks one globally applicable DDA shape to represent the snow hydrometeor category.Cloud ice is less important in the radiative transfer for two reasons.First, the ECMWF model produces much less cloud ice than snow.For example, at a grid point in the core of Hurricane Irene, the model produces roughly 40 kg m −2 of snow but only 5 kg m −2 of cloud ice.Second, ice particles are typically smaller than snow particles and hence are less effective scatterers.For simplicity in the coarse search, the optical properties of cloud water, cloud ice and rain will be held fixed and will continue to be represented by Mie spheres.The fine search then looks for a DDA shape to represent cloud ice, and investigates the use of separate DDA shapes for convective and large-scale snow.
A further issue in the data assimilation context is the likelihood of biases in the moist physics of the forecast model, which can lead to systematic errors in modelled cloud or precipitation.Cloud-related biases have so far proven very difficult to control using bias correction schemes, partly because it is hard to find simple predictors to describe complex, situation-dependent biases (e.g.Geer and Bauer, 2010).The pragmatic solution to choosing a particle shape for data assimilation is to find the one that leads to the smallest differences between simulations and observations in a global, statistical sense, acknowledging that these choices may be compensating for other forecast model or radiative transfer model biases.

Quantifying the fit of model to observations
Typically the fit between model and observations is quantified using the mean, standard deviation or rms of FG departures d (Eq.1).However, for observations sensitive to cloud and precipitation, the standard deviation and rms are affected by the well-known "double-penalty" effect: to achieve a small rms, it is better not to forecast cloud and precipitation at all, than to forecast it at the wrong time or in the wrong place.As discussed in the Introduction, cloud and precipitation are not predictable on small scales in current NWP models so the departures d are affected by mislocation errors.Hence it would be misleading to rely solely on the standard deviation or rms of d when choosing the best DDA particle shape.Here we introduce a number of alternatives that are resistant to the double penalty effect.
Figure 4 shows histograms of FG departures at 52.8 GHz computed using the Mie sphere and two of the best DDA shapes in this study (three-bullet and sector snowflake).The majority of departures are small and are associated with clear and cloudy situations.In fact, 99.7 % of Mie sphere departures are smaller than 10 K in an absolute sense.The logarithmic y axis helps focus on the small proportion of scenes with poorly simulated precipitating situations, i.e. those that cause large FG departures.A positive departure indicates cases where brightness temperatures are lower in the simulation that in the observation.Since the scattering signal dominates at 52.8 GHz, these are generally situations where either the forecast model generates more snow than is observed or the radiative transfer model simulates excessive scattering.An unbiased model would produce a symmetric histogram because there would be an equivalent set of negative departures coming from situations where there is more snow in the observations than in the model.Following this reasoning, the Mie simulations must be producing too much scattering.This can of course be more easily inferred from Fig. 1.The point of Fig. 4 is to make it easier to determine which of the DDA simulations is best: the three-bullet appears to produce slightly too much scattering at this frequency and the sector snowflake too little.To quantify this, the skewness can be used.The skewness has strong sensitivity to outliers (see e.g.Wilks, 2006) but this is a desirable property when we are looking for the large but infrequent errors associated with snow-scattering situations.The skewness is −2.6, +2.4 and 7.7 for the sector snowflake, three-bullet and Mie sphere in Fig. 4. On this measure, and for this frequency, the threebullet rosette is best.
An alternative way to avoid the double penalty issue is to compare histograms of observed and simulated brightness temperature (e.g.Doherty et al., 2007).Figure 5a shows these histograms at 52.8 GHz.It is again obvious that the Mie sphere simulations are badly wrong, because they produce brightness temperatures as low as 185 K when observations never go lower than about 220 K.The three-bullet is better, in that it generates TBs no lower than 220 K, but compared to observations there are too many occurrences in bins at 222.5 K and 227.5 K.The sector snowflake does not produce TBs lower than 230 K, which is perhaps equally undesirable.
To measure the consistency between two distributions in a statistical sense, it is common to use either the Chi squared or Kolmogorov-Smirnov tests (χ 2 or K-S, see e.g Wilks, 2006) or the Kullback and Leibler (1951, K-L) divergence.However, the K-S and K-L tests are unsuitable because they are insensitive to the parts of the histogram with small populations.Further, the χ 2 test does not work well when there are small or non-existent populations in the observed histograms.Hence, the following statistical measure is proposed: This statistic penalises discrepancies between histograms using the log of the ratio of populations in each bin (i.e the number simulated divided by the number observed).It is similar to the K-L divergence but crucially it does not weight the penalty in each bin by the number of observations in each bin.The measure in Eq. ( 8) becomes infinite when either of the bin populations is zero.To prevent this, empty bins have been assigned a population of 0.1.This number can be tuned to give a greater or smaller penalty to situations with empty bins; 0.1 seemed to give a good balance.
Figure 5b shows the log ratios corresponding to the three histograms in Fig. 5a.The unphysically low TBs produced by the Mie simulations are penalised by log ratios of +1 to +2 in the bins below 220 K.In the bins at 222.5 K and 227.5 K, the three-bullet produces too many occurrences so each bin is penalised with a positive log ratio.In contrast, the sector snowflake predicts no occurrences in this range so it is penalised with negative log ratios around −1.In the more common range of brightness temperatures in this channel (230 K to 265 K), observations and simulations agree well, so the penalties are small.There is also a slight overestimate in the number of simulated brightness temperatures in the highest populated bin at 267.5 K, which results in a log ratio of +0.7, but this cannot be associated with snow scattering since it affects all three simulations equally.To complete the statistic, the absolute value of the log ratio is summed across all bins and divided by the number of bins in which observations occur.This statistic produces values of 0.27, 0.23 and 1.54 for the sector snowflake, three-bullet and Mie sphere in Fig. 5.
A final measure of fit is based on maps of mean FG departures like those shown in Fig. 6.Here, the excessive scattering produced by the Mie sphere is indicated by a band of positive departures along the Intertropical Convergence Zone (ITCZ), i.e. in areas where deep convection is frequent.There is a smaller band of positive departures along the SH storm tracks at around 40 • S. The sector snowflake completely eliminates these features, whereas the three-bullet rosette still produces slightly too much scattering in the tropics since some areas of positive bias remain.This can be quantified by calculating an rms across all latitude-longitude bins in which a mean FG departure has been computed.In this example, the values of the rms are 0.17 K, 0.22 K and 0.45 K for the sector snowflake, three-bullet and Mie sphere.
Figure 6 exposes an issue with the DMSP-F17 SSMIS observations in the 50 GHz temperature sounding channels.Improving the snow-scattering model has enhanced the visibility of a band of negative departures across the NH at 30 • N.With the Mie sphere (Fig. 6a) these negative departures were visible mainly on the E sides of the N Pacific and N Atlantic, regions of climatologically little deep convection, i.e. areas less affected by the problems with the Mie sphere.With the sector or three-bullet (Fig. 6b and c) a band of around −0.2 K encircles the globe at this latitude.This is probably due to the previously mentioned issues with the SSMIS measurements, i.e. solar heating or solar intrusion effects.However, with a magnitude around 0.2 K, the bias at 30 • N can have no real effect on the FG departure histograms or the histogram fit, where the statistics respond to differences in brightness temperatures of order 10 K (Figs. 4 and 5).Using a variety of ways to measure the fit to observations makes the study more robust against issues like this.   5 Results over ocean

Coarse search
Results are presented for the Mie sphere and the four DDA shapes that produce the most realistic simulations: three-bullet and six-bullet rosettes and dendrite and sector snowflake.As will be seen from the statistics, these shapes fall either side of the ideal fit, scattering either slightly too much or slightly too little (as illustrated in Fig. 4).Also examined is the thin plate, which is the next available DDA shape producing more scattering than the three-bullet (Fig. 2).The poor results from the thin plate serve to illustrate that the appropriate choice of DDA shape (and/or size distribution) is crucial.Figure 7 summarises the statistics of fit across all channels.Panel (a) shows the rms of FG departures as a percentage of those from the Mie sphere simulations.The DDA thin plate results are substantially worse than the Mie sphere and have been allowed to go off scale.The four best DDA shapes provide reductions of 30-40 % in the 50.3GHz and 52.8 GHz channels and smaller reductions at 37 GHz and 53.6 GHz.However, in the higher frequencies, particularly in the 183 ± 7 channel, the six-bullet, sector and three-bullet increase the rms of the FG departures.These shapes produce more scattering than the Mie sphere, but that is a good thing at these frequencies (see Fig. 1), so the increase in rms must come from the double penalty issue.
The skewness of the departure histogram is shown in Fig 7b .The most obvious feature is the positive skewness of the Mie sphere departures in most channels from 19 GHz to 53.6 GHz and the negative skewness at 150 GHz and 183 GHz.In other words, scattering is excessive at low frequencies and insufficient at high frequencies.It is still tricky to find a DDA shape that has minimal skewness, i.e. an appropriate amount of scattering, at both high and low frequencies.The six-bullet has little skewness at 52.8 GHz and 53.6 GHz, but negative skewness at 150 GHz and 183 GHz.In contrast, the three-bullet does alright at 150 GHz and 183 GHz, but produces too much scattering, i.e. positive skewness, at 52.8 GHz and 53.6 GHz.Some DDA shapes are poor at all frequencies: the thin plate always produces excessive scattering; the dendrite too little.The sector snowflake is slightly under-scattering as has already been illustrated at 52.8 GHz (Sect.4.2) but it provides consistent results across the frequencies.
It is initially strange to see that snowflakes and rosettes produce more scattering than the Mie sphere at high frequencies (e.g.Fig 7b).Referring back to the bulk scattering properties in Fig. 3, snowflakes and rosettes produce lower SSA and extinction than the Mie sphere.However, the explanation is in the asymmetry parameter: Mie spheres produce far stronger forward scattering than any of the DDA shapes.If most scattering is in the forward direction, much more radiation from warmer lower levels will get through to the sensor; in other words scattering is less effective at causing brightness temperature depressions.This explains why the DDA shapes are able to avoid the problem of excessive scattering at lower frequencies while still providing enough brightness temperature depression at higher frequencies.It is because they do not generate such intense forward scattering as the Mie solution.
The histogram discrepancy statistic (h, Eq. 8) is shown in Fig. 7c.The Mie sphere, the thin plate and the dendrite are the worst by this measure, which penalises both excessive scattering (thin plate, Mie sphere at lower frequencies) and insufficient scattering (dendrite, Mie sphere at higher frequencies).The three-bullet, six-bullet and sector snowflake produce discrepancies less than 0.5 across most of the frequency range, with no obvious "best" shape.This means that all three produce reasonably physical distributions of brightness temperature.
One feature of note in Fig. 7c is the behaviour of the histogram discrepancies at 10 GHz.All DDA shapes are marginally worse than the Mie sphere.If the results are affected at 10 GHz, this suggests that the Mie sphere was generating scattering from snow hydrometeors at unphysically low frequencies.Much more work would be required to investigate properly, but the likelihood is that unphysical scattering from snow at 10 GHz was compensating for another bias in the model.
Figure 7d shows the rms of latitude-longitude mapped biases.As for Fig. 7a, the values are given as a percentage of the Mie sphere rms FG departures.Again this helps to normalise the biases according to the brightness temperature errors in each channel.The 10 GHz channels stand out by this measure: monthly mean biases are order 30 % of the rms of FG departures.In other words, bias is quite large compared to the signal in these channels; the large uncorrected biases are one main reason preventing operational assimilation of the 10 GHz channels at ECMWF.Ignoring the 10 GHz channels, the Mie sphere, the thin plate, the dendrite and the three-bullet all produce biases greater than 20 % in some channels.Figure 6 has already illustrated the situation at 52.8 GHz, where excessive scattering from the Mie sphere causes large biases in the ITCZ.At higher frequencies the three-bullet is worse than the Mie sphere, despite appearing a strong candidate for "best" particle shape according to the other measures.The bias maps for these higher frequencies (e.g.Fig. 8; others not shown) show that biases in the tropics are successfully minimised by the relatively strong scattering from the three-bullet, but this amount of scattering is excessive in the midlatitudes.The sector snowflake is a compromise which produces slightly too little scattering in the tropics and slightly too much in midlatitudes.There are certainly limits to the "one-shape-fits-all" strategy.In Table 2, statistics have been aggregated across all channels.This is done by computing the mean across all channels of the statistics shown in Fig. 7.An exception was the skewness from Fig. 7b, where an rms across all channels is a more appropriate way of aggregating the data.By these measures, the sector snowflake, three-bullet and six-bullet are all much better than the Mie sphere.Though the sector snowflake produces slightly too little scattering in tropical convection it gives consistently good results by all four measures of fit.

Fine search
The fine search considers three categories of frozen hydrometeor.Cloud ice, previously simulated using a Mie sphere, is instead simulated using a DDA shape, and the snow hydrometeor category is split into a convective part and a large-scale part according to which model parameterisation produced the precipitation.Table 3 lists the shapes that were tried in each hydrometeor category, listed roughly in order of their scattering ability.To keep the number of combinations within practical limits, only two cloud ice shapes were tried, yielding a total of 24 separate experiments.As illustrated by Fig. 8, the sector snowflake produces excessive scattering at 183 GHz in the midlatitudes, but insufficient in the tropics.Hence, the six-bullet and dendrite shapes were tested as alternatives in the large-scale snow category with the intention of reducing scattering at midlatitudes.Conversely, the thin plate and block column were tried in the convective snow category with the intention of increasing scattering in the tropics.For cloud ice, sector and dendrite snowflakes were tried.Using snowflakes to represent cloud ice might sound unphysical, but exploratory tests found that representing cloud ice with the thin hexagonal plate caused excessive scattering (i.e.too-low TBs) in midlatitude frontal cloud.Rather, the best results for cloud ice were to be found with low-scattering shapes like the snowflakes.
It is hard to visualise the results of the search in three dimensions, so Fig. 9 ranks the results on a scale where the worst fit in each of the four statistics is normalised to one.There are 26 experiments included: all 24 combinations plus the Mie sphere and sector snowflake experiments from the coarse search.The average ranking across all four scores is given in Table 4.Only a few illustrative experiments are identified in the figure.The Mie sphere (diamond symbol) is by far the worst in terms of skewness and histogram fit and it is among the worst in terms of rms and mapped bias.The sector snowflake from the coarse search (square symbol) is highly ranked in all statistics except histogram fit.In fact in the average ranking the sector snowflake comes joint third out of the 26.The two fine-search combinations that beat it use six-bullet for large-scale snow and sector for convective snow.These two combinations have the same average rank and the only difference between them is the cloud ice, represented either by dendrite or sector snowflakes.This suggests that the choice of particle shape is less important for cloud ice than for precipitation.A triangle on Fig. 9 identifies the six-bullet/dendrite/sector combination.However, in terms of fit, there is only a marginal advantage over the sector experiment from the coarse search.
The final combination identified on Fig. 9 uses the sector snowflake for all three frozen hydrometeor categories (star symbol).This is best in terms of the skewness statistic but it is not particularly good in terms of rms and mapped bias.This again illustrates that where a single statistic is optimised, others will often degrade.This is further justification for basing the conclusions on more than one statistic.
For modelling convective snow, the attempt to increase the amount of scattering by using thin plate or block shapes was not successful.These shapes produce most of the worst-ranked experiments in Table 4.In contrast, some of the higher-ranked experiments represent convective snow using the dendrite snowflake, which is in general the least-scattering shape.This suggests that if scattering is lacking in the tropics, it is the large-scale snow category rather than the convective category that needs attention.Perhaps the scattering properties of large-scale snow need to be different between the tropics and the midlatitudes.In these attempts to further improve the modelling of scattering, it has been very hard to do better than the sector snowflake experiment from the coarse search.Almost all the available improvement over the Mie sphere has been gained by going to the DDA sector shape in the coarse search.Additional refinements bring very little further benefit; this helps justify the strategy outlined in the Introduction of looking for a simple scheme that can be well tuned, rather than getting lost in a complex approach that is hard to tune or validate objectively.Improvements over the "one-size-fits all" approach will be found eventually, but they will require substantial further work.

Sensitivity to assumptions and inaccuracies in radiative transfer
The main fixed assumption in this study has been the tropical  Marshall and Palmer (1948) size distribution had been used instead, the results would have been much worse, with all the DDA shapes producing far too much scattering.As an example, Fig. 10 shows simulated and observed TBs from the centre of Hurricane Irene.Applying Marshall-Palmer to the sector snowflake reduces TBs (i.e. it substantially increases the amount of scattering) compared to the Field et al. results, and brings the simulations further from the observations.The Field et al. distribution emphasises the very small sizes in the distribution, with a consequent reduction in the numbers of large particles, and hence a reduction in the amount of scattering compared to the Marshall-Palmer distribution.This shows that conclusions on the "best" particle shape are entirely dependent on the chosen size distribution.However, at least the chosen Field et al. size distribution produces more physically plausible results with DDA shapes than does the Marshall-Palmer.
We did not evaluate the midlatitude version of the Field et al. (2007) size distribution, but sensitivity tests with the sector snowflake showed that changing from the tropical to the midlatitude version reduces brightness temperatures.In convection and frontal systems, reductions are of order 5 K at 90 GHz and higher frequencies.The sector snowflake with the tropical size distribution already generates slightly too much scattering in the midlatitudes, so going to the midlatitude size distribution would make things worse unless we also changed the particle shape.
For the future, we could consider using ensembles of particle shapes.One advantage, demonstrated by Kulie et al. (2010), is the ability to blend together the optical properties of more than one shape.We have seen that the optical properties from available Liu (2008) shapes can fall either side of the best fit to observations (e.g.Fig. 5) and a blend of the two best shapes might give a better fit.A second advantage (e.g.Baran and Labonnote, 2007) is that the ensemble weighting need not be constant with size, and small particles could be represented more realistically by pristine hexagonal prisms and larger shapes by snowflakes.We might also have got better results had we included models for ice aggregates (e.g.Petty and Huang, 2010).Finally, we have only considered randomly oriented particle shapes, but preferential orientation might have some effect on the brightness temperatures, though at 183 GHz and below the effects may be obvious only in limited areas and are unlikely to be larger than around 5 K (e.g.Prigent et al., 2001;Davis et al., 2005).
In addition to the settings that directly affect the bulk scattering properties of frozen hydrometeors, there are many other uncertainties in cloud and precipitation radiative transfer.Uncertainties in particle shape and size distribution affect rain hydrometeors too.Further, Bennartz and Greenwald (2011), among others, have raised concerns about the accuracy of the solver for scattering radiative transfer and the accuracy of the plane parallel approximation, in other words the lack of sub-grid and 3-D cloud and precipitation structure.
The solver for scattering radiative transfer in RTTOV-SCATT is the delta-Eddington method.As a variant of the two-stream solution, the delta-Eddington method might appear crude in comparison to the many other solvers available (e.g.Thomas and Stamnes, 1999).Indeed four-stream rather than two-stream solvers have been recommended for improved accuracy in scattering calculations in the Community Radiative Transfer Model (Bennartz and Greenwald, 2011).However, the delta-Eddington has been shown to be accurate (e.g.Smith et al., 2002;Kim et al., 2004).A reverse Monte Carlo solver was experimentally implemented in RTTOV-SCATT to compare to the accuracy of the delta-Eddington, but it was found that differences between the two solvers were small in the NWP context.The use of the delta-Eddington solver is not an important source of error.
Sub-grid variability is represented in RTTOV-SCATT using the effective cloud fraction of Geer et al. (2009a, b).This is a computationally efficient, "first order" solution that represents the model grid-box brightness temperature as a weighted average of the completely clear and completely cloudy brightness temperature.If computational efficiency were irrelevant, it might be preferable to use the multiple independent column approach (ICA) to describe the effects of sub-grid cloud and precipitation variability.In this approach, the grid box is divided into multiple sub-columns, and the cloud and precipitation is distributed among those sub-columns according to cloud overlap rules.The DDA sector snowflake simulations of Sect.5.1 were repeated using the ICA approach with 20 sub-columns (20ICA).The cloud and precipitation overlap scheme of O'Dell et al. (2007) was used to fill the sub-columns.Without going into detail, results were not too different compared to the normal RTTOV-SCATT, and there were both degradations and improvements.As an example, the rms of FG departures was changed, for good or ill, by no more than 7 %.As originally shown by Geer et al. (2009a), the effective cloud fraction used in RTTOV-SCATT is a reasonable approximation to the ICA, at least in the context of data assimilation.
The problem with all plane-parallel radiative transfer, including the ICA, is that it does not represent the slanting path of the radiation through the atmosphere.For slant paths, emission from the sides of clouds can be as important as emission from the tops (e.g.Weinman and Davies, 1978;Roberti et al., 1994).A typical microwave imager zenith angle is 53 • , so microwave imagers are particularly susceptible to these effects.O'Dell et al. (2007) and Bennartz and Greenwald (2011) have found that slant path errors can be as much as 20 K in microwave imager channels in cases where (to simplify a little) the instrument's field of view is dominated by cloud sides rather than cloud tops.There are two obvious situations where this may occur: in maritime cumulus and convection.In convection, precipitation shafts may form only a small part of the horizontal domain, but viewed obliquely, these shafts become a much more important part of the radiative transfer.Further improvements in scattering radiative transfer in convective situations may just as likely come from representing 3-D issues as from further attention to the snow microphysics.
A final concern is the observation resolution.Brightness temperature histograms such as Fig. 5a are affected by the size of the instrument's field of view.Extreme values of TB such as those associated with convection are often very localised, so the larger the field of view, the less likely it is to observe an extreme TB.This study has made comparisons using superobs in 80 km squares to roughly match the effective resolution of the model's cloud and precipitation.
To check the sensitivity of the results to the resolution, the coarse-search experiments were repeated using the full native resolution of TMI and SSMIS (i.e.not using superobs) and a T1279 (roughly 16 km) model resolution for 4 days of Hurricane Irene during August 2011.Again, the sector snowflake was the best-performing experiment.This suggests that the results are robust.

Results over land
Results over land surfaces are presented separately because they are quite different from those over ocean.Figure 11 shows the statistics of fit for the coarse-search experiments over land.With land surface emissivities in the range 0.7 to 1, cloud and precipitation emission in the lowermost troposphere is less important and the greatest atmospheric signal comes in the higher frequencies from frozen hydrometeors.Over land, both the Mie sphere and the lower-scattering DDA shapes (e.g.sector snowflake, three-bullet, six-bullet) produce too little scattering, resulting in negative skewness (Fig. 11c) and poor values of the histogram fit (Fig. 11d).Instead it is the strongly scattering thin plate that gives good results by these measures.However, the thin plate causes large increases in rms errors in the higher frequencies (Fig. 11b) and in the mapped biases (Fig. 11e) because of excessive scattering in the midlatitudes (no figure shown). 1 0 v 1 0 h 1 9 v 1 9 h 2 2 v 3 7 v 3 7 h 5 0 .3 5 2 .8 5 3 .6 9 2 v 9 2 h 1 5 0 1 8 3 ± 7 1 8 3 ± 3 1 8 3 ± 1 1 0 v 1 0 h 1 9 v 1 9 h 2 2 v 3 7 v 3 7 h 5 0 .3 5 2 .8 5 3 .6 9 2 v 9 2 h 1 5 0 1 8 3 ± 7 1 8 3 ± 3 1 8 3 ± 1 Channel [GHz] Channel [GHz] Mie  The effective cloud overlap of Geer et al. (2009a, "C av ") has not previously been tested over land surfaces, so two additional experiments were run with the sector snowflake.One, the same as presented in the over-ocean results, used the 20ICA approximation.The other one used the old "C max " cloud overlap over land surfaces.This takes the largest cloud fraction in the model profile to represent the effective cloud fraction.The results from C av and the 20ICA were quite consistent (not shown) indicating that the Geer et al. cloud overlap is a reasonable approximation to 20ICA results over land, just as it is over ocean.Surprisingly, however, the C max cloud overlap produced good results over land with the sector snowflake and in fact it was the best performer in the coarse search.This is clear to see in Fig. 11.Using C max always produces higher effective cloud fractions than the C av approach.A higher effective cloud fraction means an increase in the weight given to the cloudy column in RTTOV-SCATT and thus a greater influence of scattering on the simulated brightness temperature.Using C max helps compensate for the lack of scattering produced by the sector snowflake.The success of the sector snowflake and C max cloud overlap has practical benefits, even if the physical realism of the C max approximation is questionable.It means it is possible to get good results across land and ocean surfaces using the sector snowflake for all snow, as long as the cloud overlap is varied from C av to C max according to whether the surface is ocean or land.This is straightforward to implement technically in the ECMWF system, and it will likely be adopted for future operational implementation (Baordo et al., 2013) The fine search over land was carried out using C av because it is the truest representation of the 20ICA results (even if Sect.5.3 questions whether 20ICA is the ideal reference, given that 3-D effects may be important).The fine search is summarised in Fig. 12 but it is not examined in much detail apart from mentioning that it indicates that the real issue over land is the treatment of snow in convection.Use of thin-plate or block column shapes for convective snow produces substantial improvements in skewness and histogram fit without degrading the mapped bias.All the best experiments used thin plate or block column for convective snow.The best combination was six-bullet (largescale snow)/dendrite (cloud ice)/block (convective snow), indicated by a star.The only difference from the winner over ocean is the use of block rather than sector snowflake for convective snow.This suggests that the large-scale snow and cloud ice can be modelled in the same way over ocean and land but convective snow requires very different treatment.
One possible explanation for the discrepancy between land and ocean surfaces is forecast model bias.Figure 13 shows brightness temperature histograms for the 183 ± 1 GHz channel on SSMIS, which is sensitive to upper tropospheric moisture (e.g.Buehler and John, 2005) and to deep convection (e.g.Hong et al., 2005).It should not be affected by issues with the land surface emissivity.Low brightness temperatures, which will be taken to signify deep convection, are around twice as common over land as over ocean (e.g.TB < 235 K in 0.6 % of observations over ocean compared to 1.2 % over land).The simulations show the opposite and this is true even in the forecast model itself (for example, integrated snow water path > 2 kg m −2 in 1.6 % of simulations over ocean compared to 1.2 % over land).However, it is still possible that there are errors in the radiative transfer model that only manifest themselves over land surfaces, or maybe there are real physical differences between land and ocean, such as in the microphysics of snow and graupel.Further work is needed.

Conclusions
Simulating the bulk optical properties of snow hydrometeors using Mie spheres and the Marshall-Palmer size distribution leads to unphysically high amounts of scattering in middle frequencies (30-50 GHz) and too little scattering at high frequencies .Changing the density model or the size distribution can improve results at some frequencies but it is hard to avoid degrading the results at other frequencies.Using discrete dipole results is a better choice but the problem remains as to how to choose an appropriate particle shape (or shapes) and a size distribution.
The ECMWF data assimilation system provides a framework in which modelled clouds and precipitation can be used to drive a radiative transfer model and hence to compare simulated brightness temperatures to their observed equivalents.The statistics of the departures between model and observations can be used to optimise the choice of size distribution and particle shape for snow.The chosen statistics were the rms and skewness of the departures, a statistic to quantify the discrepancies between simulated and observed histograms of brightness temperatures, and the rms of mapped departure biases.The latter three statistics are resistant to the double penalty problem, and help to indicate whether changes in the rms have occurred for good or bad reasons.
Because of the lack of predictability of clouds and precipitation at the smaller scales in NWP models, the error budget of any comparison between model and observations is often dominated by the imperfect shape, size and intensity of modelled cloud and precipitation features, leading to rms errors of 20-40 K in brightness temperature terms.This makes it hard to objectively justify the less radiatively significant changes to a radiative transfer model, but equally it allows the use of relatively imprecise radiative transfer modelling in the observation operator.Nevertheless, the errors from using Mie spheres to model snow hydrometeors were obvious in this context, and they appear to have been the largest remaining source of error in the ECMWF all-sky assimilation of microwave observations.
Compared to using the Marshall-Palmer distribution and the Mie sphere to represent snow particles, the tropical version of the Field et al. (2007) size distribution with the Liu (2008) sector snowflake can reduce rms errors by 40 % in the 50 GHz channels and by smaller amounts in other channels.Simulated brightness temperatures are improved by up to 70 K in deep convective situations.The over-and underprediction of scattering at different frequencies, characteristic of the Mie spheres, has mostly been removed.Any Liu shape with similar bulk optical properties to the sector snowflake (e.g.three-and six-bullet rosettes) produces fairly similar results.
The story is somewhat different over land surfaces, possibly due to different bias characteristics in the forecast model over ocean and over land.The forecast model appears to produce less frequent deep convection over land as over ocean, contradicting what is seen in the observations, which have roughly twice as much convection over land than over ocean (note that this applies to the sample of observations in this study and is not intended as a general statement).Good results over land are again achieved with a sector snowflake, but it is necessary to boost the amount of scattering by using a larger effective cloud fraction (precisely, by using C max rather than C av ).This helps increase the amount of scattering affecting the simulated brightness temperatures.
Attempts were made to further improve the fit between model and observations by using DDA shapes for cloud ice and by splitting the snow hydrometeor type into a convective and a large-scale category.However, there was little additional benefit compared to the simple approach of using the DDA sector snowflake to represent all snow.More complex approaches will surely bring improvements in the future, but it has been hard to do this so far.The "one-size-fitsall" sector-snowflake model is straightforward to implement and it provides major benefits over the Mie sphere, so it will be adopted at ECMWF and will be the default in RTTOV-SCATT.
Practically, the discrete-dipole shapes produce better results not because of some broad, tunable change in scattering properties at all frequencies, but because scattering declines more rapidly at the lower frequencies compared to the Mie sphere.This advantage appears to come from a combination of two main differences between Mie sphere and nonsphere bulk optical properties (Figs. 2 and 3): (a) the faster decline in extinction and SSA towards lower frequencies and (b) the fact that the discrete dipole shapes produce much less forward-peaked scattering than do the Mie spheres at high snow water contents.
If we consider the work of improving cloud and precipitation radiative transfer as a multi-dimensional optimisation problem, we have explored one dimension here: that of the particle shape.Other "dimensions" have been held fixed simply to make the problem tractable: the particle size distribution, the possibility that the amount of cloud and precipitation generated by the forecast model is incorrect, the treatment of sub-grid variability and 3-D radiative transfer.Further, it may be hard to generalise because the hydrometeor categories used in the ECMWF model are different to those in other models or retrieval schemes.Hence it is fair to question whether the Field et al. (2007) size distribution with the sector snowflake will always be the best "one-size-fits-all" method of representing snow.Nevertheless, among the sensitivity tests that were carried out over the ocean, no better configuration could be found.The results were reasonably good across different weather situations in the tropics and midlatitudes and across all microwave channels from 10 GHz to 183 GHz, which is a major improvement over the Mie sphere.The "multi-dimensional search" is probably going in the right direction.
To make further improvements in scattering radiative transfer in cloud and precipitation, it is clearly necessary to continue optimising the choice of size distribution and particle shapes, and to look again at properly representing the different optical properties cloud ice, convective and large-scale snow.But that is hard, as shown in this study.3-D radiative transfer effects can also affect brightness temperatures by 10 to 20 K. Dealing with the 3-D radiative transfer issue may also be important for future improvements, especially in convection.
Finally, the new optical properties have improved snowscattering radiative transfer to the point that it has become possible at ECMWF to start operational assimilation of 183 GHz water vapour sounding channels in all-sky conditions.This brings real forecast benefits (Geer, 2013).the different hexagonal prism shapes.The mistake has been inherited by the current work but it only affects the hexagonal prisms and it does not affect our conclusions on the qualities of those shapes.
In the case of the Liu (2008) rosettes it is also possible to derive m(D) analytically, but the mathematics are inconvenient and they do not lead to an exact power law solution.Also, the description of the construction of these shapes is not fully enough specified (e.g.whether the bullets intersect in the centre point or are just stuck together).The sector snowflake does have a power law solution but there is no practical way of analytically computing a mass-density relationship for the complex dendrite snowflake shape.
Practically, a better way to derive the coefficients a and b is to make a functional fit to the mass of the actual shapes used in the DDA computations.The mass is reported in the Liu database as a function of particle size.However, the way these shapes have been implemented leads to small steps and jumps in the mass as a function of dimension; in practice the mass is not an exact power law function of dimension.
As an example, Fig. B1 shows the mass of the sector snowflake, both as reported from the Liu database and as modelled with a variety of power law relationships: the Kulie et al. (2010) fit, an alternative fit by the current authors, and the analytical solution based on the description of this shape in Liu (2008).All of these are quite similar, but all fail to fit the actual values from the Liu database at very small particle dimensions where it diverges from a simple power law relationship.However, these particle sizes are the least important to the radiative transfer.Also there is an example of a "step" feature in the Liu database at around 400 microns, which also cannot be exactly represented by the power law.However, the practical effect on scattering properties is small.

Figure 1 .
Figure 1.Observed and simulated over-ocean brightness temperatures on 25 August 2011 in the region of Hurricane Irene.Simulations are generated from ECMWF forecasts at T1279 horizontal resolution with snow hydrometeors represented either by Mie spheres or by DDA sector snowflakes.Channels at 10 GHz (h = horizontal polarisation) and 37 GHz (v = vertical polarisation) come from an overpass of TMI at 21:18 UTC.Channels at 52.8 GHz and 150 GHz (both h polarisation) come from an SSMIS overpass at 22:12 UTC.The domain runs from 84 • W to 70 • W and from 22 • N to 34 • N; grid spacing is 2 • .Land-contaminated 10 GHz observations on the Florida coast are not used in quantitative comparisons.

Figure 2 .
Figure 2. Bulk optical properties at 52.8 GHz and 253 K, as a function of snow water content.Optical properties of the Liu (2008) shapes have been computed with the Field et al. (2007) size distribution.SomeLiu (2008) shapes have similar bulk scattering properties to others and are ignored: the five-bullet rosette is similar to the six-bullet; the short column is similar to the thick plate.

Figure 4 .
Figure 4. Histograms of first guess departures [K] for the 52.8 GHz channel on SSMIS using different scattering models.Bin size is 2.0 K.The sample is composed of all superobs in June 2012 passing the basic quality control tests described in Sect.2.2.

Figure 5 .
Figure 5. How the measure of histogram fit is constructed, using the 52.8 GHz channel on SSMIS as an example.(a) Histograms of observed and simulated brightness temperature [K] using different scattering models.Bin size is 2.5 K; Empty bins have been filled with the value 0.1.(b) Log of the ratio of histograms (simulation divided by observation).

Figure 6 .
Figure 6.June 2012 mean FG departures in the 52.8 GHz channel on SSMIS for three choices of snow optical properties.Means are computed in bins of 5 • by 5 • in latitude and longitude.

Figure 7 .
Figure 7. Measures of fit between model and observations over ocean surfaces: (a) The rms of FG departures, normalised by the Mie sphere results; (b) skewness of FG departures; (c) brightness temperature histogram fit; (d) rms of mapped biases.Lower numbers indicate better agreement between model and observations; in the case of the skewness, large negative or positive numbers are bad.

Figure 9 .
Figure9.Measures of fit in the fine search, averaged across all channels, ranked from best (1) to worst (26) and then normalised so that the worst equals one.The positions of a few important experiments are identified by symbols.

Figure 10 .
Figure 10.Influence of the size distribution on the Liu (2008) sector snowflake results, for a single case in the centre of Hurricane Irene.The thin solid line shows results with the Field et al. (2007) size distribution; the thin dashed line shows the Marshall-Palmer results.The thick dashed line shows the observations.

Figure 11 .
Figure 11.Measures of fit between model and observations over land surfaces: (a) The rms of FG departures, normalised by the Mie sphere results; (b) skewness of FG departures; (c) brightness temperature histogram fit; (d) rms of mapped biases.

Figure 12 .
Figure12.Measures of fit in the fine search over land, each ranked from best to worst and normalised so that the worst equals 1.

Figure 13 .
Figure 13.Normalised histograms of 183 ± 1 GHz brightness temperatures in June 2012 over land and ocean.Model results are generated using the sector snowflake experiment from the coarse search.

Figure B1 .
Figure B1.Particle mass as a function of size for the Liu sector snowflake.The x axis scale is in microns for convenience but the fits reported on the legend are based on SI units.

Table 2 .
Measures of fit between model and observations in the coarse search over ocean, summarised across all channels.To summarise the rms, histogram fit and rms of mapped biases, a mean is computed across all channels.For the skewness, an rms is computed across the channels.In all cases lower numbers indicate better agreement between model and observations.

Table 3 .
Frozen hydrometeor shapes used in the fine search, ranked roughly from lowest to highest scattering.

Table 4 .
Fine-search combinations ranked according to their average position in the statistics of fit.