Physical characteristics of frozen hydrometeors inferred with parameter estimation

Frozen hydrometeors are found in a huge range of shapes and sizes, with variability on much smaller scales than those of typical models or satellite observations. Neither models nor in-situ measurements can fully describe this variability, so assumptions have to be made in applications including atmospheric modelling and radiative transfer. Here parameter estimation has been used to optimise six different assumptions relevant to frozen hydrometeors in passive microwave radiative transfer. This covers cloud overlap, convective water content and particle size distribution (PSD), the shapes of large-scale snow and 5 convective snow and an initial exploration of the ice cloud representation (particle shape and PSD combined). These parameters were simultaneously adjusted to find the best fit between simulations from the European Centre for Medium-range Weather Forecasts (ECMWF) assimilation system, and near-global microwave observations covering the frequency range 19 GHz to 190 GHz. The choices for the cloud overlap and the convective particle shape were particularly well constrained (or identifiable) and there was even constraint on the cloud ice PSD. The practical output is a set of improved assumptions to be used in version 10 13 of the Radiative Transfer for TOVS microwave scattering package (RTTOV-SCATT), taking into account newly available particle shapes such as aggregates and hail, and additional PSD options. The parameter estimation explored the full parameter space using an efficient assumption of linearly additive perturbations. This helped illustrate issues such as multiple minima in the cost function, and non-Gaussian errors, that would make it hard to implement the same approach in a standard data assimilation system for weather forecasting. Nevertheless, as modelling systems grow more complex, parameter estimation is 15 likely to be a necessary part of the development process.


Introduction
Clouds and precipitation are some of the most uncertain processes in the earth system, leading to systematic errors in models (e.g. Klein et al., 2009;Forbes et al., 2016) and big uncertainties in climate change predictions (e.g. Zelinka et al., 2020). 20 The parametrisations that represent hydrometeors 1 in global models rely on physically-informed heuristic abstractions, such as the representation of convection by upward and downward mass fluxes (e.g. Tiedtke, 1989). They also rely on compressing 1 Here, the words hydrometeor and particle will be used interchangeably to describe liquid and frozen water particles in the atmosphere by the updraft and downdraft mass fluxes, which are assumed to occupy a fixed 5% of the grid box. Cloud water and ice can be detrained into the large-scale scheme, allowing convection to form anvil cloud.

Observation operator
The observation operator for satellite radiances in the IFS is RTTOV (Saunders et al., 2018). Within this, RTTOV-SCATT (Bauer et al., 2006) provides the all-sky microwave capabilities. The all-sky radiance simulated by RTTOV-SCATT is the 95 weighted combination of radiances from two independent column calculations, one that is completely clear, I clear and one with horizontally homogeneous cloud and precipitation, I cloud . The cloudy column is weighted by an effective cloud fraction C to approximate the nonlinear effect of sub-grid variability of cloud and precipitation, which is sometimes known as the beamfilling effect (e.g. Geer et al., 2009a;. Hence the all-sky radiance is: This is then converted to black-body equivalent brightness temperature (TB, in K) which is the usual, more intuitive variable to represent satellite radiance observations. The effective cloud fraction C simplifies what could be a complex 3D arrangement of cloud and precipitation. Two options for deriving C will be tested in the parameter search (Sect. 3.1).
In Eq. (1) the clear column accounts only for the surface interaction and gaseous absorption. In the cloudy column, scattering radiative transfer is represented using the delta-Eddington approach (Joseph et al., 1976;Kummerow, 1993;Bauer et al., 2006). 105 Although simple compared to more exact scattering methods, this is fast enough for use in a weather forecasting system and still accurate compared to reference models.
One simplification in the scattering radiative transfer is the treatment of polarisation. Each polarisation is considered independent and hydrometeor optical properties are invariant with polarisation. This neglects the polarising effects of scattering from particles, which can transfer energy from one polarisation to another, and causes optical properties to vary as a function 110 of polarisation and viewing angle. RTTOV v13 will have a simplified treatment of polarised scattering from preferentially oriented frozen hydrometeors; this reduces errors in the simulated vertical-horizontal polarisation difference by 10 -15 K ( . This is not used in the main study; however the parameter estimation is updated at the end, taking account of the new polarisation scheme, to provide a self-consistent final configuration for RTTOV v13 (Sect. 5.4).
The surface interaction is treated as specular reflection. Over ice-free ocean surfaces, this is handled by the FASTEM-6 115 surface emissivity model (Kazumori et al., 2016), which includes a simplified treatment for diffuse scattering. Over sea-ice and over land, an all-sky dynamic emissivity retrieval operates at lower frequencies to supply the surface emissivity for the 50 GHz and 183 GHz sounding channels (Baordo and Geer, 2016). The TELSEM surface emissivity database (Aires et al., 2011) is used for the other channels of SSMIS over land surfaces and also as a backup when the retrieval is not possible. In the passive monitoring experiments used here, a change to physical assumptions in the radiative transfer can affect the dynamic emissivity retrieval. In practice, the effect is negligible and even in the convective graupel experiment (the biggest change tested, see Table 2 later) the maximum change in surface emissivity is tiny at around 1e −6 .
The scattering solver requires the bulk optical properties of each layer of atmosphere. Hydrometeor optical properties (extinction, single-scattering albedo, and the asymmetry parameter) are supplied from lookup tables as a function of frequency, temperature, hydrometeor type and hydrometeor water content. Changes to the microphyscial assumptions explored in this 125 work were made by changing these lookup tables, which are generated by an offline tool available within RTTOV . This integrates single-particle optical properties over the assumed particle size distribution, given the water content. It offers a choice of hydrometeor representations based on Mie spheres or non-spherical particles from the Liu (2008) and ARTS  databases. Version 13 includes new options for particle size distributions that were added in support of the current work. From the available options, four sets of parameters relating to large-scale snow, convective snow and cloud 130 ice were added to the parameter estimation as will be described in Sect. 3.
As discussed, the IFS has 4 large-scale prognostic hydrometeors, a prognostic cloud fraction, and a diagnostic large-scale precipitation fraction. Additionally, diagnostic mass fluxes for snow and rain represent convection, which is assumed to take 5% of the gridbox area. This is the complete set of cloud and precipitation information provided to RTTOV-SCATT, although some conversions are needed. To standardise the treatment with the convection fluxes, all six hydrometeors (rain, snow, cloud 135 ice, cloud water from the large-scale scheme; rain and snow from the convection scheme) are extracted from the model as fluxes (kg m −2 s −1 ) and then converted to mass mixing ratios (kg kg −1 ) before input to RTTOV-SCATT (see Geer et al., 2007, appendix B). Strictly, this ratio is defined as the mass of the hydrometeor with respect the mass of moist air. However, due to the lack of a widely-used terminology to specify whether hydrometeor mixing ratios are based on dry or moist air, this quantity is referred to throughout as just "mixing ratio". The assumptions involved in obtaining the mixing ratio from the flux make the 140 convective mixing ratio particularly uncertain (Geer et al., 2017a), but rather than adjust the assumptions in the model, here the parameter estimation is used to estimate a scaling factor to adjust the convective snow mixing ratio directly -further details in Sect. 3.
In the IFS, RTTOV-SCATT has been used in a four-hydrometeor configuration until now (rain, snow, cloud water, cloud ice) with the convective rain and snow mixing ratios added together with the large-scale rain and snow respectively; this provides 145 the control setup. Going forward, the aim is a five-hydrometeor configuration, adding a new hydrometeor type specifically for convective snow, which is no longer combined with the large-scale snow. Within RTTOV-SCATT this hydrometeor is loosely referred to as 'graupel' but strictly it is for convective snow, making no assumptions on the particle habit. The most important causes of uncertainty in the radiative transfer are these microphysical assumptions, along with the representation of sub-grid heterogeneity (Bennartz and Greenwald, 2011;) -hence the attempt to improve these aspects 150 through parameter estimation.

Observations
This study uses observations from the Special Sensor Microwave Imager / Sounder (SSMIS, Kunkee et al., 2008). This is a conical-scanning microwave sensor, meaning that it has an approximately fixed zenith angle and polarisation across its swath. There are 24 channels on the instrument of which 13 have been selected for use in this study (Table 1). Although SSMIS has 155 additional calibration issues (e.g. Bell et al., 2008) compared to later microwave imagers such as the GPM microwave imager (GMI, Draper et al., 2015) these are mostly corrected by bias correction and it has been extensively used at ECMWF for all-sky assimilation and parameter estimation (Geer et al., 2017b;Geer and Baordo, 2014). The advantage of using SSMIS is that, with its temperature sounding channels around 50 GHz, it is a single instrument that covers the majority of frequencies currently in use or planned for use in all-sky assimilation , with the main exception being frequencies above 190 GHz 160 which will not be available from space until the launch of ICI in 2024. For the IFS, the observations are superobbed onto an 80 km by 80 km grid (Geer and Bauer, 2010) which averages together at least 9 original SSMIS observations. This standardises the resolution across frequencies and helps match the effective resolution of cloud in the forecast model (note the simulations are made using the model grid point nearest to the centre of the superob).
For the parameter estimation, SSMIS observations are taken from the F-17 satellite over a 10 day study period from 13 -165 22 June 2019 inclusive. The selection of observations aims to include as much data as possible, even in situations where the observations are not well-enough represented to be used in active assimilation. The philosophy is that unless difficult situations are included as a constraint in the parameter estimation, the quality of radiative transfer in those situations will never become good enough to allow them to be used. This accepts the risk that unrelated systematic errors, such as poor surface emissivity, will lead to compensating biases in the parameter estimation.  Figure 1 illustrates the number of observations available in the study period, using a logarithmic colour scale. These are accumulated in the same 10 • by 10 • latitude-longitude bins that will be used in the cost function for the parameter estimation (Sect. 4). Observations are only excluded in areas of orography higher than 800m, areas of sea-ice, in mixed scenes where the grid point land-water mask is neither completely land or completely ocean, and over land where both the dynamic emissivity retrieval and the atlas backup failed. In any one bin, observations must be all ocean or all land, favouring the larger number.

175
Typically around 3000 observations are available per bin over ice-free oceans and around 1000 observations over land areas; an approximate total of 1 million observations is used. The main losses are due to sea-ice, the orography check and the restriction on mixed scenes.

Microphysical and macrophysical options
This section will outline the 6 dimensions of the parameter search, summarised in Table 2. This is a search among discrete 180 possibilities, so the total number of possible configurations is the product of the dimension sizes: 2 × 3 × 2 × 6 × 6 × 8 = 3456.
Along each dimension, Table 2 orders the options according to the all-channel global mean change in simulated brightness temperature when that parameter is varied. The main impact of frozen hydrometeors at microwave frequencies is to scatter radiation, which generally leads to lower brightness temperatures at the top of the atmosphere. Hence the dimensions are ordered from the least scattering options (which increase brightness temperatures compared to the control) to the most scattering 185 (which decrease them). The biggest increase in scattering of any of the options is the move to a convective graupel particle shape, which reduces brightness temperatures globally by 0.187 K. Such global mean changes are obviously small, but this is because cloud and precipitation are localised processes. In heavy precipitation there are still changes of up to around 100 K, as will be illustrated later; see also e.g. Geer et al. (2021) for further illustration of the very strong sensitivity of simulated brightness temperatures to microphysical assumptions.

190
Although it would be tempting to add more dimensions to the parameter search, it soon becomes cumbersome and would ultimately become impossible due to the curse of dimensionality. The range of options already encompasses research that will  Table 3 8 take many pages to summarise. More dimensions could have been included, particularly in the forecast model physics, but that is out of scope of the current study.
3.1 Cloud overlap over land

195
As described in Sect. 2.2, sub-grid variability of cloud and precipitation has a huge effect on the simulated brightness temperatures. This is represented by the effective cloud fraction C in the two-column approach used by RTTOV-SCATT (Eq. (1)).
Over ocean, C is computed as a hydrometeor-weighted vertical average of the hydrometeor fractions (C av ) which is a reason-able approximation to more accurate multiple-independent column models (Geer et al., 2009a, b); this remains fixed in the parameter estimation. But over land surfaces, C has previously been set to the maximum cloud fraction in the profile, C max , not 200 because this is the most physically correct approach but to counter systematic errors: Geer and Baordo (2014) took advantage of the fact that C max generates artificially high scattering TB depressions to compensate the apparent lack of convective cloud and precipitation over land in the IFS, as seen at frequencies from the microwave to the infrared (e.g. Geer et al., 2019). The hope is that the parameter search can find a configuration that does not require a physically unjustified representation of cloud overlap over land.

205
The choice of effective cloud fraction over land affects mainly deep convective situations. In the IFS, convective locations are represented by a convective core occupying 5% of the model gridbox; the anvil is represented by detrainment into the large-scale cloud scheme, which typically results in a thin capping ice cloud with a cloud fraction that is close to 1. In the computation of C av , the convective core typically dominates the mass-weighted average and hence C av can be as low as 0.05 in convection. By contrast the maximum cloud fraction C max is typically found in the anvil cloud and is typically closer to 1.

210
Hence, using C max in Eq.
(1) gives more weighting to the cloudy column and generates much deeper TB depressions. Figure 2b shows that going to C av increases the mean simulated brightness temperature in many land-surface areas. These locations are similar to those affected by increasing convective snow amount ( Fig. 2 c,d; see next section) confirming that it is mainly convective situations affected.

215
The physical parameters of convective snow are not well-known and may explain many of the largest errors in the observation operator. These uncertainties are explored by dimensions 2, 3 and 4 of Table 2. As mentioned earlier, the convective snow mixing ratio is treated as uncertain and allowed to vary using a scaling factor; this is dimension 2. As discussed in the previous section, the observational evidence from all-sky infrared assimilation (Geer et al., 2019) would be consistent with the convective mixing ratio being too low, particularly over land. However there are also reasons to suspect it may be too high. As described in 220 Sect. 2.2 the convection parametrisation represents snow as a bulk mass flux and additional uncertainty arises in the conversion to mixing ratio, which makes assumptions on the PSD and the fall speed distribution. The assumed fall speeds are suspected of being too low for convective particles (Geer et al., 2017a, section 3.3). This would lead to convective snow mixing ratios being overestimated by about 50%. Hence the second dimension has options both to increase or decrease the mixing ratio by this amount. Figure 2c and d show the impact on mean brightness temperatures and reveal the locations where significant 225 convection occurred during the 10 day study period: around the ITCZ, in the SH storm-tracks, and at midlatitudes over land surfaces.
A second possible change is the size distribution of convective particles. The control uses the Field et al. (2007, F07) tropical (T) PSD for all frozen precipitation. However, this was derived from in-situ measurements of tropical anvils and midlatitude stratiform clouds so it may not be appropriate for the large frozen particles within the convective cores, which are not accessible 230 to in-situ measurement from aircraft. Indirect evidence from hail reaching the ground and from remote sensing suggests that the most intense convection can generate particles up to at least the size of golf balls (e.g. Zipser et al., 2006). For more typical  cases, there is still uncertainty as to the amount of large (e.g. cm-sized) particles in the convective cores (e.g. Elsaesser et al., 2017). Hence the Marshall and Palmer (1948) PSD is provided as an alternative in the parameter search; as implemented here this boosts the number of cm-size particles by at least 2 orders of magnitude compared to the F07 PSD, and hence gives more  scattering (see Geer et al., 2021, for more details). The effects of this are not illustrated as they are broadly similar to increasing the convective snow mixing ratio by 50%.
The final potential adjustment for convective snow is the particle shape. Note that the particle shape also implies a mass-size relation and hence this affects the bulk scattering properties through adjustments both to the size distribution and the singlescattering characteristics (Eriksson et al., 2015;Geer et al., 2021). The Liu (2008) sector snowflake is used in the control, 240 since this was the best particle to represent the sum of convective and large-scale snow in the parameter search of Geer and Baordo (2014). However it is an idealised smooth and solid particle that will be inadequate to represent scattering in future sub-mm applications (Fox et al., 2019). The ARTS database offers new options of ice aggregates and graupel and these may be more physically representative . Hence in the fourth dimension of the search, particle types have been chosen to sample different levels of scattering and also some very different physical representations: the ARTS column 245 aggregate 4 generates less scattering than the Liu sector, whereas the ARTS block aggregate, column and gem graupel generate increasingly more scattering. The options also include the Liu (2008) 3-bullet rosette, which was found by Geer and Baordo (2014) to give the best results when convective snow was represented separately to large-scale snow. It is striking that the particle shape assumptions illustrated in Fig. 2e and f have a larger impact on the simulated TB than a 50% change in the mass mixing ratio (as seen in Fig. 2c and d).
250 Figure 3 illustrates the frequency-resolved global mean change in brightness temperature coming from changes in convective particle shape. Graupel generates significantly increased scattering down to 19 GHz and has its biggest effect around 92 GHz.
The other chosen particles do not change the results so much at frequencies below 92 GHz, but there is plenty of variability in the higher frequencies. For example, the ARTS column and ARTS block aggregate would reduce TBs by similar amounts at 92 GHz, but they have different effects at 150 GHz. Hence there is a spectral signature associated with the snow particle shape 255 that should help resolve any ambiguity with changes in mixing ratio. For further illustration of the spectral signatures of these particles, see Geer et al. (2021); in heavy cloud situations these choices can change the simulated brightness temperature by up to 150 K.

Options for large-scale snow
The representation of large-scale snow is not perfect but probably better than that of convection. Hence only a single dimension 260 is explored, that of particle shape, as shown in Table 2. The Field et al. (2007) tropical PSD was retained for large-scale snow in all experiments, mainly to reduce the search space, but also because this PSD is well-established in this context (e.g. Fox et al., 2019). Among the particle shape options, the ARTS sector snowflake is a separate option to the Liu sector snowflake from the control, and it generates warmer mean brightness temperatures (see Geer et al., 2021, this is due to differing mass-size relations; the single-particle optical properties are near identical). The ARTS column and ARTS block aggregate generate more 265 scattering and lead to colder brightness temperatures. Figure 2g and h show the geographical effect of changing the large scale snow particle shape. This highlights areas of significant large-scale precipitation during the 10-day study period, including the SH storm-tracks, the W Pacific and the Atlantic coast of N America. Although there is some overlap with the spatial pattern of convection seen in Fig. 2c and d, there are also distinct differences -for example in N America, convection dominates over the midwest and large-scale precipitation over 270 the Atlantic coast. This suggests that the convective and large-scale physical parameters should be independently identifiable through their spatial signatures. And as with the convective snow particle choices, there is also spectral differentiation between the large-scale snow shapes that should further improve identifiability (see Geer et al., 2021).

Options for cloud ice
The current representation of ice cloud in RTTOV-SCATT uses a Mie sphere combined with a gamma PSD and is physically 275 unrealistic (Geer and Baordo, 2014). Ice cloud is expected to have a small effect on radiances at 183 GHz (e.g. Hong et al., 2005;Doherty et al., 2007) but in the control configuration it has almost no effect (not shown). However there is sensitivity at 183 GHz and below  so it is hoped to find a baseline representation of ice cloud that can be further improved once ICI has been launched. Table 3 summarises the ice cloud options in the parameter search and further details can be found in Geer et al. (2021). The 280 highest-scattering option is the combination of the Field et al. (2007) PSD and Liu (2008) dendrite, similar to the configuration which gave the best results for ice cloud in the fine search of Geer and Baordo (2014), but using the midlatitude option (M), rather than tropical option, to reduce the number of large particles. Figure 2j shows this configuration decreases brightness temperatures in many areas compared to the control. However, preliminary experimentation suggested there was too much scattering in this configuration. It is only with the new particle shapes and PSDs in RTTOV-SCATT that lower-scattering  options have become available . From the ARTS particle database, the ARTS column aggregate and the Evans aggregate generate less scattering. Also a range of alternative PSDs help reduce the number of large particles and hence reduce the overall amount of scattering. As shown in the table, progressively less scattering can be generated using the McFarquhar and Heymsfield (1997, MH97) and Heymsfield et al. (2013 PSDs. Early attempts at the parameter search suggested that even less-scattering options might be required, so two ad-hoc PSDs were created using the all-purpose modified 290 gamma distribution (Petty and Huang, 2011, MGD); these are known as the MGD 2e4 and MGD 1e4 PSDs, with the number referring to the lambda parameter in the MGD. The configurations using these PSDs make only small changes compared to the control (e.g. Fig. 2i). Figure 2j shows that ice cloud effects can be seen in convective and large-scale areas alike, so ice cloud has its own distinct spatial signature; there are also some frequency signatures (such as differences between the MH97 and H2013 options, thought 295 to be due to different temperature dependences in the PSDs, not shown). Hence there is hope that ice cloud parameters will be identifiable in this work.

Search
There are 3456 different parameter combinations available and it was not feasible to evaluate every combination directly, since 300 each would require a passive monitoring run of the IFS. Instead the cost function was evaluated approximately at each point in the search space. This was done using an assumption of linearly additive perturbations, which will be justified a-posteriori.
This meant that only one experiment needed to be run for each of the parameter options described in Sect. 3, in addition to the control -a total of 22 experiments. An alternative to running the full IFS would have been to archive the necessary geophysical inputs to RTTOV-SCATT for around 10 6 observations, but significant storage and computation resources would still have been 305 needed.
The assumption of linearly additive perturbations is expressed mathematically as: Here, T is the simulated brightness temperature, i the index across the approx. 10 6 observations and j the index across the 13 selected channels of SSMIS. k is the index across 6 search dimensions and v k the index within each search dimension. v is a 310 vector containing the six v k indices. T control is the brightness temperature simulated in the control experiment. As an example, In the parameter search, the effect of these choices on the brightness temperature was estimated as follows, using Eq. (2):

315
Note that for simplicity the i and j indices have been dropped here. Figure 4 illustrates the validity of the linear assumption for the "intermediate" configuration. This is a scatter plot based on the left and right sides of Eq. (4) after subtracting T Control . Although the estimate falls up to 20 K below the 1:1 line for the largest perturbations, the estimate has errors less than 2 K for smaller perturbations. The standard deviation of the difference between the estimated and actual perturbation is 0.21 K, suggesting that in this case the assumption is valid to reasonable 320 accuracy in most situations. Similar levels of accuracy have been seen with other configurations.

Cost function
In order to choose a unique "best" configuration, there must be a single objective way to measure the fit between the model and the observations. Following the data assimilation approach, this is provided by a cost function based on the departure d between the observation and the simulated equivalent: Here b i,j is an observation bias correction estimated using Variational Bias Correction (VarBC, Dee, 2004;Auligné et al., 2007) within the weather forecasting framework, which will be discussed shortly. In a data assimilation cost function, the equivalent "observation term" is usually based on the square of the departure normalised by the observation error, following an assumption of Gaussian errors (see e.g. Geer, 2021). The best fit to observations, known as the analysis, is obtained where the cost function 330 has a minimum.
The square of the departure is a poor choice of cost metric when it comes to cloud and precipitation. One reason is that modelled and observed clouds are not always in the same place, so it is often possible to reduce the squared error by removing the simulated cloud or precipitation -this is known as the double penalty effect. As explored by Geer and Baordo (2014), more appropriate cost functions for cloud and precipitation parameter estimation could include the mean and skewness of 335 departures, or measures of fit between the observed and simulated PDF of brightness temperatures, using measures similar to the Kullback-Leibler divergence (Kullback and Leibler, 1951). The PDF divergence approach was not pursued here because it is hard to combine with geographical binning (see below) which is used to help make the parameters more identifiable. A vastly larger dataset would be needed to ensure good sampling of the brightness temperature PDFs in each geographical bin.
In order to have a single cost function to minimise, a combination of the mean and skewness of the departures was chosen 340 for the current parameter estimation. Section 5.2 explores alternative choices. The mean and skewness are computed in 10 • by 10 • latitude-longitude bins (Fig. 1). Then the absolute value in each geographical bin is combined, giving the following cost function: Bin l (·) represents the geographical subset of observations in bin l, Mean(·) computes the mean of a set of departures and Skew(·) computes the skewness (with zero being a perfectly symmetric PDF of departures). The inner sum is over geographical bins and the outer sum is over the m selected SSMIS channels j. The aim of the parameter search is to find the set of physical choices v that minimises the cost function J(v).
The form of cost function in Eq. (6) has a number of advantages. A single global mean is of little use, since it allows different 350 situation-dependent biases to be summed together. A set of physical choices leading to compensating systematic errors could score as well as a set with smaller systematic error. As illustrated in Fig. 7, geographical separation helps to implement an approximate regime-based separation -for example the errors over tropical land surfaces are separated from those in frontal cloud in the storm-tracks, and cannot compensate each other. It is also helpful to compute the skewness in geographical bins, as this measure can otherwise be dominated by the large departures in the vicinity of tropical convection. As will be seen later,

355
this is a significant advance on the use of a global skewness by Geer and Baordo (2014) because it helps identify problems with the use of C max over higher latitude land surfaces, an issue that was missed in the earlier work. Normal data assimilation also includes a background term measuring the deviation from prior knowledge, and also serving to regularise the problem.
Section 6.2 discusses the issue of prior knowledge further, and regularisation is not needed as the problem is likely to be well-determined, with 6 unknowns and around 10 6 observations.

360
The VarBC bias correction, b i,j in Eq. (5), is fixed throughout the passive monitoring experiments conducted here. A large part of the modelled bias, especially for SSMIS, is thought to come from the instrument calibration, so it is necessary to include a bias correction. VarBC estimates the bias as a linear combination of a globally constant offset, a number of polynomial terms based on the instrument's scan angle, and a number of model-based predictors including the total column water vapour (TCWV), skin temperature and surface windspeed for surface-sensitive observations over ocean, and various functions of the 365 layer thickness (an average temperature, essentially) for sounding channels. None of the predictors is directly correlated with cloud-related biases, so these do not map directly onto the VarBC bias correction. In some channels the bias correction aliases some of the latitude-dependent biases, such as differences between tropical and midlatitude cloud simulations, but the effect is generally small. Ideally VarBC would be allowed to vary, but as a first approximation it is held constant.

Global parameter search
The first application of the parameter search is to find a single set of physical parameters that would provide the best fit globally, a "one size fits all" approach. Figure 5 helps to visualise the 6-dimensional cost function by showing slices along its dimensions at the control configuration and at the best configuration. At the control configuration, perturbations in most directions would increase cost: for example, all other choices of convective snow particle would make things worse. This indicates that the 375 configuration provided by Geer and Baordo (2014) was already at a reasonably stable local minimum.  Figure 5. Slices through the 6-dimensional cost function around the control configuration (dashed) and around the best configuration (solid, see also Table 4). The dimensions are ordered from least-scattering to most-scattering options. The dots indicate the parameter settings for the control and for the very best (lowest cost) configuration. The best configuration is summarised in Table 4. Here, the shape of the cost function in is very different (Fig. 5). The largest difference is the reversal in the gradient of the cost function with respect to cloud overlap. C av overlap over land is now clearly better than the C max overlap used in the control. Another big difference is that the choice of convective snow particle is now well bounded at the extremes: at low-scattering the ARTS column aggregate has high cost and is clearly inappropriate; at the 380 high-scattering end the ARTS graupel also gives a high cost. The ARTS column is marginally the best, but the ARTS block aggregate or Liu 3-bullet would be a reasonable alternative. A final big difference is the change in the gradient of the cost function with respect to convective snow mixing ratio. At the best configuration, the mixing ratio should be increased by 50%, whereas the control configuration could have been improved by reducing mixing ratio by the same amount.
There are smaller changes in the other dimensions in Fig. 5. For the large-scale snow, a small improvement is found going to 385 the ARTS sector, although the ARTS 6-bullet, ARTS plate aggregate or Liu sector (the control configuration) would produce similar results. The lower-scattering end is not bounded so it might be possible to further improve the simulations with a lessscattering configuration for large-scale snow. The least constrained dimension is the ice cloud microphysics: all configurations produce similar cost, with the gamma 1e4 PSD being only marginally better than other alternatives. Only the F07M PSD with the Liu dendrite shows a visibly increased cost; as previously mentioned this was identified in initial work as having too 390 much scattering. Table 4 summarises the best global configuration: overall this reduces cost to 0.908, compared to 0.925 at the control. This is a small improvement, but as section 5.2 will show, it is possible to do better with situation-dependent choices. in the observations, deep scattering TB depressions are seen at these frequencies only in the very most intense convection (e.g. 400 Zipser et al., 2006). The new search continues to constrain the choices to those that have relatively low scattering at the lower frequencies, while still being able to produce significant scattering at higher frequencies. This illustrates the importance of including the lower frequencies in the cost function. Figure 7 shows maps of the mean and skewness of departures in the geographical bins that make up the cost function (Eq. (6)). The best option improves the mean of binned departures over extratropical land surfaces, particularly in the NH at 405 183±3 GHz. The skewness of binned departures in the control is around +3 in many of these areas; the best option reduces skewness mostly to the range -1 to +1, at both 150 GHz and 183±3 GHz. This is a result of changing from C max to hydrometeorweighted C av cloud overlap over land surfaces. Positive skewness of departures is consistent with the simulations containing much higher scattering TB depression than is observed. Moving to the more physically correct C av approach, the representation of extratropical convection over land is much improved (this will be further illustrated in the next figure). However, some other 410 areas have got slightly worse to compensate, for example there is more negative skewness in the ITCZ over central Africa at 150 GHz and 183±3 GHz and slightly more positive skewness over the Southern Ocean at 183 GHz. Figure 8 illustrates the changes to the simulated SSMIS 150h brightness temperatures when going to the best configuration.
A day was chosen that had a broad variety of cloud and precipitation features across the Atlantic and Europe. Brightness much closer representation of the observations, following the change to the C av cloud overlap over land surfaces. However, the 420 increased brightness temperatures in convection over central Africa take the simulation further from the observations, exposing the problem of missing convection that the C max overlap had previously helped to hide. The best configuration (Table 4) also includes an increased convective mixing ratio and the move to ARTS column to represent convective snow, both of which significantly increase scattering (Table 2, Fig. 7); however this is not enough to fully compensate, ultimately leading to negative skewness over tropical land areas (Fig. 7).

425
Over the ocean, Fig. 8 illustrates a general slight decrease in brightness temperatures, and hence an increase in scattering.
Here, the cloud overlap is not changed but the other four changes in the best configuration (Table 4)     The cost given is that of the chosen configuration but calculated with the standard cost function (half mean, half skewness -Eq. (6))

Sensitivity tests 430
The global best solution from Fig. 5 is fairly stable. Figure 9 shows slices through the cost function around the best solutions that can be obtained with the convective snow mixing ratio fixed to one of the three possible options. Along the dimensions representing cloud overlap, convective snow particle shape, and ice cloud particle shape and PSD, the cost function has a stable form and produces similar solutions in all cases. This suggests that all three choices are reasonably insensitive to the convective snow mixing ratio. However, for convective snow PSD and the large-scale particle, the shape of the cost function, and the best 435 choice, varies depending on the convective snow mixing ratio. In particular the large-scale particle shape has two possible minima -one at the ARTS sector and one at the control (Liu sector) configuration. The latter is found, along with a change to the convective PSD, when the convective snow mixing ratio is halved.

Similar figures can be constructed by holding the choices fixed across any other dimension (not shown). The cost function
can be returned to something resembling its shape at control (as seen in Fig. 5) by retaining C max hydrometeor-weighted 440 cloud overlap over land or by pushing the large-scale snow to either of the two most scattering shapes. Hence at the broadest scales the cost function has two stable areas giving reasonable solutions, one around the control and one around the new best configuration. Varying the convective snow PSD or the ice cloud particle shape has little effect on the shape of the cost function at solution, which remains close to the new best, further confirming its stability. Finally, setting the convective snow particle shape to the extreme possibilities can generate very different cost function shapes, but these lead to low-quality (high cost) 445 configurations. Table 5 explores sensitivity to the choice of cost function itself. The mean and skewness find similar solutions to those found using the average of both (Table 4). On their own they identify some of the reasonable alternative choices, such as the ARTS plate aggregate for large-scale snow, that could already be identified from Fig. 5. Neither the mean or skewnessbased cost function would on its own adjust the convective mixing ratio, consistent with the lack of sensitivity to adjustments 450 in this parameter. A more unexpected aspect of the skewness-based cost function is the choice of the Liu dendrite and the  Table 5 also shows the results using a cost function based on the RMS of the geographically binned departures. This puts most parameters towards the least scattering end of the scale, and in particular it reduces convective snow by 50%. This is not expected to be a valid solution: the double penalty effect typically means that RMS error can be minimised by reducing scattering, essentially by removing cloud from the simulations. This creates a model with an incorrect PDF of brightness temperatures, as measured by the skewness or by examining the PDF directly (not shown). Indeed the cost as measured by the 460 standard cost function is increased substantially to 0.977 at this solution.

Situation-dependent parameter choices
The global "one shape fits all" approach has always had trouble getting good results in the tropics and midlatitudes simultaneously (Geer and Baordo, 2014) and the new global solution remains a trade-off with too little scattering in tropical convection and too much scattering in the storm tracks (Fig. 7h). Table 6 summarises the best configurations when the globe is instead which is not surprising given these scenes dominate the global score numerically, but there is an unexpected reduction of 50% 470 in convective mixing ratio. Compared to Fig. 2, the skewness of background departures is brought much closer to zero at lati- tudes around 60 • over ocean (not shown). This is consistent with the conjecture that particle fall speeds in the flux to mixing ratio conversion are too low (Geer et al., 2017a, and see earlier). This also demonstrates the importance of correctly modelling convection even in the midlatitudes.
The tropical convective ocean solution is very different to the global solution. There is a 50% increase in the convective snow 475 mixing ratio and a move to the Marshall and Palmer (1948) PSD. Both of these changes significantly increase the scattering from convection. Further, the ice cloud is represented with the Liu dendrite, which is the most scattering option. The solution over land is different again, boosting scattering in convective scenes by retaining the C max cloud overlap and again using the Liu dendrite for cloud ice. However, the convective and large-scale snow habits themselves move to slightly less scattering options.
These solutions both make a substantial increase in the amount of scattering coming from convection, again illustrating that the 480 simulations do not provide enough scattering in tropical areas over both land and ocean. These results will be explored further in the discussion (Sect. 6.1). Using the situation-dependent choices from Table 6, the global cost function could be reduced to 0.886, more than doubling the improvements that were made with the best global solution (Table 4). However, the remaining cost is still high, suggesting many further improvements need to be found.

485
One drawback of developing models using parameter estimation is that they need to be updated when other components change. This configuration is otherwise similar to the global best found in Sect. 4, but with slightly different microphysical choices from within the poorly constrained areas of the cost function (Fig. 5). This configuration was important for testing the other developments towards v13. In particular, it was used by  to develop a representation of polarised scattering and to tune the required polarisation ratio parameter. In areas affected by strong scattering, the new approach to polarisation increases brightness temperatures by up to around 7 K in vertically-polarised channels and decreases them by the same amount 495 in horizontally-polarised channels. Since all the high frequency channels of SSMIS are horizontally polarised (≥ 150 GHz; see Table 1), this alters the reference against which the results in the current work were derived.
Two other changes were made to RTTOV-SCATT in the late stages of development for v13. First, there was a bugfix to the downward scattering source terms, leading to changes of up to around 1 K in simulated brightness temperatures in some situations. Second, there was a change to the minimum permitted effective cloud fraction, which gave changes mostly within 500 2 K but occasionally up to around 8 K (for more information on both, see Saunders et al., 2020). Mean changes were an order of magnitude smaller than those examined in typical perturbations in the parameter search (not shown) but, along with the polarisation update, it was important to take these modifications into account in a final re-run of the parameter estimation.
To come up with the final configuration for v13, the three late changes were imposed as additional deltas in Eq.
(2), relying again on the linear additivity of the changes, which is not significantly affected compared to what is illustrated in Fig. 4. None

505
of the changes was treated as optional, since they all make the radiative transfer more physically correct, and the tunable polarisation ratio parameter could not be constrained in the current work. The additional 3 deltas increased the control cost slightly (to 0.932, Table 7). The convective mixing ratio was held constant, given it is undesirable for future developments to build in a tuning factor in this way, and given the lack of sensitivity to this parameter (Sect. 5.2). Table 7 gives the final chosen configuration, again just a slight perturbation on the best configuration found earlier, and still within the uncertainty zone of 510 Figs. 5 and 9.
The final chosen configuration is the best configuration in the parameter estimation, but with an imposed change for cloud ice. The true minimum was found at the control settings for cloud ice (a Mie sphere). Any of the higher-scattering options for cloud ice now generates significantly increased cost, and although the shape of the cost function is similar to Fig. 5f, the gradient is increased (not shown). Hence it is not possible to force the choice of the preferred Heymsfield et al. (2013) PSD any 515 longer. However, the control configuration for ice cloud is undesirable for future sub-mm applications since it uses the outdated Mie sphere approach. A compromise was possible, and the use of the MGD 1e4 PSD with the ARTS column aggregate has been imposed since, after this change, the cost function is only slightly increased. This section discusses whether the parameter estimation has provided physical suggestions that may be generally applicable, or whether the results just reflect tuning to fit the biases of the forecast model used as the reference, and/or the limitations of the radiative transfer model. The biggest improvement was to change the effective cloud fraction over land, which is not surprising as it is the more physically correct approach (Geer et al., 2009a). However, there may be a number of more novel suggestions that could be generally applicable. The first of these relates to convective cloud in the tropics. When the search was confined to land areas (Table 6) it retained the C max cloud fraction, despite the fact it is less physical. It is hence worth asking if the specified 5% convective fraction is too small to represent intense tropical convection systems, noting that even in cloud-resolving models this is an uncertain parameter (e.g. Varble et al., 2014). Radar cross-sections can show vertically spreading high-reflectivity areas in the upper parts of these clouds (e.g. Hong et al., 2005, their Fig. 1). Hence a hypothesis, based also on results from cloud-resolving models, is that the frozen upper parts of the convective cores become more widely 530 spread and merge together in a gradual transition to anvil cloud (Geer et al., 2017a). The basic representation of convection in the IFS and RTTOV-SCATT, as a narrow vertical core and a thin high anvil, may therefore be an overly severe truncation of reality.
Further justification for the hypothesis on the structure of intense tropical convection is seen over tropical ocean surfaces in the situation-dependent search ( Table 6). The results over ocean are less affected by model bias: here the more physically-535 correct C av cloud fraction is used, and the all-sky infrared biases suggest only a slight over-estimation of convective cirrus in frequency or coverage (Geer et al., 2019). But the parameter estimation still boosts scattering by increasing the convective mixing ratio, by increasing the size of the convective particles using the MP48 PSD, and by making the cirrus more scattering too (using the F07M PSD and the dendrite snow particle). All of these could suggest the need to represent either broader convective cores in the frozen parts of the cloud, or nearly equivalently, deeper convective anvils, containing higher mixing 540 ratios of larger and more scattering particles.
Another interesting point is the level of constraint on the convective snow particles. At the most scattering end, the ARTS graupel particle is decisively rejected. This is always a bad option, independent of the choice of cloud overlap over land ( Fig. 5) and independent of convective mixing ratio (Fig. 9), because it generates too much scattering at lower frequencies.
However, the global search finds solutions (such as the Liu 3-bullet or ARTS column) that strongly increase scattering at 545 150 GHz and 183 GHz, and by nearly as much as the ARTS graupel (Fig. 3). The parameter estimation prefers a spectral signature of scattering that is strong around 183 GHz but weak around 50 GHz. This eliminates a wide class of dense particles  and similarly rejects the least-scattering particles, particularly the ARTS column aggregate, which strongly reduces scattering around 183 GHz (Fig. 3). This spectral signature is hence strong enough to be mostly independent of both the sub-grid structure and the mixing ratio of the convective cloud.

550
The detectable influence of a particle model in this work is its scattering signature, as a function of frequency and geographical location, not directly its physical properties. The scattering signature comes from the bulk scattering properties that arise from the combination of particle size distribution and the single-particle scattering properties. In the approach used here the chosen particle model affects both the single-particle scattering signature and the PSD, since the particle choice also defines the mass-size relation. Changing the mass-size relation changes the shape of the PSD (e.g. Eriksson et al., 2015;Geer et al., 2021).

555
In the mass-size relation the ARTS graupel has an exponent of 2.96; the Liu 3-bullet and the ARTS column have exponents 2.37 and 2.05 . This suggests the scattering signature may be as sensitive to the exponent in the mass-size relation as to the actual particle shape. This would help explain how a habit that is physically unlikely in convection, such as the ARTS column, has ended up being chosen to represent convective snow.
all the PSDs that were explored here as representations of ice cloud (McFarquhar and Heymsfield, 1997;Field et al., 2007;Heymsfield et al., 2013). Note the rejection was much clearer in the final search for RTTOV v13 in Sect. 5.4 than in the main search, presumably because the scattering generated by ice cloud had been boosted by the parametrisation for polarised scattering. The PSDs from the literature were trialled with the very least-scattering particles available to RTTOV-SCATT, i.e.
the ARTS column aggregate and the Evans snow aggregate (see Geer et al., 2021). Even so, the observations could only be 565 matched by creating a modified gamma distribution (MGD) with the key property that it generated orders of magnitude fewer centimetre-sized particles than the PSDs from the literature (see Geer et al., 2021). Therefore, although frequencies around 183 GHz are not expected to have strong sensitivity to smaller cloud ice particles, they still impose a constraint on the number of large particles in the PSD.
It is difficult to compare directly to other studies. A main reason is that the characterisation of frozen particles, into large-scale 570 snow, convective snow, and cloud ice, is specific to certain types of global model. It is hard to map onto other representations of ice particle, particularly those in which no forecast model is available to categorise the hydrometeor type (e.g. Ekelund et al., 2020). Further, many studies do not attempt to fit such a range of frequencies or such a complete sample of observations (either in frequency or location) so they are able to find solutions that would not work here. Finally it is worth restating that the results are valid on the 80 km by 80 km smoothing/superobbing grid used for all-sky assimilation in the IFS; hence results may 575 be different from a much finer-scale cloud-resolving viewpoint. However the following examples based on regional models illustrate that there is plenty of common ground between them and the conclusions from the current parameter estimation. Guerbette et al. (2016) used RTTOV-SCATT in combination with a tropical local area model, to simulate all-sky microwave observations in the tropics around 183 GHz. Their study applies to snow, but with a focus on tropical cyclones, so this is interpreted as mapping onto the convective snow category. They tried all the Liu (2008) particles available to RTTOV-SCATT, 580 finding the best fit to observations using the Liu block hex column and the F07T particle size distribution. However they did not include any lower frequencies to help constrain their particle choice. If they had done, the choice of the Liu block hex column would not have worked. At frequencies of 183 GHz and below, this particle is roughly the highest scattering of any particle available in RTTOV-SCATT. At 50 GHz, for example, it is much more scattering even than the ARTS graupel . However, what is consistent in their results compared to the current study is the need to generate as much scattering as 585 possible in tropical convection at 183 GHz.
Fox (2020) performed a closure study between a regional forecast model and airborne observations over northern Europe covering the frequency range 89 GHz to 874 GHz, using the ARTS model for radiative transfer. The regional forecast model had a hydrometeor category known as "cloud ice" which does not distinguish between suspended and falling particles, e.g.
between ice cloud and snow. Using the Field et al. (2007) tropical PSD, the ARTS column aggregate gave the best results. The 590 study also showed that the ARTS block aggregate produced too much scattering. Although there is no direct equivalent to the combined ice category, the choice of mid-to-less-scattering particles is still consistent with the parameter estimation finding less-scattering particles in the IFS snow and ice cloud categories (along with a less-scattering PSD in the ice cloud category).
Many other studies have explored the best representation of frozen particles using non-spherical particles but there is no space for an exhaustive comparison (among them Kulie et al., 2010;Eriksson et al., 2015;Duruisseau et al., 2019;Fox et al., 595 2019; Ekelund et al., 2020). This discussion has at least shown it is possible, with careful consideration, to find consistency with other studies. However, it is important to consider the precise hydrometeor definitions, the cloud situations that are targeted, and even the level of constraint (whether one frequency or a wider spectrum, whether a globally representative study or one focusing on a smaller area). Since many of the results of this study can also be explained in a physically plausible way, these conclusions probably can be applied more widely than just to the modelling framework in which they were derived. 600

Parameter estimation methodology
This section explores the current parameter estimation in the context of wider efforts to learn models from data. Data assimilation, along with machine learning, can be derived from Bayesian principles (e.g. Bocquet et al., 2020) where knowledge is represented quantitatively using probability distribution functions (PDFs): this gives a fundamental framework for how observations can be combined with prior knowledge to provide improved knowledge of physical states and processes. The process 605 of learning model parameters through the Bayesian framework is known as parameter estimation. In weather and ocean forecasting, the history of parameter estimation is long but sparse (e.g. Yu and O'Brien, 1991;Dee and Da Silva, 1999;Ruiz et al., 2013). Cloud and precipitation parameter estimation has been tried in simpler contexts (e.g. Norris and Da Silva, 2007;Posselt and Vukicevic, 2010;Posselt and Bishop, 2012) and within more complete atmospheric models (Ollinaho et al., 2013;Posselt, 2016;Ruckstuhl and Janjić, 2020;Kotsuki et al., 2020). However it has not yet reached a level of development to have been 610 used in operational weather forecasts. The ultimate goal has been demonstrated in an idealised context by Zhu and Navon (1999): a combined state and parameter estimation where the improved model parameters benefit the quality of the forecasts, even at longer forecast ranges when the impact of improved initial conditions becomes small.
Parameter estimation remains challenging. Parameters may not be well-constrained by the available observations (an issue known as identifiability, see e.g. Dee and Da Silva, 1999). When there are model components that are not subject to parameter 615 estimation but have compensating errors, the estimated parameters will not be fully realistic (e.g. Ruiz and Pulido, 2015).
Especially in cloud and precipitation, the uncertainty PDFs are highly non-Gaussian and quite often multimodal (Posselt and Bishop, 2012). Equivalently, cost functions have multiple minima, making for a difficult optimisation problem. The complexity of the cost function in the current parameter estimation is clear from Figs. 5 and 9. In the search for globally constant parameters, the parameter estimation had to jump from the existing 'control' configuration to an area of parameter space that 620 gave substantially better fits to observations. Hence this work further confirms the problem of multiple minima; this would make it hard to apply typical data assimilation frameworks, which use gradient descent methods that could get stuck in a local minimum.
When the cost function is complex and has multiple minima, it important to search the parameter space broadly. The current work has explored the whole of the search space. In larger systems this would become infeasible due to the curse of dimen-625 sionality. Here, the space has been kept small (a total of 3456 discrete parameter combinations) and an assumption of linear addition of perturbations allowed the entire space to be approximated on the basis of only 22 simulation experiments, one for each parameter setting, plus the control. The linear assumption has limited validity in the presence of cloud and precipitation processes that are known to be nonlinear, but with careful design and within the scope of the current study it gave reasonable results, to a maximum of only 20% error in the reconstructed brightness temperature in one example. The method could be 630 extended to much larger search spaces, allowing parameters from the forecast model to be included; however the exhaustive evaluation of the cost function (Eqs. (4), (5) and (6)) would itself become a burden and the number of evaluations would have to be reduced by using a more efficient method such as Markov Chain Monte-Carlo or other known approaches (e.g. Gelman et al., 2013). If the linear addition assumption breaks down (which can be tested a-posteriori) then the underlying nonlinear experiments could be re-run iteratively to help improve the solution, a litle bit like the incremental 4D-Var approach (Courtier 635 et al., 1994).
Applying Bayesian methods in the form of data assimilation, it is normal to use a Gaussian error assumption, or equivalently a squared-error cost function (see e.g. Geer, 2021). This seems to be unsuitable for cloud and precipitation parameter estimation and hence many different cost metrics (or equivalently, uncertainty representations) have been proposed, including both higher and lower order statistical moments (e.g. means, skewness and kurtosis), threat scores, correlations and PDF di-640 vergence (Geer and Baordo, 2014;Duan et al., 2017;Schneider et al., 2017;Ruckstuhl and Janjić, 2018). In the current work, the geographically-binned mean and skewness have been used, and the RMS error is shown to produce a false "minimum scattering" solution, in other words trying to remove cloud from the simulations. To avoid this, the mean targets broad systematic errors of the kind that can affect the model climate. The skewness provides higher sensitivity to errors in extreme cloud and precipitation situations and is a more feasible alternative to computing the full divergence between simulated and observed 645 PDFs (Geer and Baordo, 2014). The statistics are aggregated across latitude-longitude bins and across a range of microwave frequencies from 19 GHz to 183 GHz. Spatial binning provides regime-sensitivity and helps make the parameters relating to large-scale snow, convection and ice cloud more separately identifiable. The spectral dimension helps identify large-scale and convective snow particle shapes through their distinct spectral signatures.
It is also interesting to compare the current work to the ideal Bayesian approach (e.g. Gelman et al., 2013;Posselt and Bishop, 650 2012; Posselt, 2016;Geer, 2021) even if that is likely to be infeasible in a high-dimensional system. The first point of difference is the treatment of prior knowledge, which in a Bayesian approach is represented by a PDF, and in variational methods this corresponds to the background term in the cost function. Although substantial expert knowledge exists, for example the way ice particle habits vary with supersaturation and with temperature (e.g. Bailey and Hallett, 2009), this is hard to transfer into a quantitative form like a PDF or cost function. Further, as discussed by Geer and Baordo (2014), the aircraft data on which 655 much of this knowledge is based can be limited to particular areas of the globe and certainly it is limited to particular types of cloud; this means current knowledge may not be universally appropriate when trying to fit observations on a global basis.
Here, the lack of an explicit prior means that every position in parameter space could be equally likely. This corresponds to the common Bayesian approach of "uniform priors", which is a reasonable starting point when prior knowledge is hard to quantify. However, the choice of discrete options presented to the parameter search does encode substantial prior knowledge.
For example, completely unrealistic particle and PSD choices have already been discarded. Where significant uncertainty exists, such as in the modelling of convective ice cloud, there are three dimensions of variability available and a broader range of options.
This severe restriction of the search space is perhaps the biggest difference compared to the ideal approach. To find the true best model configuration, every parameter should be infinitely variable and every model component should be open to 665 improvement, for example by changing the equations on which they are based. However, if the need for a new search point becomes obvious from the cost function, for example a 25% increase in convective mixing ratio, it is still possible to add one and hence to refine the solution further. It is also possible to add new dimensions when they become necessary. Other aspects also map onto the broader field of Bayesian methods: the linear assumption is being used in a similar way to which a statistical emulator is often used to speed up a parameter estimation (e.g. Duan et al., 2017). Finally, although the current study is based 670 on an ad-hoc cost function, that cost function implies a PDF that could give a clue to the shape of the true, and non-Gaussian, likelihood PDF that could be employed in a more formally Bayesian approach.

Conclusion
This work presented a parameter search across 6 different microphysical and macrophysical assumptions associated with frozen cloud and precipitation in a scattering radiative transfer model. The search used a cost function based on the misfit (or departure) 675 between observations and simulations based on geophysical fields from the Integrated Forecasting System (IFS). Precisely, this was based on the geographically binned mean and skewness of departures, in equal weights. A 10 day period of Special Sensor Microwave Imager Sounder (SSMIS) observations was used, providing around 10 6 observations with near global sampling (ice-covered polar regions and high orography were the main exclusions), and frequency coverage from 19 GHz to 190 GHz.
The practical output was the recommended default configuration for version 13 of the Radiative Transfer for TOVS microwave 680 scattering module (RTTOV-SCATT), released in late 2020 (Saunders et al., 2020). A main goal for the future is to optimise the parameters of the moist physics of the forecast model alongside those in the radiative transfer, to reduce the effect of compensating model errors, but for the moment these were held constant.
The new RTTOV-SCATT defaults (Table 7) are intended to provide a reasonable global "one shape fits all" configuration for microwave and sub-mm scattering radiative transfer. New microphysical settings are provided for three frozen hydrometeor 685 types. Cloud ice is moved from a Mie sphere to an ARTS large column aggregate, in order to provide more accurate simulations for the future ICI instrument, which will operate at frequencies high enough to be sensitive to the microstructure of cloud ice.
However, standard ice PSDs from the literature have been rejected as they produce too many large (cm-sized) particles, which causes excessive scattering even around 183 GHz. A newly proposed exponential PSD is used to avoid this. The problem is likely to be the representation of separate snow and ice cloud categories, which is typical in global models like the IFS, whereas 690 the PSDs based on in-situ observations may represent both. The large-scale snow hydrometeor is moved from the Liu (2008) sector snowflake to a more physically representative ARTS large plate aggregate, which again should give better simulations at high frequencies. Finally, the new convective snow hydrometeor uses the Field et al. (2007) tropical PSD with the ARTS column. Though the habit is physically unlikely, it provides the best fit to the observed scattering signature of minimal scattering at 50 GHz and maximum scattering at 183 GHz. Apart from this one aspect, the settings can be justified physically and there is some consistency with other studies.
One issue that was clearly linked to a model bias was the move to the hydrometeor-weighted effective cloud fraction (Geer et al., 2009a) over land as well as ocean surfaces. This has always been the RTTOV-SCATT default, but Geer and Baordo (2014) used the maximum cloud fraction to boost scattering in tropical convection over land. In the new parameter estimation, the use of geographically-binned skewness statistics gave a more appropriate representation of the extratropics in the cost 700 function. This highlighted that the maximum cloud overlap causes large errors in convection outside the tropics. Going to the hydrometeor-weighted overlap reduced simulation errors in these conditions by up to around 50 K, and was the single largest improvement compared to the previous (control) configuration.
The new global configuration reduces the cost function from 0.925 to 0.908, highlighting that many errors remain. A possible way forward is to use situation-dependent parameter choices. Here an example was to split the data into land, tropical ocean, and 705 the rest of the ocean, and performing three separate parameter estimations. This allowed the global cost to be further reduced to 0.886. However, some of the chosen settings were unexpected, such as the retention of the less-physical maximum cloud fraction for land areas, and the choice of the most scattering available representation for cloud ice over land and tropical ocean.
These highlighted potential issues with the representation of tropical convection. A hypothesis is that the modelling should better represent the broadening of convective cores in the frozen upper parts of tropical convection, or nearly equivalently, 710 a deeper anvil composed of larger and more scattering particles. Hence it is important to find better ways to represent the structure of deep convection in the tropics.
Compared to the earlier work of Geer and Baordo (2014), the parameter search methodology is a step forward in two ways: a total of 6 parameters were adjusted simultaneously, and the quality of the parameter choices was measured using a single cost function. This latter makes the method a form of parameter estimation, which can be compared to a more ideal 715 Bayesian approach (e.g. Posselt and Bishop, 2012;Posselt, 2016;Geer, 2021) in which the PDF of prior knowledge is updated based on observations, given the knowledge of the likelihood (conditional) PDF of those observations. Prior knowledge is not explicitly represented in the current work but it is encoded in the range and choice of parameters for optimisation. The likelihood PDF is implicitly represented in the choice of cost function. The search space has been kept deliberately small (3456 possible combinations) so that it could be searched in its entirety. Even so, because the experimentation was done in a weather 720 forecasting system, it was too costly to explore each combination exactly. Instead an approximate method was used, assuming that parameter perturbations could be tested separately (in 22 experiments) and then linearly added to explore the whole space.
The many links to the Bayesian approach should help improve the method in future, taking further ideas from a wide body of work on Bayesian inverse methods.
An important question in parameter estimation is whether the parameters are identifiable. Here, it has been shown that 725 multiple physical parameters of large-scale snow, convective snow, and cloud ice can be separately identified. The particle model for convective snow is particularly well constrained, being mostly insensitive to variations in the mass mixing ratio.
The constraint on the particle model comes from the strong frequency variations in the bulk scattering properties of frozen particles. This "spectral signature" of scattering is generated by the combination of the particle size distribution, the mass-size signature is constrained by the inclusion of frequencies from 19 GHz to 190 GHz in the study. Particularly, the inclusion of lower microwave frequencies (e.g. 50 GHz) helps eliminate many proposed particle models that generate far too much scattering at these frequencies, such as the ARTS gem graupel. However it is important to be clear that the particle habit itself is not being precisely identified, but rather the full "hydrometeor model" which is controlled by the choice of particle habit in the current framework but includes also an assumed mass-size relation ). An important future step will be to 735 better separate the PSD shape, the mass-size relation, and the particle habits in the parameter estimation.
There is an aspiration to do more cloud and precipitation parameter estimation within the data assimilation systems used at weather forecasting centres. The current work illustrates some difficulties: First, data assimilation usually makes the assumption of Gaussian errors, but the suggestion here is that they are far from Gaussian; Second, the cost function has multiple minima, making it hard for conventional techniques to find the right solution and requiring a more complete search of the parameter 740 space than is normally possible; Third, it is difficult to quantify prior knowledge in the parameters. However, these issues might be addressed by including an ever wider selection of cloud and precipitation-sensitive observations in the cost function, a strategy discussed by Geer et al. (2017a) as "microphysical and macrophysical closure". The Ice Cloud Imager (ICI), when it is launched in 2024, will add new frequencies from 243 GHz to 664 GHz which will better constrain smaller frozen particles.
Existing passive radiance observations in the infrared and visible parts of the spectrum will add much information on cloud 745 fraction and the microphysics of ice and water cloud. Active measurements from cloud radar and lidar will be particularly useful in terms of the 3D structure of clouds, for example bringing further constraints on the convective core fraction and anvil coverage that has caused uncertainty in the current work. It would also be extremely useful to include field site and campaign data on clouds and precipitation directly in the parameter estimation, rather than using indirect representations such as PSDs from the literature. This additional information could make any type of parameter estimation easier: First, by making the 750 cost function less complex (by reducing ambiguity between parameters, and by eliminating physically implausible choices); Second, by reducing the need to quantify prior knowledge. Hence it is worth continuing to develop cloud and precipitation parameter estimation both inside and outside existing data assimilation systems, which will be greatly helped by an ever wider set of observations from space and ground-based instrumentation.
Code and data availability. This work used development versions of RTTOV leading up to v13.0, which is available to registered users from https://nwp-saf.eumetsat.int/. Data assimilation experiments were adaptations of version 46r1 of the ECMWF IFS which is accessed through the national hydrometeorological services of its member states. The large amount of data generated by these experiments is not viable for long-term archiving but publicly available and curated versions are available through the ERA-5 project.
Author contributions. The sole author did the work